黏弹性流体数值模拟实践

本研究采用交错网格有限差分法离散求解控制方程，称之为LCR方法。为方便对比，本文也直接求解Oldroyd-B模型，并称之为CR方法。本构方程对流项求解方法可采取Kurganov–Tadmor（KT）格式、MINMOD格式与中心差分（CD）格式，求解方法可以参考我们的近期研究文献。

[1] Kurganov A , Tadmor E . New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection–Diffusion Equations[J]. Journal of Computational Physics, 2000, 160(1):241-282.

[2] Yu Bo, Kawaguchi Y . Direct numerical simulation of viscoelastic drag-reducing flow: a faithful finite difference method[J]. Journal of Non-Newtonian Fluid Mechanics, 2004, 116(2-3):431-466.

[3] Zhang Wenhua, Zhang Hongna, et al. Comparison of turbulent drag reduction mechanisms of viscoelastic fluids based on the Fukagata-Iwamoto-Kasagi identity and the Renard-Deck identity[J]. Physics of Fluids, 2020, 32(1): 013104.

上述图片给出了Re=2857，Wi=1或5，β=1/9流动工况下的部分计算结果，槽道尺寸为4π×2×π或者8π×2×2π，网格分辨率为256×128×128。其中LCD代表在对数域实施中心差分格式，PMCD代表在物理域实施修正的中心差分格式，PCD代表在物理域实施中心差分格式。若LCR方法采用PCD离散格式，需要将物理域构象张量的插值结果映射到对数域。同理，CR方法采用LCD离散格式，需要将对数域构象张量的插值结果映射到物理域。

计算结果表明，在物理域直接离散情况下，CR方法可以模拟出黏弹性流体的湍流减阻现象，因为与牛顿流体相比各种离散格式下摩阻系数均有所减小，黏弹性微观结构得到拉伸使得湍流运动得以显著抑制。在对数域直接离散情况下，LCR方法难以模拟出湍流减阻现象。在物理域插值再映射到对数域情况下，LCR方法可以得到湍流减阻现象，但减阻效果较弱。另外，在对数域插值再映射到物理域情况下，CR方法得到的湍流减阻效果也有所减弱。

瞬时摩阻系数随时间的演化过程

 y⁺=15处x-z平面上牛顿流体与黏弹性流体（Wi=1）湍流瞬时流向速度场

Q方法提取的牛顿流体与黏弹性流体（Wi=5）湍流涡结构，Q=0.4Qrms

y-z平面上瞬时流向涡云图，Wi=1

y-z平面上瞬时黏弹性微观结构拉伸场

无论是对于LCR方法还是CR方法，物理域插值都能够更好地模拟出黏弹性流体湍流减阻现象，对数域插值方式下减阻剂的减阻效果较弱，LCR-LCD方法甚至几乎无法得到湍流减阻效果，减阻率仅为1%。另外，研究发现在对数域直接进行对数构象张量插值的LCR方法在顶盖驱动等流动中出现了难以收敛的现象。因此，我们认为对数域插值方式下LCR方法收敛性较差导致了该方法难以模拟出湍流减阻现象，而利用物理域插值可以部分地改善这一问题。但是LCR-PCD方法下的湍流减阻效果仍达不到CR-PMCD方法，说明LCR方法收敛性差还具有其他原因，例如本研究模拟中非稳态项使用的是Adams-Bashforth格式，该格式使用上一时层（n-1）与当前时层（n）非稳态项插值得到中间时层（n+0.5）非稳态项。LCR方法下直接在对数域实施该格式也可能会涉及与界面插值类似的插值域问题。本研究尝试在一阶显示格式下实施LCR-PCD方法，尽管能够进一步提高湍流减阻率，但是得到构象张量场越界现象十分严重，说服力较差，因而没有在本文中进行展示。

众所周知，守恒型方程的计算收敛性强于非守恒型方程。LCR方法在离散对数构象张量控制方程时通常采用守恒型方程来增强计算收敛性。需要强调的是对数构象张量控制方程本身不具有直接的物理意义，对数构象张量通量在控制容积界面处的守恒也不具有物理意义，这一操作无法保证原需要求解的本构方程对流项的守恒性。

上述涉及的对数构象张量界面插值问题、非稳态项插值问题与对流项守恒性问题似乎都难以直接地得到很好解决。例如本研究提出的界面对数构象张量物理域插值方法只能在简单的CD格式中实施，在对LCR方法稳定性有重要帮助的KT格式、MINMOD格式下实施会产生难以预料的结果。因此，我们未来考虑将本构方程分解为特征值方程与特征矩阵方程。求解特征值方程时使用LCR方法，但是非稳态项在物理域计算然后作为整体转换到对数域，从而解决上述三个问题。值得注意的是，三个特征值方程为互不相关的标量方程，可以准确而直接地实现上述转换，即直接除以各自特征值本身。

作者简介：张文华，本硕博毕业于中国石油大学（北京），微信号390819095，邮箱[e-mail]，研究方向为黏弹性流体湍流减阻、计算流体力学与计算传热学 ，目前就职于中山大学中法核工程与技术学院。