L3 · 第 4 章

FENE-P 模型怎样用有限链长给伸展和伸长黏度设上限

Sandbox Physics · 2026-10-29

线性弹簧可以无限伸长，所以 Oldroyd-B 的伸长应力在阈值以上无界增长。真实的链有轮廓长度，拉到接近全长时弹簧会变得极硬。FENE-P 模型把这一点写进构象方程：伸展越接近上限，回复得越快。结果是伸长黏度停在一个有限的平台上，剪切黏度随剪切率下降，而弹性能、外力做功和耗散之间仍有一条精确的平衡式。

有最大长度的弹簧

真实高分子链的端到端距离不能超过它的轮廓长度。FENE 弹簧的张力在接近这个长度时发散。Peterlin 近似把张力里依赖长度的因子换成它在平均构象处的值，得到一个只依赖 $\operatorname{tr}\mathbf A$ 的因子 $f$ 。 $L^2$ 是 $\operatorname{tr}\mathbf A$ 的上限，即最大伸长的平方。

平衡时 $\operatorname{tr}\mathbf A=3$ ， $f=1$ ，模型与 Oldroyd-B 相同； $\operatorname{tr}\mathbf A$ 接近 $L^2$ 时 $f$ 发散，弹簧回复的速率 $f/\lambda$ 可以追上任何拉伸率。 $L^2\to\infty$ 时回到 Oldroyd-B。

\overset{\nabla}{\mathbf A}=-\frac{1}{\lambda}\left(f\mathbf A-\mathbf I\right),\qquad \boldsymbol\tau_p=G\left(f\mathbf A-\mathbf I\right),\qquad f=\frac{L^2-3}{L^2-\operatorname{tr}\mathbf A}

Morozov 与 Spagnolie 的式 (1.86)—(1.87) 用应力写出同一模型，其中因子为

1+\lambda\operatorname{tr}\boldsymbol\tau_p/(\eta_pL^2)

；代入应力与构象的关系即得此式。Datta 等人表 II 采用这里的写法。

来源：A. Morozov and S. E. Spagnolie · S. S. Datta et al.

伸长黏度停在平台上

平面伸长的稳态与第三章一样分量分离，只是弛豫率乘上了 $f$ 。三个对角分量各自写出后， $f$ 又依赖它们的和，所以要解一个自洽方程；用二分法在 $f>2\mathrm{Wi}$ 的范围里求根即可。

强拉伸时 $A_{xx}$ 很大，分母 $f-2\mathrm{Wi}$ 必须很小，所以 $f$ 只比 $2\mathrm{Wi}$ 稍大：弹簧自动变硬，直到回复率正好抵消拉伸率。由 $f(A_{xx}-A_{yy})=2\mathrm{Wi}(A_{xx}+A_{yy})$ 这一恒等式，伸长黏度可以写成下式。 $\mathrm{Wi}$ 很大时 $A_{xx}$ 趋近 $L^2$ ，伸长黏度趋于上限。

McKinley 的式 (7) 给出单轴伸长中拉伸应力的上限 $2\eta_p\dot\varepsilon L^2$ ，第 5.4 节相应的伸长黏度上限是 $3\eta_s+2\eta_pL^2$ ；平面伸长的上限把 3 换成 4。

A_{xx}=\frac{1}{f-2\mathrm{Wi}},\quad A_{yy}=\frac{1}{f+2\mathrm{Wi}},\quad A_{zz}=\frac1f,\qquad \eta_P=4\eta_s+2\eta_p\left(A_{xx}+A_{yy}\right)\xrightarrow{\ \mathrm{Wi}\to\infty\ }4\eta_s+2\eta_pL^2

稳态时

f

由

f\left(L^2-\operatorname{tr}\mathbf A\right)=L^2-3

决定。

来源：G. H. McKinley · A. Morozov and S. E. Spagnolie

剪切中的变稀

稳态剪切中 $A_{yy}=1/f$ ，高分子剪切应力为 $G\,\mathrm{Wi}/f$ ，剪切黏度因此等于 $\eta_s+\eta_p/f$ 。剪切越快，线团越长， $f$ 越大，黏度越小。Oldroyd-B 的剪切黏度不变，FENE-P 的剪切黏度随剪切率下降。

把剪切的稳态分量代入自洽条件，得到一个关于 $f$ 的三次方程。 $L^2$ 越大，同样的 $\mathrm{Wi}$ 下 $f$ 越接近 1，变稀越弱。

\eta=\eta_s+\frac{\eta_p}{f},\qquad f^2\left(f-1\right)=\frac{2\,\mathrm{Wi}^2}{L^2}

由稳态剪切的

A_{yy}=1/f

、

A_{xy}=\mathrm{Wi}/f^2

、

A_{xx}=1/f+2\mathrm{Wi}^2/f^3

代入自洽条件得到。

来源：A. Morozov and S. E. Spagnolie

弹性能、做功与耗散

FENE-P 有一个弹性自由能 $\psi$ ，平衡时为零，其他构象下为正。跟着流体看，它的变化率等于流动对高分子做的功率 $\boldsymbol\tau_p:\mathbf L$ 减去耗散 $\mathcal D$ ；耗散写成本征值 $a_i$ 的平方和，永不为负。

