具有间断解的PTT型黏弹性流体热对流数值模拟*

张栏， 蒋子超， 陈玉惠， 张仪， 杨耿超， 姚清河

1.中山大学航空航天学院，广东 深圳 518107

2.Robotics and Mechatronics， University of Twente， Enschede 7522 NB， the Netherlands

热对流是由温度差驱动的流体运动，是自然界中普遍而重要的现象，广泛存在于海洋、大气、地核以及恒星和行星的内部。Rayleigh-Bénard（RB）热对流是最常见、最典型的热对流模型之一。近年来，国内外学者针对牛顿流体的RB 热对流进行了广泛研究（Bodenschatz et al.，2000；石峯等，2008；周全等，2012；贺啸秋等，2022），但对黏弹性流体的RB 热对流研究仍较少。事实上，研究黏弹性流体的RB 热对流对于化学工业和地质研究都有重要意义（Hayat et al.，2008；Ogawa，2008）。相对牛顿流体而言，黏弹性流体包含了复杂的流变动力学特性，其实验模拟和解析求解都存在较大困难。近年来，随着计算机性能和计算流体力学数值方法的发展，黏弹性流体数值模拟的精度与分辨率都取得了突破性提升，对其进行直接数值模拟也逐渐成为了该领域重要的研究手段（Abu-Ramadan et al.，2003；Goswami et al.，2021；Kuron et al.，2021；Tseng，2021）。

黏弹性流体的本构模型选择对于其仿真有着关键影响，当前被广泛采用的本构模型主要有Oldroyd-B 模型（Oldroyd，1950；李勇等，2019）、FENE-P 模型（Bird et al.，1995）、Phan-Thien-Tanner（PTT）模型（Chen et al.，2021，2022）等。其中，PTT 模型由Phan-Thien 和Tanner 根据Lodge-Yamamoto 聚合物网络理论推导得到（Thien et al.，1977；Phan-Thien，1978），相较于Oldroyd-B 模型，PTT本构模型更适合描述聚合物熔体和浓缩溶液的流变特性。通过调整模型材料参数，PTT模型可模拟多种高分子聚合物的流变性质（Quinzani et al.，

1994；Azaiez et al.，1996；Schoonen et al.，1998；Shin et al.，2007）。对于挤出胀大等黏弹性流体效应，基于PTT 模型的仿真结果具有较好的预测能力（Edeleva et al.，2021）。

近年来已有基于PTT 模型的黏弹性流体的热对流及其他流动仿真，特别是低Rayleigh 数下RB热对流的数值模拟研究。例如Park et al.（2004，2009，2018）使用基于逐格反演方法（grid-by-grid inversion method）的有限体积法，Zheng et al.（2022）使用基于高阶迎风中心格式的有限差分法（FDM，finite difference method）均实现了PTT 型黏弹性流体RB 热对流的数值模拟，对其流动特性进行了详细探讨。目前研究所采用的数值方法普遍难以处理间断现象，需要考虑高数值稳定性的新型离散格式。

方程非线性项引发的间断解问题一直是计算流体力学中的关键问题，围绕方程非线性项的离散，已有许多重要研究方法，总变差不增（TVD，total variation diminishing）格式为其中的关键成果之一。TVD 格式最早由Harten（1983）提出并被应用于可压缩流动的FDM 求解，由于其能够准确地捕捉激波的位置、在高分辨率下仍可以维持对间断解的数值稳定性，从而逐渐成为了主流离散格式。近年来TVD 格式被广泛应用于磁流体动力学（Yuan，2019）、地热储层工程（Oldenburg et al.，2000）和溃坝水力学（刘玉玲等，2010）等领域内，在其基础上构造的其他高精度无波动格式也取得了诸多应用成果（Čada et al.，2009；Kemm，2011；Dubey，2013）。但是，将TVD 格式应用于黏弹性流体RB对流数值模拟的研究仍较少。

从间断解充分发展的黏弹性流体RB 热对流具备不同流态的时间信息和发展过程的流态变换，研究该问题对阐明PTT 型黏弹性流体RB 热对流的流动机理具有重要意义。为解决PTT 本构模型在RB 对流中可能会引起间断和数值振荡问题，本研究将基于TVD 离散格式和PTT 本构模型，开发一种黏弹性流体RB 热对流FDM 模拟方法。该方法能够应用间断初始条件进行长时间模拟，并通过与稳定周期解的对比验证其有效性。

1 数值模型及算法实现

1.1 控制方程

采取布辛涅司克近似模拟该对流问题，控制方程（Zheng et al.，2022）为

