李亚飞，周懿，胡钺，高政

非牛顿流体在不规则管道中的流动特性研究＊

李亚飞，周懿，胡钺，高政

（船舶动力工程技术交通行业重点实验室，湖北 武汉 430063；武汉理工大学 能源与动力工程学院，湖北 武汉 430063）

非牛顿流体被大量应用于工业生产中，如利用剪切稀化流体应用的酱类食品生产运输、高分子聚合物在工业生产中的管道运输及钻井用洗井液完井液等。相比传统的牛顿流体，关于非牛顿流体在管道内的流动的研究还有很大发展空间，因此对其机理的研究具有重要的现实意义。借助多物理场仿真软件COMSOL Multiphysics，对非牛顿流体在不规则管道内的流动现象进行了数值模拟研究，分析了非牛顿流体在管道内的速度场分布、剪切应变速率分布以及黏度分布情况。结果表明，在同一不规则管道中，管内流动的非牛顿Carreau流体的剪切应变速率与黏度关系符合其本构方程特性，相比牛顿流体表现出了明显较差的流动性，且速度场的分布一直保持一致；实验组非牛顿流体的剪切应变速率量级比对照组牛顿流体小很多，解释了两种流体在速度场分布的差异性；验证了Carreau流体构型中剪切应变速率与黏度的负相关关系。

非牛顿流体；管道；流动特性；幂律流体

非牛顿流体是剪切力与剪切应变率之间不是线性关系的流体。非牛顿流体的应用十分普遍，涉及采用剪切增稠液体的装甲、食品工厂设计果酱类食品的运输加工及高分子聚合物在化工生产中的运输等诸多方面。这类非牛顿流体在工业生产方面的应用研究，直接影响工业生产效率及国家经济，因此对非牛顿流体在管道中的流动机理具有较大的研究价值。

对非牛顿流体的研究方法包括实验法、解析解法与数值解法。实验法直接，但成本高且实验普遍性不高。解析解法理论上为最理想的研究方法，通过建立合适的微分方程，使用纯数学方法得出方程的精确解；但对于非牛顿流体复杂的流动情况来说，求解的难度过大。数值解法则是应用计算机将物理场离散化，之后将流体微分方程组转化为代数方程并求出各个节点上的参数值，属于一种近似解法。由于数值解法容易获得且能保证足够的求解精度，其已经成为研究非牛顿流体问题最为常用的方法。

本文选取一种注塑喷嘴模型作为研究对象，以COMSOL Multiphysics这一软件作为计算平台，采用有限元法对非牛顿流体流经不规则管道时表现出的流动特性进行数值模拟研究。在不规则腔体内，对照研究非牛顿Carreau流体与幂律牛顿流体的流动特性，讨论非牛顿流体的剪切应变速率与黏度对于速度大小与速度场分布的影响。

1 管内流动数学模型

1.1 控制方程组

一切的流体流动过程，都以三个基本的物理学原理为基础：质量守恒定律，牛顿第二定律与能量守恒定律。将这些物理学用于构建流动模型，将会导出一组方程，即连续性方程、动量方程与能量方程。这些方程是上述物理学原理的数学描述，本文不讨论传热，因此不引入能量守恒方程：

质量守恒定律（连续性方程）：

式（1）表示瞬态三维可压流体流动的连续性方程，本文所分析的流体流动处于稳态且不可压缩，密度不会随着时间的变化而改变，因此流体流动的数学描述为：

动量守恒方程：

式（3）（4）（5）是对于任何流体都成立的动量守恒方程，是微元体内流体动量对于事件的变化率等于外界作用于该微元体上的各种力的和，简称动量方程，也称纳维斯托克斯（N-S）方程。

1.2 本构方程

本构方程是反映物料宏观性质的数学模型，又被称为流变状态方程或是流变方程。在流变学中，本构方程是在某些假定条件下，对流体或弹性体的材料力学行为的数学描述，可用来区分流体类型。本构方程与连续性方程、运动方程一起构成封闭的方程组，用于求解流体的流动特性。

