L4 · 第 7 章

四辊流的驻点怎样失去对称，数值方法又改变了多少

Sandbox Physics · 2026-10-29

用一个定常、对称的力推动没有惯性的高分子溶液，流动能否自己打破对称、开始振荡？Thomases 与 Shelley 的四辊模型给出了肯定的数值答案：Weissenberg 数约 5 以上，驻点处的应力细带会偏向一侧，再往上流动持续振荡并混合流体。后来的一项研究认为，这一对称破缺至少部分来自数值计算里加入的应力扩散。课程把两边的证据摆出来，也给出判断一个数值结果是否可信的办法。

四个涡与一个驻点

计算区域是边长 $2\pi$ 的周期方盒。定常的外力 $\mathbf f$ 产生四个涡，中心是一个驻点：流体沿一个方向流入，沿另一个方向流出，和第三章的平面伸长一样。没有高分子时，速度正好是外力的一半。

时间单位 $\tau_f$ 取牛顿流动在驻点处应变率的倒数，Weissenberg 数就是弛豫时间与它之比。高分子与溶剂的耦合 $\beta$ 随 Weissenberg 数调整， $\beta\mathrm{Wi}=1/2$ 保持不变，也就是高分子与溶剂黏度之比固定为 0.5，以便对照 Arratia 等人的十字槽实验。应力方程里还有一项扩散 $\nu_p\Delta\mathbf S$ ，下文会专门讨论它。

-\nabla p+\Delta\mathbf u=-\beta\nabla\cdot\mathbf S+\mathbf f,\quad \nabla\cdot\mathbf u=0,\qquad \overset{\nabla}{\mathbf S}+\mathrm{Wi}^{-1}\left(\mathbf S-\mathbf I\right)=\nu_p\Delta\mathbf S,\qquad \mathbf f=\left(2\sin x\cos y,\,-2\cos x\sin y\right)

Thomases、Shelley 与 Thiffeault 的式 (1)、(2) 与 (6)；构象张量在文中记作 S。

来源：B. Thomases · B. Thomases and M. Shelley

驻点上的应力细带

驻点附近速度近似为 $\alpha\,(x,-y)$ ，高分子感受到的是有效 Weissenberg 数 $\alpha\mathrm{Wi}$ 。高分子应力会反过来减慢驻点处的流动：Weissenberg 数很小时 $\alpha$ 趋于 2/3，Weissenberg 数增大时还会更小。

Thomases 与 Shelley 2007 年不加应力扩散计算这个模型，网格最多每边 4096 点。有效值达到 1/3 时，应力沿流出方向的剖面出现尖点，对应 $\mathrm{Wi}\approx0.5$ ；有效值达到 1/2 时应力无界增长，对应 $\mathrm{Wi}\approx0.9$ 。后者正是第三章平面伸长的阈值。

高分子被拉长的区域因此收缩成沿流出方向的一条细带，越来越窄、越来越强。换成 FENE-P 后应力有界，但尖点状的结构仍然存在。

\mathbf u\approx\alpha\,(x,\,-y),\qquad \mathrm{Wi}_{\mathrm{eff}}=\alpha\,\mathrm{Wi}

驻点附近的局部线性流动，Thomases 与 Shelley 2007 年的写法。

来源：B. Thomases and M. Shelley

对称破缺与振荡

加上 $\nu_p=10^{-3}$ 的应力扩散后，Thomases、Shelley 与 Thiffeault 先求出对称的定常解，再加一个小扰动： $\mathrm{Wi}=4.8$ 时扰动衰减， $4.9$ 时扰动增长（图 2）。他们用 $256\times256$ 的网格，并报告分辨率加倍后动力学不变。

约 4.9 到 9 之间，失去对称的流动经过一段振荡后重新稳定在一个不对称的定常态。约 9 以上，一个涡的中心脱离固定位置，流动持续振荡，并开始混合四个象限的流体。 $\mathrm{Wi}=10$ 时速度谱有两个周期，约 330 和 55 个 $\tau_f$ （图 6(d)）；2009 年的短文报告的是 317 和 59，并说对称约在 $t\approx200$ 时破缺。Weissenberg 数为 12 和 15 时运动是周期的，约 20 以上变为非周期。

两篇文章都用聚合物应变能的平衡式检查计算：能量 $E$ 取 $\tfrac12\operatorname{tr}(\mathbf S-\mathbf I)$ 的积分，它的变化由外力做功、黏性耗散和弛豫决定。

\dot E+\mathrm{Wi}^{-1}E=-\mathcal D+\mathcal W,\qquad E=\frac12\int\operatorname{tr}\left(\mathbf S-\mathbf I\right)\,d\mathbf x

