非牛顿流体 — CAE用语解释

分类：词汇表 | 2026-01-15

CAE visualization for non newtonian fluid - technical simulation diagram

非牛顿流体

🙋

老师，CFD中如何处理非牛顿流体？

非牛顿流体的理论基础

基本概念和控制方程

🙋

非牛顿流体与牛顿流体在根本上有什么不同？仅仅是粘度不恒定吗？

🎓

根本区别在于剪切应力与剪切速率的关系是线性还是非线性。对于牛顿流体，粘度μ是常数，关系式为

$$ \tau = \mu \dot{\gamma} $$

。而对于非牛顿流体，表观粘度η是剪切速率

$$ \dot{\gamma} $$

的函数，表示为

$$ \tau = \eta(\dot{\gamma}) \dot{\gamma} $$

。例如，番茄酱在没有外力作用时不会流动（有屈服应力），而游泳池中的阻力比清水大，这是因为添加了增稠剂使其成为拟塑性流体。

🙋

"表观粘度"变化的机制是什么？分子水平上发生了什么？

🎓

以高分子溶液为例，在低剪切速率时，高分子链相互缠绕会大幅阻碍流动。当剪切速率升高时，链被拉伸并沿流动方向取向，从而阻力减少。这就是拟塑性（剪切变稀）。相反，在悬浊液中，高剪切区域的颗粒可能"卡住"产生更大阻力，这是膨胀型行为。数值上，Carreau模型用概率论方法处理链的拉伸。

🙋

控制方程仍然是Navier-Stokes方程吗？还是需要修改？

🎓

连续方程和动量守恒法则的形式相同，但应力张量τ的本构关系（Constitutive Equation）完全不同。假设不可压缩性，动量方程为

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{g} $$

。对于牛顿流体，

$$ \boldsymbol{\tau} = 2 \mu \mathbf{D} $$

（D是变形速率张量），但非牛顿流体中这个简单关系不成立。这个本构关系的建模是一切的出发点。

🙋

有哪些著名的本构关系模型？最简单的是什么？

🎓

幂律（Power-law）模型最基础：

$$ \eta(\dot{\gamma}) = K \dot{\gamma}^{n-1} $$

。其中K是粘度系数[Pa·s^n]，n是幂指数。n<1表示拟塑性（剪切变稀），n>1表示膨胀型流体，n=1则回到牛顿流体。例如，0.3%羧甲基纤维素水溶液n≈0.6，而浓玉米淀粉悬浊液的n>1表现出膨胀特性。但该模型在极低或极高剪切速率下会失效。

非牛顿流体的数值计算方法

离散化和求解器设置

🙋

在非牛顿流体CFD分析中，由于粘度在各处变化，离散化需要特殊处理吗？

🎓

是的，关键问题是"表观粘度"是速度梯度（剪切速率）的函数，导致强烈的非线性。在有限体积法中，如何在单元界面处插值粘度η_face至关重要。简单的算术平均不够准确，需要在界面处重新计算剪切速率，或采用逆风方法。此外，剪切速率

$$ \dot{\gamma} = \sqrt{2 \mathbf{D}:\mathbf{D}} $$

的计算涉及速度梯度张量的不变量，需要在单元中心和界面保持一致性。

🙋

求解器的收敛似乎会变差，有没有除牛顿法外的有效解法？

🎓

伪时间推进法（Pseudo-transient）很有效。即使是定常分析，也可以添加非物理的阻尼项来求解演化方程，直到收敛到定常状态。OpenFOAM的`simpleFoam`非牛顿对应求解器采用了"分裂解法"：粘度在外层迭代中更新，而内层的速度-压力迭代（SIMPLE）中粘度保持不变。对于n非常小（如0.3）的幂律模型，可以用"延续法"：从牛顿流体（n=1）开始，逐步降低n的值。

🙋

对于有屈服应力的Bingham流体，未屈服区域（固体行为）与屈服区域的边界如何处理？

🎓

这是最困难的问题。Bingham模型中，只有应力超过屈服应力τ_y的区域才会流动。数值方法常采用"粘度正则化"，将伪固体区域视为具有极高粘度的流体。例如，Papanastasiou正则化模型为：

$$ \eta(\dot{\gamma}) = \mu_p + \frac{\tau_y}{\dot{\gamma}} \left[ 1 - \exp(-m \dot{\gamma}) \right] $$

。参数m越大越接近理想Bingham流体，但数值刚性也会增加。通常从m=1000 1/s左右开始。

非牛顿流体的实务应用

工作流程和检查清单

🙋

开始进行非牛顿流体分析时，首先需要确定什么？

🎓

绝对需要的是"实测流动曲线数据"。要明确剪切速率的范围。例如，圆管内流的壁面剪切速率可概算为

$$ \dot{\gamma}_w = \frac{8V}{D} $$

（V是平均流速）。用流变仪测量的τ-γ̇数据拟合到幂律、Carreau、Herschel-Bulkley等模型。如果拟合的决定系数R²低于0.99，需要质疑模型的适用范围，或在不同范围内分段使用不同模型。

🙋

网格划分有什么注意事项？特别是壁面附近？

🎓

壁面处剪切速率最高，粘度变化最剧烈，所以需要捕捉粘度剧变区域（如Bingham流体的屈服面附近）。牛顿流体的y+<1准则不能直接使用。建议先以牛顿流体模拟，查看得到的剪切速率分布，然后在粘度急变区域（如圆管中心剪切速率为0处）加密网格。特别是圆管中心剪切速率变为零时，要确保粘度不发散，模型的正则化能正常工作。

