張欄， 蔣子超， 陳玉惠， 張儀， 楊耿超， 姚清河

1.中山大學(xué)航空航天學(xué)院，廣東 深圳 518107

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

熱對流是由溫度差驅(qū)動的流體運動，是自然界中普遍而重要的現(xiàn)象，廣泛存在于海洋、大氣、地核以及恒星和行星的內(nèi)部。Rayleigh-Bénard（RB）熱對流是最常見、最典型的熱對流模型之一。近年來，國內(nèi)外學(xué)者針對牛頓流體的RB 熱對流進(jìn)行了廣泛研究（Bodenschatz et al.，2000；石峯等，2008；周全等，2012；賀嘯秋等，2022），但對黏彈性流體的RB 熱對流研究仍較少。事實上，研究黏彈性流體的RB 熱對流對于化學(xué)工業(yè)和地質(zhì)研究都有重要意義（Hayat et al.，2008；Ogawa，2008）。相對牛頓流體而言，黏彈性流體包含了復(fù)雜的流變動力學(xué)特性，其實驗?zāi)M和解析求解都存在較大困難。近年來，隨著計算機性能和計算流體力學(xué)數(shù)值方法的發(fā)展，黏彈性流體數(shù)值模擬的精度與分辨率都取得了突破性提升，對其進(jìn)行直接數(shù)值模擬也逐漸成為了該領(lǐng)域重要的研究手段（Abu-Ramadan et al.，2003；Goswami et al.，2021；Kuron et al.，2021；Tseng，2021）。

黏彈性流體的本構(gòu)模型選擇對于其仿真有著關(guān)鍵影響，當(dāng)前被廣泛采用的本構(gòu)模型主要有Oldroyd-B 模型（Oldroyd，1950；李勇等，2019）、FENE-P 模型（Bird et al.，1995）、Phan-Thien-Tanner（PTT）模型（Chen et al.，2021，2022）等。其中，PTT 模型由Phan-Thien 和Tanner 根據(jù)Lodge-Yamamoto 聚合物網(wǎng)絡(luò)理論推導(dǎo)得到（Thien et al.，1977；Phan-Thien，1978），相較于Oldroyd-B 模型，PTT本構(gòu)模型更適合描述聚合物熔體和濃縮溶液的流變特性。通過調(diào)整模型材料參數(shù)，PTT模型可模擬多種高分子聚合物的流變性質(zhì)（Quinzani et al.，

1994；Azaiez et al.，1996；Schoonen et al.，1998；Shin et al.，2007）。對于擠出脹大等黏彈性流體效應(yīng)，基于PTT 模型的仿真結(jié)果具有較好的預(yù)測能力（Edeleva et al.，2021）。

近年來已有基于PTT 模型的黏彈性流體的熱對流及其他流動仿真，特別是低Rayleigh 數(shù)下RB熱對流的數(shù)值模擬研究。例如Park et al.（2004，2009，2018）使用基于逐格反演方法（grid-by-grid inversion method）的有限體積法，Zheng et al.（2022）使用基于高階迎風(fēng)中心格式的有限差分法（FDM，finite difference method）均實現(xiàn)了PTT 型黏彈性流體RB 熱對流的數(shù)值模擬，對其流動特性進(jìn)行了詳細(xì)探討。目前研究所采用的數(shù)值方法普遍難以處理間斷現(xiàn)象，需要考慮高數(shù)值穩(wěn)定性的新型離散格式。

方程非線性項引發(fā)的間斷解問題一直是計算流體力學(xué)中的關(guān)鍵問題，圍繞方程非線性項的離散，已有許多重要研究方法，總變差不增（TVD，total variation diminishing）格式為其中的關(guān)鍵成果之一。TVD 格式最早由Harten（1983）提出并被應(yīng)用于可壓縮流動的FDM 求解，由于其能夠準(zhǔn)確地捕捉激波的位置、在高分辨率下仍可以維持對間斷解的數(shù)值穩(wěn)定性，從而逐漸成為了主流離散格式。近年來TVD 格式被廣泛應(yīng)用于磁流體動力學(xué)（Yuan，2019）、地?zé)醿庸こ蹋∣ldenburg et al.，2000）和潰壩水力學(xué)（劉玉玲等，2010）等領(lǐng)域內(nèi)，在其基礎(chǔ)上構(gòu)造的其他高精度無波動格式也取得了諸多應(yīng)用成果（?ada et al.，2009；Kemm，2011；Dubey，2013）。但是，將TVD 格式應(yīng)用于黏彈性流體RB對流數(shù)值模擬的研究仍較少。