其中u，p，T，τ分别为黏弹性流体的速度、压力、温度和偏应力张量。基准温度T0=(T1+T2)/2，T1为上边界温度，T2为下边界温度，ρ0为T=T0时黏弹性流体的密度。常数α为热膨胀系数，g为重力加速度，k为热传导系数，Cp为比热容，ej为重力加速度方向上的单位向量。

偏应力张量τ为二阶对称张量，通常可以分解为牛顿溶剂贡献和聚合物贡献

其中τs= 2μsD为黏性偏应力张量，μs为溶剂黏度，应变率张量D=.τp为弹性偏应力张量，对于PTT本构模型有

其中λ为弛豫时间，μp为聚合物黏度，ϵ为分子应变率，ξ为分子滑移率。

将控制方程中各量进行无量纲化，

其中H为特征尺度，ΔT为最大（边界）温差，Uc=为特征速度，κ为热扩散系数，ν为运动学黏度。

无量纲化的控制方程为

其中Ra=为瑞利数，We=κλ/H2为魏森贝格数，Pr=(μ0Cp)/k为普朗特数，均为黏弹性流体RB 热对流问题中重要无量纲参数。β=μs/μ0为比黏度，即溶剂黏度与总黏度的比值，总黏度μ0=μs+μp，Ma=为马赫数，表示横波速度与流体特征速度的比值。

本文采用有限差分法对上述方程进行数值模拟，其中温度场由显格式求解，压力场基于SIMPLE算法求解。

1.2 适用于间断解的离散方法

首先将方程（2）～（3）分离源项拟线性化并转化成守恒形式得到

其中W=[u，v，τ11，τ12，τ22]T，u和v分别为x和y方向上的速度分量，τij(i，j= 1，2)为弹性偏应力的应力分量，流通量矢量定义为

相应地，源项

应用局部特征方法，方程(5)可以表示为

其中雅可比矩阵

A2具有类似的形式。

时间上对时间的偏导数项采取一阶向前差分格式离散，其他项均显式处理。方程（1）～（4）在时间上的半离散格式为

空间上的离散采用均匀网格，使用二阶中心差分格式离散黏性项、能量方程的扩散项，使用二阶迎风TVD差分格式（Harten，1983）离散拟线性项

格式黏性函数Q(z)定义为

本文取ζ=[0.05，0.5].

定义

本文采用minmod函数作为限制器

其中

2 数值验证

为了验证本文的数值方法的离散精度与收敛性，本文将基于理论解析解对数值结果进行精度校验，考虑具有额外源项f1，f2的控制方程

其中额外源项f1=[fu，fv，fτ11，fτ12，fτ22]T，f2=fT，且

可证明，在宽高比为1∶1的方腔内，方程(6)～(8)的解析解

其中ω=kπ/L，η=10-5.选取参数β= 0.2，We= 0.1，Ra= 1 480，Pr= 0.7，ϵ= 0.1，ξ= 0.05，取时间步长Δt= 10-5，在不同空间步长Δx下，本文所采用的求解算法所得的数值解与式(9)的解析解间2 范数误差（L2-error）如图1所示。

图1中，Δx选取了1/1 024 ～ 1/32间对数坐标上的等距点，其对应的2 范数误差与Δx近似成幂函数关系，其指数约为1.71.由此论证了本文所采用的TVD格式在空间上的收敛性与近似2阶精度。

3 结果分析

本文RB 热对流几何模型采用宽高比为2∶1 的充满黏弹性流体的方腔（见图2）。

图2 RB热对流几何模型与温度边界设置Fig.2 Geometric model of RB convection

无量纲化后该PTT 型黏弹性流体RB 热对流的计算域为(x，y)：[0，2]×[0，1].速度边界条件为四壁均为无滑移边界条件；应力边界条件为诺伊曼边界条件；温度边界条件为上下壁面为等温边界条件，左右壁面为绝热边界条件，即

速度、应力初始条件均为0，初始温度场为间断场：上边界和下边界分别为0.5 和-0.5，除边界外均为T0.

各个常数取值分别为β= 0.2，We= 0.1，Ra=1 480，Pr= 0.7，ϵ= 0.1，ξ= 0.05.构建空间步长为1/64的均匀网格，取Δt= 10-3，对该模型进行时长550的数值模拟，其动能（Ek）发展如图3所示。

图3 （a）全局动能随时间的变化情况；（b）动能波动周期对比及单个周期内的5个代表性时间步（t1～t5）Fig.3 (a) Time evolution of global kinetic energy,（b） Comparison of kinetic energy fluctuation periods and five representative time steps （t1～t5）within a single period