Thomases、Shelley 与 Thiffeault 的式 (4)—(5)；耗散与做功的具体写法见原文。

来源：B. Thomases · B. Thomases and M. Shelley

应力扩散带来的疑问

应力扩散让方程在细带处保持光滑，代价是改变了物理。Gupta 与 Vincenzi 指出， $\nu_p=10^{-3}$ 对应的 Schmidt 数约为 $10^3$ ，比真实高分子溶液小三个数量级。他们用类似四辊的胞状外力、Deborah 数 10 做对照：加这样的扩散时，流动长时间几乎冻结，随后打破胞状对称；用不加扩散、能保持正定的格式时，流动很快变得混乱，却一直保持外力的对称。

他们的结论是：在他们的设置中，对称破缺是扩散造成的假象；对应变与涡量混在一起的外力，这种影响可能弱一些。Thomases、Shelley 与 Thiffeault 在附录 B 中按 FENE-P 标定了扩散系数，但没有报告阈值随扩散系数怎样变化。

两组计算的外力形式和数值格式都不同，目前还不能断定四辊模型的对称破缺有多少是物理的。这是一个尚未解决的问题。Datta 等人的综述也提醒，Oldroyd-B 在高 Weissenberg 数下几乎从来不是定量准确的。

来源：A. Gupta and D. Vincenzi · B. Thomases · S. S. Datta et al.

怎样检查一个数值结果

一个数值上的阈值或周期，要经过几项独立的检查，才能当作模型的性质。第一，网格加密一倍，结果不变。第二，时间步长减半，结果不变。第三，改变应力扩散系数，至少再取一个更小的值；能做到的话，换一个不需要人工扩散的格式。第四，全程监视构象张量的最小本征值，它必须保持为正；对数构象格式能自动保证这一点。第五，换成 FENE-P 这样有界的模型，看结论是否保留。

还有一项便宜的检查：累计能量平衡式两边的差。差值不随时间增长，说明时间积分没有悄悄注入或吞掉能量。课程配套的四辊实验把这些量显示在页面上，每次运行都记录网格、时间步长、扩散系数和扰动种子。它是一个比原文粗的离散，所以这些检查正是为它准备的。

来源：S. S. Datta et al. · Sandbox Physics FR26

完整例题：把论文的数换成弛豫时间

时间以 $\tau_f$ 计，弛豫时间是 $\mathrm{Wi}\,\tau_f$ 。 $\mathrm{Wi}=10$ 时，周期 330 和 55 分别是 33 和 5.5 个弛豫时间；2009 年报告的对称破缺时刻 200 约为 20 个弛豫时间。慢的那个周期比高分子自身的记忆长得多，说明它是整个流动的集体行为。

扩散在一个弛豫时间内抹平的距离约为 $\sqrt{\nu_p\,\mathrm{Wi}\,\tau_f}$ 。 $\nu_p=10^{-3}$ 、 $\mathrm{Wi}=10$ 时是 $\sqrt{10^{-2}}=0.1$ ，而方盒边长为 $2\pi\approx6.3$ 。论文在 Weissenberg 数大于 10 时取 $\nu_p=0.01/\mathrm{Wi}$ ，乘积 $\nu_p\mathrm{Wi}$ 固定为 0.01，这个距离始终是 0.1。

阈值的精度：4.8 衰减、4.9 增长，把阈值限定在 $4.85\pm0.05$ ，约 1%。这是在一个网格、一个扩散系数下的结果；上一节的检查是把它推广到“模型的阈值”之前要做的事。

来源：B. Thomases · B. Thomases and M. Shelley

到实验里试一试

找对称破缺，并检查它

先选“对称的应力细带”预设（Weissenberg 数 4），再选“对称破缺后归于稳定”（Weissenberg 数 6）。观察不对称量的曲线：一个衰减到零附近，另一个在对数坐标上沿直线增长吗？
改为 Weissenberg 数 10，记录探针速度的振荡周期，以 $\tau_f$ 为单位，与论文的约 330 和 55 比较。写下网格大小和应力扩散系数。
保持 Weissenberg 数 10，依次只改一项：网格加密、时间步长减半、应力扩散增大。哪一项改变了结果？再看拉伸的最小特征值：在默认的 64 乘 64 网格上从静止开始时，应力带形成期间它会低于零，说明网格分辨不了应力带；在 128 乘 128 网格上它保持为正。

对照均匀平面伸长的阈值

用 Oldroyd-B 在 Weissenberg 数 0.45 和 0.55 下运行，确认伸长流中的伸展只在后者无界增长。这就是四辊驻点处有效 Weissenberg 数达到 $1/2$ 时发生的事。