從間斷解充分發(fā)展的黏彈性流體RB 熱對流具備不同流態(tài)的時間信息和發(fā)展過程的流態(tài)變換，研究該問題對闡明PTT 型黏彈性流體RB 熱對流的流動機理具有重要意義。為解決PTT 本構(gòu)模型在RB 對流中可能會引起間斷和數(shù)值振蕩問題，本研究將基于TVD 離散格式和PTT 本構(gòu)模型，開發(fā)一種黏彈性流體RB 熱對流FDM 模擬方法。該方法能夠應(yīng)用間斷初始條件進(jìn)行長時間模擬，并通過與穩(wěn)定周期解的對比驗證其有效性。

1 數(shù)值模型及算法實現(xiàn)

1.1 控制方程

采取布辛涅司克近似模擬該對流問題，控制方程（Zheng et al.，2022）為

其中u，p，T，τ分別為黏彈性流體的速度、壓力、溫度和偏應(yīng)力張量?；鶞?zhǔn)溫度T0=(T1+T2)/2，T1為上邊界溫度，T2為下邊界溫度，ρ0為T=T0時黏彈性流體的密度。常數(shù)α為熱膨脹系數(shù)，g為重力加速度，k為熱傳導(dǎo)系數(shù)，Cp為比熱容，ej為重力加速度方向上的單位向量。

偏應(yīng)力張量τ為二階對稱張量，通?？梢苑纸鉃榕ｎD溶劑貢獻(xiàn)和聚合物貢獻(xiàn)

其中τs= 2μsD為黏性偏應(yīng)力張量，μs為溶劑黏度，應(yīng)變率張量D=.τp為彈性偏應(yīng)力張量，對于PTT本構(gòu)模型有

其中λ為弛豫時間，μp為聚合物黏度，?為分子應(yīng)變率，ξ為分子滑移率。

將控制方程中各量進(jìn)行無量綱化，

其中H為特征尺度，ΔT為最大（邊界）溫差，Uc=為特征速度，κ為熱擴散系數(shù)，ν為運動學(xué)黏度。

無量綱化的控制方程為

其中Ra=為瑞利數(shù)，We=κλ/H2為魏森貝格數(shù)，Pr=(μ0Cp)/k為普朗特數(shù)，均為黏彈性流體RB 熱對流問題中重要無量綱參數(shù)。β=μs/μ0為比黏度，即溶劑黏度與總黏度的比值，總黏度μ0=μs+μp，Ma=為馬赫數(shù)，表示橫波速度與流體特征速度的比值。

本文采用有限差分法對上述方程進(jìn)行數(shù)值模擬，其中溫度場由顯格式求解，壓力場基于SIMPLE算法求解。

1.2 適用于間斷解的離散方法

首先將方程（2）～（3）分離源項擬線性化并轉(zhuǎn)化成守恒形式得到

其中W=[u，v，τ11，τ12，τ22]T，u和v分別為x和y方向上的速度分量，τij(i，j= 1，2)為彈性偏應(yīng)力的應(yīng)力分量，流通量矢量定義為

相應(yīng)地，源項

應(yīng)用局部特征方法，方程(5)可以表示為

其中雅可比矩陣

A2具有類似的形式。

時間上對時間的偏導(dǎo)數(shù)項采取一階向前差分格式離散，其他項均顯式處理。方程（1）～（4）在時間上的半離散格式為