在本文仿真中使用的幂律流体的本构方程式如下：

式（6）中：为黏稠系数，表示物料的黏稠程度，Pa·sn；为幂律流变指数（简称幂律指数），为无量纲量，表示非牛顿流体的流动特性偏离牛顿流体的程度（=1时为牛顿 流体）。

1.3 非牛顿流体管道流动模型

剪切应变速率描述的是流体的剪切流动，定义为单位时间的剪切应变变化：

式（7）中：为剪切应变。

值得注意的是，剪切应变速率常与速度梯度混淆。实际上二者是不同的概念。速度梯度是流体的速度对空间坐标的导数，用/来表示。在数学上，二者的数值有时相等，这是因为一般速度梯度符合：

但二者的物理意义并不相同，且数值上有时并不相等（如流体在同轴圆筒之间的流动，此时有角速度的影响）。

2 物理模型及计算方法

在规则管道内，流体的速度场分布较为均匀，各流动特性参数的变化梯度较小。本文在这一部分选取一不规则管道腔体进行管道流动特性仿真，以更清楚地展示流体的流动特性。另选取适当的非牛顿流体与牛顿流体进行变量控制，构成实验组与对照组。通过对比两种流体的仿真结果，说明非牛顿流体不同于牛顿流体的流动特性。

2.1 物理模型

2.1.1 几何模型

本组仿真的计算区域如图1所示，为一个呈圆柱对称的不规则弯曲管道腔体。流体由腔体的圆形入口处流入，中途流经一截面当量直径较小的喉口，最后由腔体的环形出口处流出。腔体的总长度为42 mm，圆形入口直径为18 mm，环形出口宽度为12 mm，其余具体几何尺寸如图2所示。A—B为管道进口段；F—G为管道出口段；B—C—D—E—F 与 G—H—I—J—A为管道壁面。使用此不规则腔体进行数值模拟的原因：在腔体喉口处由于其的截面当量直径相对较小，流体经过时会出现较大的压力梯度，进而显著地影响速度、剪切率与动力黏度。所以在喉口当量直径最窄处选取一个截面，计算分析两种流体各项流动参数在这个截面上的分布，以更精确地说明流体的流动特性。截面具体位置如图3所示。喉口处截面的位置与图2中R12与R3圆心连线重合。

图1 不规则管道腔体示意图

图2 不规则管道腔体的具体尺寸（单位：mm）

图3 喉口处截面的位置（单位：mm）

2.1.2 流体性质

2.1.2.1 非牛顿流体性质（实验组）

本组仿真计算所使用的非牛顿流体为线性聚苯乙烯在1-氯苯中的溶液（A solution of linear polystyrene in 1-chloronaphthalene），其被广泛用于塑料制品制造中。该溶液在图1中不规则腔体内的流动过程，可视为对塑料制品生产中注塑工序的溶液流经注塑喷嘴工段过程的模拟。该聚苯乙烯溶液属于非牛顿Carreau流体，其黏度eff与流体剪切应变速率相关，具体关系式如下：

式（9）中：0为零剪切应变速率时的黏度；inf为无限大剪切应变速率时的黏度；为松弛时间，s；为模型指数。

本溶液的相关非牛顿Carreau流体模型参数值与仿真需要的其他物性如表1所示。对照组牛顿流体的物性参数值如表2所示。

2.1.2.2 牛顿流体性质（对照组）

根据控制单一变量的原则，本组仿真所使用的对照组牛顿流体只进行黏度相关系数的调整，其他物性参数与实验组非牛顿流体保持一致。对照组牛顿流体使用幂律流体模型，幂律指数=1。当幂律指数=1时，流体的黏度eff实际上固定为1 Pa·s。

表1 线性聚苯乙烯在1-氯苯中的溶液的物性参数值

参数名参数值 零剪切应变速率时黏度μ00 无限大剪切应变速率时黏度μinf /（Pa·s）166 松弛时间λ/s1.73×10-2 Carreau模型指数n 0.583 密度ρ/（kg/m3）450

表2 对照组牛顿流体的物性参数值