🙋

没有实验数据的情况下如何验证结果的可信度？

🎓

没有实验也要做"基准测试"与理论解或已知相关式对比。例如，圆管内幂律流体的发达流压力损失可以理论计算（使用广义雷诺数Re_g）。网格依赖性研究中，将网格数加倍后压力损失或中心流速的变化应小于1%。还要检查求解器残差下降时，整体力平衡（流入动量与压力梯度+壁面剪切力）是否守恒。

非牛顿流体的软件对比

Ansys, Abaqus, COMSOL等的实现

🙋

在Ansys Fluent中处理非牛顿流体时，有哪些标准可用的模型？

🎓

Fluent内置了丰富的非牛顿模型：幂律、Carreau、Herschel-Bulkley（带屈服应力）、Cross、Bird-Carreau等。用户定义函数（UDF）可实现自定义模型。实际应用中，高分子熔融树脂的注射成形仿真常使用"Modified Cross Model"。设置时需对剪切速率设置下限和上限裁切值（如1e-6 1/s和1e6 1/s）防止数值发散。求解器通常采用"Pressure-Based"，梯度重建方案用"Green-Gauss Node Based"比"Least Squares Cell Based"更稳定。

🙋

COMSOL Multiphysics以"多物理场"著称，在非牛顿流体耦合中有什么优势？

🎓

COMSOL的强项是能够完全耦合地求解温度变化导致的粘度变化或化学反应导致的组成变化与流体流动。"非等温非牛顿流动"接口中，粘度模型可结合Arrhenius型温度依赖性

$$ \eta(T) = A \exp(E_a / RT) $$

与剪切速率依赖性一起定义。对于涂布工艺中的自由表面问题，可将Level Set法或Phase Field法与表面张力及非牛顿性完全耦合求解。用户定义的粘度函数也可直接用MATLAB Live Link编写。

🙋

结构分析软件Abaqus中能否处理非牛顿流体？

🎓

可以，特别在"型腔充填分析"和"橡胶材料超弹性+粘弹性分析"中很重要。Abaqus有CFD模块，但通常在Standard或Explicit中将"粘性"材料模型定义相关参数。例如，树脂注射成形仿真（采用Abaqus/Explicit的CEL法）常用7系数的Cross-WLF模型作业界标准。该模型同时依赖剪切速率和温度：

$$ \eta(\dot{\gamma}, T) = \frac{\eta_0(T)}{1+(\eta_0 \dot{\gamma} / \tau^*)^{1-n}} $$

。材料数据通常从Moldflow等专用软件导入。

非牛顿流体的故障排除

常见错误和解决方案

🙋

分析时出现"检测到负的剪切速率"的错误。剪切速率本来不应该是负数？

🎓

这是数值舍入误差或发散过程中的非物理速度场导致的。单元内变形速率张量D的第2不变量计算结果可能出现微小的负值。两个对策：1) 在剪切速率计算中加入小正数ε进行裁切：

$$ \dot{\gamma} = \sqrt{2\mathbf{D}:\mathbf{D} + \epsilon^2} $$

（ε=1e-6 1/s等）。2) 从根本上重新检查初始条件和边界条件，在流速急剧变化处加密网格，提高速度梯度计算精度。Fluent可以设置"Shear Rate Limit"，这是最简便的方法。

🙋

幂律模型中n很小（如0.2）时，求解器残差振荡，根本不收敛。怎么办？

🎓

n小时，低剪切速率区域的粘度会非常高，方程的刚性急剧增加。解决步骤：1) **初始化**：先用n=1（牛顿流体）得到完全收敛的解。2) **延续法**：以该解为初值，逐步降低n至0.9, 0.8, ...，多个求解步骤连续进行。多数求解器的"参数研究"功能可自动化。3) **放松系数**：大幅降低速度和压力的欠松弛系数（≤0.2）。4) **模型调整**：实际上，在低剪切区粘度收敛到有限值的"Carreau模型"等更合适。

🙋

Bingham流体分析时，本不应屈服的区域（栓塞区）数值上有微小流动。这是允许误差吗？

🎓

这是"正则化误差"。前面提过的Papanastasiou模型中，由于m是有限值，在剪切速率接近零的区域粘度不会无穷大，出现微弱的伪蠕变流。评估误差方法是：栓塞区最大流速V_plug与代表流速V的无量纲比

$$ V_plug / V < 0.01 $$

即可接受。如果要减少误差可增大正则化参数m，但会导致求解器收敛变差。工程上的判断是：这种微小流动是否影响工程结论（压力损失、混合效率等）。

🙋

实现了用户定义粘度函数后，计算成本异常上升。原因是什么？

🎓

主要有两个原因：1) **函数内的条件分支（if语句）**：按"剪切速率是否超过阈值"分支会破坏矢量化，大幅降低CPU缓存效率。尽可能用光滑函数（如tanh）来避免分支。2) **高成本的数学函数**：exp()、log()、pow()计算代价高，特别是pow(x, n)在n为非整数时内部执行exp(log(x)*n)。可预先建立查询表（Look-up Table），用插值求粘度，性能可提升10倍以上。Ansys Fluent的UDF中，`DEFINE_PROPERTY`宏用`real shear_rate = sqrt(2.0*C_DOT(c,t)*C_DOT(c,t));`计算剪切速率本身也有成本。

非牛顿流体 — CAE用语解释