这条平衡式给数值解一个独立的检查：在一段时间内累计功、储存的能量和耗散，三者之差应接近零。Thomases、Shelley 与 Thiffeault 的式 (4)—(5) 对 Oldroyd-B 写了同类平衡式，只是用的能量是 $\tfrac12\operatorname{tr}(\mathbf S-\mathbf I)$ 的积分。

\psi=\frac{G}{2}\left[-(L^2-3)\ln\frac{L^2-\operatorname{tr}\mathbf A}{L^2-3}-\ln\det\mathbf A\right],\qquad \frac{D\psi}{Dt}=\boldsymbol\tau_p:\mathbf L-\mathcal D,\qquad \mathcal D=\frac{G}{2\lambda}\sum_i\frac{\left(fa_i-1\right)^2}{a_i}

课程自行推导：对

\psi

求导，代入构象方程，并用不可压缩条件

\operatorname{tr}\mathbf L=0

。

L^2\to\infty

时

\psi

变为

\tfrac G2\left(\operatorname{tr}\mathbf A-3-\ln\det\mathbf A\right)

。

来源：B. Thomases

数值解要守住的条件

构象张量必须正定，FENE-P 还要求 $\operatorname{tr}\mathbf A<L^2$ 。离散误差可能让最小本征值变成负数，或让迹越过上限，此后的应力就没有物理意义。Morozov 与 Spagnolie 第 6 节提醒，Oldroyd-B 模拟在高 Weissenberg 数下常常失败，这在文献中称为高 Weissenberg 数问题。

常见的对策有三种。一是对构象张量取矩阵对数再求解，即对数构象格式，它自动保持正定；二是给应力方程加一点扩散，抹平过陡的梯度；三是改用 FENE-P 这类有上限的模型。Datta 等人的综述讨论了第一种，并指出 Oldroyd-B 在高 Weissenberg 数下几乎从来不是定量准确的。第二种办法的代价会在第七章看到。

来源：A. Morozov and S. E. Spagnolie · S. S. Datta et al.

完整例题：平面伸长中的两种链长

取 $\beta=0.5$ ， $\mathrm{Wi}=3$ 。Oldroyd-B 在这里没有稳态，伸展会指数增长。FENE-P 取 $L^2=400$ ，二分法解出 $f=6.003$ ，只比 $2\mathrm{Wi}=6$ 大 0.003。于是 $A_{xx}=333.6$ ， $A_{yy}=0.0833$ ， $A_{zz}=0.167$ ，迹 333.9，用掉了上限的 83.5%。伸展比 $\sqrt{333.9/3}=10.5$ ，最大可能值是 $\sqrt{400/3}=11.5$ 。

伸长黏度 $\eta_P=4\times0.5+2\times0.5\times(333.6+0.083)=335.7\,\eta_0$ ，上限为 $4\times0.5+400=402\,\eta_0$ 。

改为 $L^2=100$ ： $f=6.012$ ， $A_{xx}=83.6$ ，迹 83.9，用掉 83.9%，伸展比 5.29，伸长黏度 $85.7\,\eta_0$ ，上限 $102\,\eta_0$ 。链长的平方缩小到四分之一，伸长黏度也大致缩小到四分之一；用掉上限的比例几乎不变。

能量检查：取 $\lambda=1\,\mathrm{s}$ 、 $G=0.5\,\mathrm{Pa}$ ， $L^2=400$ 时流动对高分子做功的功率是 $(\tau_{xx}-\tau_{yy})\dot\varepsilon=3.00\times10^{3}\,\mathrm{W/m^3}$ ，按本征值算出的耗散与它相同。稳态下弹性能不再变化，做的功全部耗散掉。

来源：A. Morozov and S. E. Spagnolie

到实验里试一试

让链被拉到全长

选“链被拉到全长”预设（FENE-P，Weissenberg 数 3， $L^2=100$ ）。运行结束时，伸长流中的伸展应接近例题的 5.29，“用掉的有限伸长上限的最大比例”应接近 0.84。
把最大伸长的平方改为 400，预测伸展和伸长黏度各变成多少，再对照例题。看详细读数中“功减储存减耗散的最大残差”，确认它远小于做功本身。
选“恒定应力自动刹车”预设。给定应力时，伸长黏度一旦变大，同样的应力就只能维持较小的拉伸率。观察伸长一侧的 Weissenberg 数怎样随时间下降。

练习与解答