参数名参数值 稠度指数m1 幂律指数n1 最低剪切应变速率/s-10.01 密度ρ/（kg/m3）450

2.2 边界条件

进口边界条件：进口为图2中A—B段。设定进口压力in，取值从10 kPa开始，以40 kPa的梯度递增至210 kPa，模拟不同进口压力下实验组与对照组流体的流动特性。

出口边界条件：出口为图2中F—G段，设定出口压力为0。

壁面边界条件：壁面为图2中B—C—D—E—F与G—H—I—J—A段。设定无滑移边界条件，即壁面处的流体速度为0。

无黏性应力条件：·[（▽+（▽）T）]。

其他：流体不可压缩；所有管道区域内实验组非牛顿流体处于层流状态，进口压力较大时，对照组牛顿流体会出现湍流。

2.3 求解器介绍

本次模拟使用的COMSOL Multiphysics是一款功能强大的多物理场仿真软件，可以对多个领域的物理过程进行模拟计算。

本文研究的是非牛顿流体在管道内的流动特性，因此选用的是COMSOLMultiphysics中的Laminar flow求解器模块。这个模块内置了许多流动的基本物理参数，可以在这些参数的基础上结合前文所述的数学模型定义更多仿真所需要的变量。

2.4 数值计算方法

本文涉及的数值计算使用的是有限元分析法（Finite element method，FEM）。它的基本思路为：一个物体或系统被分解为由多个相互联结、简单、独立的点组成的几何模型。在这种方法中这些独立的点的数量是有限的，因此被称为有限元。由实际的物理模型中推导出来的平衡方程式被使用到每个点上，由此产生了一个方程组，这个方程组可以用线性代数的方法来求解。

3 仿真结果分析

3.1 两种流体的速度场对比

当进口压力分别为10 kPa、90 kPa、170 kPa时，从实验数据分析可知，随着进口压力的增大，牛顿流体的速度场分布情况发生了较为明显的变化。进口压力为10 kPa时，流体速度在腔体入口段达到最大值，喉口处的中心速度也达到了局部峰值，喉口段之后的流体速度渐缓。进口压力增大时，腔体入口段的速度场分布情况差异不大，但喉口处的局部速度峰值与入口段的速度最大值之间的梯度差异会逐渐减小，最终在喉口处达到全局速度最大值。在进口压力为90 kPa时，速度入口段的速度最大值约为 57.669 m/s，喉口处速度峰值则在52.111 m/s附近；进口压力为170 kPa时，速度入口段的速度最大值约为82.369 m/s，喉口处速度峰值则已经达到85.698 m/s。此外，喉口后部的扩散段速度场分布也随着进口压力增大而出现较为明显的差异。当进口压力为90 kPa时，喉口处的速度峰值区域已经开始变得狭长，峰值点更接近下方壁面。峰值区域狭长化和峰值点向下方壁面移动的程度随着进口压力进一步上升而强化，这一特点从进口压力达到170 kPa之后就非常明显。另一方面，流场内流体速度的绝对值也随着入口压力的增加出现明显的上升。速度峰值由10 kPa入口压力时的2.746 m/s上升至90 kPa入口压力时的57.669 m/s，在170 kPa时则是82.369 m/s。

不同进口压力下非牛顿流体的速度场分布如图4所示。从图4可以看出，在进口压力为10 kPa时，非牛顿流体的速度场分布与牛顿流体的分布情况相似；但随着进口压力的增加，非牛顿流体的速度场分布并不会出现明显变化，这与牛顿流体不同。此外，非牛顿流体的速度绝对值也远小于牛顿流体，进口压力为10 kPa时速度峰值约为0.015 m/s，进口压力为90 kPa时速度峰值约为0.157 m/s，进口压力为 170 kPa时速度峰值约为0.346 m/s。很明显，实验组的非牛顿流体流动性相比对照组的牛顿流体要差得多。

造成这种流动性差异的原因可以通过比较两种流体的剪切应变速率和黏性来得出。非牛顿流体剪切应变速率的分布情况如图5所示。