由图3 可见，PTT 黏弹性流体的RB 热对流在充分发展后呈现出周期性流动的特性，动能随着时间的波动周期为10.4，与Zheng et al.（2022）中由稳定周期解作为初始场的计算结果10.2 接近一致。由动能的量级估计，流动开始发展的时间为232.3且358.7 后达到稳定周期解。为进一步阐释单周期内流动的性质，取相邻动能最大值点内的5个代表性时间步，其温度与速度场分布如图4所示。

图4 一个周期内的温度场和速度场分布（由上至下分别对应t1，t2，t3，t4，t5）Fig.4 The distribution of velocity (a) and temperature (b) at the five time points(t1 to t5 from top to bottom) in one period

图4 中，t1= 494，t5= 504.4 时，全局动能最大；t2= 499.2 时，方腔中心线底壁附近有新涡出现；t3= 499.5 时，全局动能最小；t4= 500.9 时，温度场分布近似纯热传导。由图4（a）可见，在达到全局动能最大值(t=t1)后，流动沿着方腔中心线向上流动；至t=t2时，在方腔垂直中心线附近的底壁附近开始出现2 个小涡，然后迅速发展至t3时形成的四涡结构；时间在t3到t4之间涡逐渐合并，并在t4时再次过渡为两涡结构。从方腔垂直中心线附近底壁出现小涡开始，流场的流向沿着方腔中心线向上流动，到t3时四涡结构形成，整个流场转变到沿着方腔中心线向下流动。由图4（b）可见，在t1到t3之间，温度场基本没有发生变化，但速度场却在持续减弱，最大速度从t1的0.08 变化到t3的0.006。在t3～t4时间内，温度场发生明显变化。在t4时温度分布几乎可以近似在y方向上呈线性分布，接近于无对流的纯热传导的稳态温度分布。从t4到t5，全场速度逐渐增大并在t5处达到极值，温度场也随之偏离t=t3时的分布。

区别于牛顿流体，黏弹性流体热对流除温差驱动力外还存在着非均匀分布的内应力。在该数值模拟中，应力场的仿真结果如图5所示。

图5 一个周期内的弹性偏应力场分布（由上至下分别对应t1，t2，t3，t4，t5）Fig.5 The distribution of extra-stress τ11, τ12 and τ22 at the five time points(t1 to t5 from top to bottom) in one period

由图5 可知，在t=t1后，速度梯度极值不断增大，同时弹性偏应力分量随之增大。在t2时，τ11和τ22的最大值出现在(x，y) =(1，0.9) 和(x，y) =(1，0.3)附近，而从t1向t2的过程中，τij的散度增大，使得u和v由于∂τ11/∂x和∂τ22/∂y而增大。从t2到t3阶段τij达到极大值后向减小开始变化，但减小的幅度很轻微；在t3到t4阶段τij减小幅度很大，并于t4时τij达到极小值。此阶段速度场的增大加剧了τij的衰减。从t3到t4阶段，τ11和τ22的正负沿垂直中心线变化。在t3时，τ11在(x，y) =(1，0.9)处为正，在(x，y) =(1，0.3)处为负，而在t5时，τ11在相同位置改变正负，在(x，y) =(1，0.9)处为负，在(x，y) =(1，0.3)处为正。

在下一个全局动能周期，也是速度场的下半周期，速度梯度由零开始逐渐向极值增大，速度梯度的增大使τij增大，全局动能最小时在垂直中心线附近的顶壁附近出现小涡，迅速扩大并侵入腔体快速转化成四涡结构，流动方向发生逆转，温度场的对流分布转变为近似在y方向上呈线性分布；τij继续衰减，τ11和τ22沿方腔垂直中心线变化，速度再次放大达到极值，温度场再次变为对流分布。

上述流动现象表明，速度场的减弱-放大循环与弹性偏应力场的衰减-增大循环之间存在相位差，相位差的存在导致了PTT 型黏弹性流体RB 热对流随时间变化的流动状态，并且在速度场的完整周期里出现反向对流现象。

4 结 论

本文采用了基于TVD 差分格式的FDM 模拟了基于PTT 型黏弹性流体的RB 热对流问题，对非线性黏弹性流体在封闭方腔内部的对流过程进行了研究，得到如下结论：

1） 本文的数值结果与解析解与文献相符，且具有二阶空间精度。

2） 数值模拟结果表明，基于TVD 格式的求解方法能够有效地模拟黏弹性流体RB 热对流问题，并得到流动从间断初始条件开始发展到稳定周期解的完整过程的瞬态信息。

3） PTT 型黏弹性流体RB 热对流的基本流动现象显示，在充分发展后会呈现出随时间变化的振荡对流特征，并具有特定流动模式以及流动反转的流动特性。