练习与解答

四辊模型中，若驻点处的 $\alpha$ 一直等于 2/3，Weissenberg 数多大时有效值达到 1/2？Thomases 与 Shelley 实际得到约 0.9，为什么更大？
展开解答
$\mathrm{Wi}=(1/2)/(2/3)=0.75$ 。实际值更大，是因为高分子应力随 Weissenberg 数增大而进一步减慢驻点处的流动， $\alpha$ 降到 2/3 以下，要更大的 Weissenberg 数才能让乘积达到 1/2。
取 $\nu_p=0.01/\mathrm{Wi}$ ，在 Weissenberg 数为 20 时，一个弛豫时间内的扩散距离是多少？在慢周期 330 内呢（按 10 时的周期估计）？
展开解答
一个弛豫时间内 $\sqrt{(0.01/20)\times20}=0.1$ 。若周期仍按 330 估计， $\sqrt{5\times10^{-4}\times330}=0.41$ ，已是方盒边长的约 7%。
要说服别人“对称破缺是模型本身的性质，而非应力扩散的产物”，至少要给出哪些证据？
展开解答
同一结论在网格加密和时间步长减半后不变；在更小的应力扩散系数下，或在不需要人工扩散的格式中仍然出现；构象张量全程正定；最好在 FENE-P 中也能看到。只在一种扩散系数下得到的阈值，只能算这一离散方法的结果。
为什么 Thomases 等人保持 $\beta\mathrm{Wi}=1/2$ ，而不是固定 $\beta$ ？
展开解答
这样高分子与溶剂的黏度之比固定为 0.5，改变 Weissenberg 数时只改变弹性，不改变溶液的黏度组成，便于对照 Arratia 等人的十字槽实验。Thomases 与 Shelley 的解释是：对一种给定的工作液，这个比值不随实验条件改变。

来源与复算范围

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 and M. Shelley, Emergence of singular structures in Oldroyd-B fluids, Physics of Fluids 19, 103103 (2007)
不加应力扩散、由四辊力驱动的 Stokes–Oldroyd-B 流。中心驻点附近，局部应变率乘 Weissenberg 数起有效 Weissenberg 数的作用；Weissenberg 数趋于零时这一应变率趋于 2/3。有效值的两个临界点 1/3 和 1/2，分别对应 Weissenberg 数约 0.5（应力出现尖点）和约 0.9（应力无界）。网格最多每边 4096 点。加上 FENE-P 惩罚后应力有界，但尖点状结构仍在。
B. Thomases and M. Shelley, Transition to mixing and oscillations in a Stokesian viscoelastic flow, Physical Review Letters 103, 094501 (2009)
式 (1)—(2)：周期盒中的 Stokes 流加 Oldroyd-B 应力，由定常四辊力驱动；聚合物与溶剂黏度之比固定为 0.5，对应 Arratia 等人的十字槽实验；应力扩散取 1e-3。Weissenberg 数为 10 时，流动约在论文单位的时刻 200 失去对称，之后出现周期约 317 和 59 的振荡；低于约 5 时流动保持被外力“奴役”的对称形态。
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 标定；文中没有报告它对阈值的影响。
A. Gupta and D. Vincenzi, Effect of polymer-stress diffusion in the numerical simulation of elastic turbulence, arXiv:1809.09648v2 (2019)
2019 年 3 月 7 日 v2；未核对期刊版本。二维 Oldroyd-B 流，由胞状外力驱动，Deborah 数 10，采用 Cholesky 分解与 Kurganov–Tadmor 格式。使用 Thomases 等人的应力扩散（Schmidt 数 1e3）时，流动长时间几乎冻结，随后打破胞状对称（图 2—4）；不加扩散时很快变为混沌，却保持外力的对称（图 5）。第 3 节：在他们的设置中，对称破缺是扩散造成的假象。第 4 节：FENE-P 结果相同；对应变与涡量混在一起的外力，这种影响可能较弱。
Sandbox Physics FR26, four-roll elastic instability: numerical model (docs/frontier/fr26/four-roll-elastic-instability-model.md, with long-run records in docs/frontier/fr26/records/)
模型版本 fr26-four-roll-elastic-instability-v1。模型说明写明：双周期的四辊单元代替了 2001 年实验中的弯曲流道，它有弯曲流线和拉伸驻点，但没有壁面、没有净流量、没有流道宽度和半径，因此不能拿流道的阈值和频谱与它比较。说明中记录了方程、伪谱离散方法和实验页面显示的诊断量。每次运行都与网格、时间步长、应力扩散和扰动种子一起保存。这些是本路线的数值结果，不是实验数据。