从实验数据中可得牛顿流体的剪切应变速率值远大于非牛顿流体的（3个数量级的差异），说明单位时间内牛顿流体的剪切应变变化量更多，表现为较高的流速；非牛顿流体则相反。

另外，如前文所述，对照组牛顿流体由于幂律指数=1，所以流体的黏度大小固定为1 Pa·s；而实验组非牛顿流体的黏度由Carreau模型决定，为一个与剪切应变速率呈负相关的变量，其具体大小与分布如图6所示。可以清楚地得出，实验组非牛顿流体的黏度远远大于对照组牛顿流体的黏度，大部分区域黏度值达140 Pa·s以上，这也可以解释非牛顿流体的流动性为何远不如牛顿流体。

图4 不同进口压力下非牛顿流体的速度场分布

图 5 不同进口压力下非牛顿流体的剪切应变速率分布

图6 不同进口压力下非牛顿流体的黏度分布

3.2 喉口截面处的流动特性对比分析

从之前对两种流体速度场的分析可以得知，速度、黏度等变量在流动区域的喉口处均会出现比较明显的梯度变化。利用流动腔体圆柱对称的特性，通过二维函数图像即可确定各项流动参数在喉口处截面上的分布。

两种流体流速沿喉口截面径向的分布曲线如图7和图8所示。

图7与图8所反映的速度大小变化情况，与前文图5反映的喉口处速度场分布相符。随着进口压力的增大，喉口处的流体流速也相应增加。图7中各条速度线的峰值，对应的是实验中牛顿流体的喉口处局部速度峰值；峰值位置的横坐标随进口压力增大而向左横移，对应的是实验时牛顿流体局部速度峰值点区域随进口压力增大的下移。图8中各速度线峰值横坐标则几乎没有出现横移的现象，也符合图4中非牛顿流体不同进口压力下速度场分布大致相同的情况。

从式（9）与式（10）可以看出，实验组非牛顿流体的速度分布与剪切应变速率分布呈正相关，剪切应变速率的分布又与黏度分布呈负相关。为了验证这两组关系，绘制实验组流体的剪切应变速率与黏度沿截面径向的分布曲线，如图9、图10所示。

从图9、图10可以看出，随着进口压力增大，实验组非牛顿流体在截面上的剪切应变速率曲线逐渐上凹，黏度曲线逐渐下凸。速度曲线与剪切应变速率曲线的变化方向相同，剪切应变速率曲线又与黏度曲线变化的趋势相反，这就验证了由式（9）与表1数据所表示的剪切应变速率与黏度之间的负相关关系。

图7 牛顿流体在截面上的速度分布

图8 非牛顿流体在截面上的速度分布

图9 非牛顿流体在截面上的剪切应变速率分布

图10 非牛顿流体在截面上的黏度分布

4 结论

本文通过 COMSOL Multiphysics平台，使用有限元法与局部优化网格技术，建立了非牛顿流体在不规则管道腔体内的流动模型。通过以上数值模型，分析了非牛顿流体在管道内的速度分布、剪切应变速率分布及黏度分布，并讨论了流体流动特性参数之间的关系。

具体总结如下：①随着进口压力的增大，对照组牛顿流体的流速不断提高，速度场的分布中出现狭长状局部峰值区域；实验组非牛顿流体流速一直非常小，表现出很差的流动性，且速度场的分布情况一直保持一致。②在高速管道流动中，非牛顿Carreau流体的剪切应变速率远小于牛顿流体，进而导致其黏度比牛顿流体的黏度大得多，表现为管道内极差的流动性与不随进口压力变化而改变的速度场分布。③通过对喉口处非牛顿流体的剪切应变速率与黏度分布曲线的分析，验证了Carreau流体构型中剪切应变速率与黏度的负相关关系。

O373

A

10.15913/j.cnki.kjycx.2020.08.009

2095－6835（2020）08－0024－05

国家自然科学基金青年科学基金项目（编号：51806157）

李亚飞（1998—），男，湖北武汉人，本科在读。

〔编辑：严丽琴〕