FENE-P 溶液 $\beta=0.1$ 、 $L^2=1000$ 。平面伸长黏度的上限是零剪切黏度的多少倍？是牛顿值的多少倍？
展开解答
上限为 $4\times0.1+2\times0.9\times1000=1800.4$ 倍零剪切黏度；牛顿值是 $4\eta_0$ ，所以约为牛顿值的 450 倍。
取 Oldroyd-B 极限的弹性能，令 $\mathbf A=\mathrm{diag}(a,1/a,1)$ 。写出 $\psi$ ，并说明它为什么不小于零。
展开解答
$\operatorname{tr}\mathbf A-3=a+1/a-2$ ， $\det\mathbf A=1$ ，所以 $\psi=\tfrac G2\left(a+1/a-2\right)$ 。由算术—几何平均不等式 $a+1/a\ge2$ ，只在 $a=1$ 时取等号。
FENE-P， $\beta=0.5$ 、 $L^2=100$ ，稳态剪切 $\mathrm{Wi}=3$ 。估计剪切黏度。
展开解答
解 $f^2(f-1)=2\times9/100=0.18$ ，得 $f\approx1.14$ 。剪切黏度 $0.5+0.5/1.14=0.94\,\eta_0$ ，只降了约 6%；同一 Weissenberg 数下伸长黏度却是零剪切黏度的 86 倍。
为什么强拉伸时 $f$ 总是只比 $2\mathrm{Wi}$ 稍大？
展开解答
稳态 $A_{xx}=1/(f-2\mathrm{Wi})$ 。强拉伸时 $A_{xx}$ 接近 $L^2$ ，是个大数，分母必须很小。物理上，弹簧不断变硬，直到回复率 $f/\lambda$ 正好抵消流动的放大率 $2\dot\varepsilon$ 。

来源与复算范围

A. Morozov and S. E. Spagnolie, Introduction to complex fluids, in Complex Fluids in Biological Systems (Springer, 2015), chapter 1
第 3 节：Weissenberg 数定义为弛豫时间乘剪切率；式 (1.27) 为广义牛顿液体。第 4.1 节：式 (1.47)—(1.51)，弹簧、阻尼器、Maxwell 液体及其记忆积分，弛豫时间等于黏度除以模量。式 (1.72)：上随体导数；式 (1.80)：Oldroyd-B；式 (1.86)—(1.87)：FENE-P。第 4.4 节：Hooke 哑铃，式 (1.98)—(1.110)。式 (1.116)：Oldroyd-B 稳态剪切的法向应力。式 (1.117)—(1.125)：环向应力与锥板法向力。式 (1.128)—(1.131)：单轴伸长，以及 Oldroyd-B 在伸长率达到弛豫时间倒数一半时的发散。书中没有启动过程的解，也没有平面伸长；这两部分由课程自行推导。
G. H. McKinley, Visco-elasto-capillary thinning and break-up of complex fluids, MIT Hatsopoulos Microfluids Laboratory report 05-P-04, submitted to Annual Rheology Reviews (2005)
式 (1)—(2)：变细液丝的应变率与 Hencky 应变。式 (4)—(7)：单轴伸长中 Oldroyd-B 的应力增长，以及有限伸长带来的上限。式 (8)：弹性—毛细变细，半径按指数下降，时间常数为三倍弛豫时间，液丝的 Weissenberg 数为 2/3。式 (11)：伸长黏度不变后线性变细。式 (15)：颈部的零维力平衡。式 (16)：牛顿液体的相似解，系数 0.0709；其后的正文指出简单零维模型的估计偏差超过 100%。第 5.4 节：有限伸长哑铃的终极伸长黏度。式 (8) 的前置系数与表 2 所列相差一个因子，即 2 的立方根，第 5.4 节把差别归于一维相似解；课程只用指数，不用前置系数。
S. S. Datta et al., Perspectives on viscoelastic flow instabilities and elastic turbulence, arXiv:2108.09841v1 (2021); Physical Review Fluids 7, 080701 (2022)
arXiv v1，2021 年 8 月 22 日；未核对期刊版本。第 II 节式 (1)：Pakdel–McKinley 判据，把流线曲率与法向应力和剪切应力之比结合起来。文中没有给出 Taylor–Couette 或锥板流的临界值。第 III 节：十字槽流动的模拟中，这一判据的局部最大值超过约 4 至 5 时，流动失去对称或变得随时间变化；式 (2) 与表 II 用构象张量写出 Oldroyd-B、FENE-P 等模型，便于检查正定性并采用对数构象格式；一个二维 Taylor–Couette 模拟在 Weissenberg 数约 10 处变为混沌。研讨会总结指出，在高 Weissenberg 数下 Oldroyd-B 几乎从来不是定量准确的。
B. Thomases, M. Shelley and J.-L. Thiffeault, A Stokesian viscoelastic flow: transition to oscillations and mixing, Physica D 240, 1602 (2011)
式 (1)—(2) 与 (6)：同一模型，Weissenberg 数不超过 10 时应力扩散取 1e-3。式 (4)—(5)：聚合物应变能的平衡式。第 3.1 节图 2：对称定常态加小扰动，Weissenberg 数为 4.8 时衰减，4.9 时增长；网格 256 乘 256，分辨率加倍动力学不变。图 6(d)：Weissenberg 数为 10 时速度谱的周期约 330 和 55。第 5 节：约 4.9 到 9 之间，失去对称的流动会重新稳定；约 9 以上一个涡旋中心脱离固定位置，持续振荡并混合流体；12 和 15 时为周期运动，约 20 以上为非周期。附录 B：应力扩散按 FENE-P 标定；文中没有报告它对阈值的影响。