空間上的離散采用均勻網(wǎng)格，使用二階中心差分格式離散黏性項、能量方程的擴散項，使用二階迎風(fēng)TVD差分格式（Harten，1983）離散擬線性項

格式黏性函數(shù)Q(z)定義為

本文取ζ=[0.05，0.5].

定義

本文采用minmod函數(shù)作為限制器

其中

2 數(shù)值驗證

為了驗證本文的數(shù)值方法的離散精度與收斂性，本文將基于理論解析解對數(shù)值結(jié)果進(jìn)行精度校驗，考慮具有額外源項f1，f2的控制方程

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

可證明，在寬高比為1∶1的方腔內(nèi)，方程(6)～(8)的解析解

其中ω=kπ/L，η=10-5.選取參數(shù)β= 0.2，We= 0.1，Ra= 1 480，Pr= 0.7，?= 0.1，ξ= 0.05，取時間步長Δt= 10-5，在不同空間步長Δx下，本文所采用的求解算法所得的數(shù)值解與式(9)的解析解間2 范數(shù)誤差（L2-error）如圖1所示。

圖1中，Δx選取了1/1 024 ～ 1/32間對數(shù)坐標(biāo)上的等距點，其對應(yīng)的2 范數(shù)誤差與Δx近似成冪函數(shù)關(guān)系，其指數(shù)約為1.71.由此論證了本文所采用的TVD格式在空間上的收斂性與近似2階精度。

3 結(jié)果分析

本文RB 熱對流幾何模型采用寬高比為2∶1 的充滿黏彈性流體的方腔（見圖2）。

圖2 RB熱對流幾何模型與溫度邊界設(shè)置Fig.2 Geometric model of RB convection

無量綱化后該PTT 型黏彈性流體RB 熱對流的計算域為(x，y)：[0，2]×[0，1].速度邊界條件為四壁均為無滑移邊界條件；應(yīng)力邊界條件為諾伊曼邊界條件；溫度邊界條件為上下壁面為等溫邊界條件，左右壁面為絕熱邊界條件，即

速度、應(yīng)力初始條件均為0，初始溫度場為間斷場：上邊界和下邊界分別為0.5 和-0.5，除邊界外均為T0.

各個常數(shù)取值分別為β= 0.2，We= 0.1，Ra=1 480，Pr= 0.7，?= 0.1，ξ= 0.05.構(gòu)建空間步長為1/64的均勻網(wǎng)格，取Δt= 10-3，對該模型進(jìn)行時長550的數(shù)值模擬，其動能（Ek）發(fā)展如圖3所示。

圖3 （a）全局動能隨時間的變化情況；（b）動能波動周期對比及單個周期內(nèi)的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 熱對流在充分發(fā)展后呈現(xiàn)出周期性流動的特性，動能隨著時間的波動周期為10.4，與Zheng et al.（2022）中由穩(wěn)定周期解作為初始場的計算結(jié)果10.2 接近一致。由動能的量級估計，流動開始發(fā)展的時間為232.3且358.7 后達(dá)到穩(wěn)定周期解。為進(jìn)一步闡釋單周期內(nèi)流動的性質(zhì)，取相鄰動能最大值點內(nèi)的5個代表性時間步，其溫度與速度場分布如圖4所示。

圖4 一個周期內(nèi)的溫度場和速度場分布（由上至下分別對應(yīng)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 時，方腔中心線底壁附近有新渦出現(xiàn)；t3= 499.5 時，全局動能最小；t4= 500.9 時，溫度場分布近似純熱傳導(dǎo)。由圖4（a）可見，在達(dá)到全局動能最大值(t=t1)后，流動沿著方腔中心線向上流動；至t=t2時，在方腔垂直中心線附近的底壁附近開始出現(xiàn)2 個小渦，然后迅速發(fā)展至t3時形成的四渦結(jié)構(gòu)；時間在t3到t4之間渦逐漸合并，并在t4時再次過渡為兩渦結(jié)構(gòu)。從方腔垂直中心線附近底壁出現(xiàn)小渦開始，流場的流向沿著方腔中心線向上流動，到t3時四渦結(jié)構(gòu)形成，整個流場轉(zhuǎn)變到沿著方腔中心線向下流動。由圖4（b）可見，在t1到t3之間，溫度場基本沒有發(fā)生變化，但速度場卻在持續(xù)減弱，最大速度從t1的0.08 變化到t3的0.006。在t3～t4時間內(nèi)，溫度場發(fā)生明顯變化。在t4時溫度分布幾乎可以近似在y方向上呈線性分布，接近于無對流的純熱傳導(dǎo)的穩(wěn)態(tài)溫度分布。從t4到t5，全場速度逐漸增大并在t5處達(dá)到極值，溫度場也隨之偏離t=t3時的分布。

區(qū)別于牛頓流體，黏彈性流體熱對流除溫差驅(qū)動力外還存在著非均勻分布的內(nèi)應(yīng)力。在該數(shù)值模擬中，應(yīng)力場的仿真結(jié)果如圖5所示。

圖5 一個周期內(nèi)的彈性偏應(yīng)力場分布（由上至下分別對應(yīng)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后，速度梯度極值不斷增大，同時彈性偏應(yīng)力分量隨之增大。在t2時，τ11和τ22的最大值出現(xiàn)在(x，y) =(1，0.9) 和(x，y) =(1，0.3)附近，而從t1向t2的過程中，τij的散度增大，使得u和v由于?τ11/?x和?τ22/?y而增大。從t2到t3階段τij達(dá)到極大值后向減小開始變化，但減小的幅度很輕微；在t3到t4階段τij減小幅度很大，并于t4時τij達(dá)到極小值。此階段速度場的增大加劇了τij的衰減。從t3到t4階段，τ11和τ22的正負(fù)沿垂直中心線變化。在t3時，τ11在(x，y) =(1，0.9)處為正，在(x，y) =(1，0.3)處為負(fù)，而在t5時，τ11在相同位置改變正負(fù)，在(x，y) =(1，0.9)處為負(fù)，在(x，y) =(1，0.3)處為正。

在下一個全局動能周期，也是速度場的下半周期，速度梯度由零開始逐漸向極值增大，速度梯度的增大使τij增大，全局動能最小時在垂直中心線附近的頂壁附近出現(xiàn)小渦，迅速擴大并侵入腔體快速轉(zhuǎn)化成四渦結(jié)構(gòu)，流動方向發(fā)生逆轉(zhuǎn)，溫度場的對流分布轉(zhuǎn)變?yōu)榻圃趛方向上呈線性分布；τij繼續(xù)衰減，τ11和τ22沿方腔垂直中心線變化，速度再次放大達(dá)到極值，溫度場再次變?yōu)閷α鞣植肌?/p>

上述流動現(xiàn)象表明，速度場的減弱-放大循環(huán)與彈性偏應(yīng)力場的衰減-增大循環(huán)之間存在相位差，相位差的存在導(dǎo)致了PTT 型黏彈性流體RB 熱對流隨時間變化的流動狀態(tài)，并且在速度場的完整周期里出現(xiàn)反向?qū)α鳜F(xiàn)象。

4 結(jié) 論

本文采用了基于TVD 差分格式的FDM 模擬了基于PTT 型黏彈性流體的RB 熱對流問題，對非線性黏彈性流體在封閉方腔內(nèi)部的對流過程進(jìn)行了研究，得到如下結(jié)論：

1） 本文的數(shù)值結(jié)果與解析解與文獻(xiàn)相符，且具有二階空間精度。

2） 數(shù)值模擬結(jié)果表明，基于TVD 格式的求解方法能夠有效地模擬黏彈性流體RB 熱對流問題，并得到流動從間斷初始條件開始發(fā)展到穩(wěn)定周期解的完整過程的瞬態(tài)信息。

3） PTT 型黏彈性流體RB 熱對流的基本流動現(xiàn)象顯示，在充分發(fā)展后會呈現(xiàn)出隨時間變化的振蕩對流特征，并具有特定流動模式以及流動反轉(zhuǎn)的流動特性。