基于THINC/QQ格式的兩相界面流動(dòng)數(shù)值模擬

趙海洋，明平劍，張文平，齊文亮

(哈爾濱工程大學(xué) 動(dòng)力與能源工程學(xué)院，黑龍江 哈爾濱 150001)

兩相界面流動(dòng)現(xiàn)象隨處可見。在實(shí)際工程計(jì)算中，海洋平臺(tái)受力分析、波浪水池、氣泡上升及毛細(xì)管中氣泡流動(dòng)問題，都需要借助于精確實(shí)用的兩相流模擬技術(shù)才能有效解決問題[1-2]。眾多學(xué)者提出并發(fā)展了一系列成熟的可用于商業(yè)計(jì)算的兩相界面捕捉方法，其中體積分?jǐn)?shù)類方法(volume of fluid，VOF)因其嚴(yán)格的質(zhì)量守恒特性得到了廣泛應(yīng)用。Xiao等[3]提出了采用雙曲正切函數(shù)重構(gòu)網(wǎng)格內(nèi)自由界面(tangent of hyperbola interface capturing method，THINC)方法，該方法看作是幾何重構(gòu)/代數(shù)混合型的體積分?jǐn)?shù)方法，通過確定單元內(nèi)表征兩相分布的雙曲正切函數(shù)顯式確定網(wǎng)格單元內(nèi)自由界面位置，通過對(duì)雙曲正切函數(shù)在網(wǎng)格間界面上積分平均得到數(shù)值流通量，其間并沒有復(fù)雜的幾何重構(gòu)。通過參數(shù)β，THINC格式能夠很好的控制自由界面的厚度在2、3層網(wǎng)格，而不會(huì)帶來很強(qiáng)的數(shù)值擴(kuò)散。Xiao等[4]將該方法從一維擴(kuò)展至多維，由結(jié)構(gòu)化網(wǎng)格到非結(jié)構(gòu)化網(wǎng)格，網(wǎng)格內(nèi)自由界面的形狀也由一開始的傾斜線段(平面)發(fā)展到THINC/QQ中考慮界面曲率的曲線(曲面)，已發(fā)展成一種適用于多種網(wǎng)格類型的能精確模擬兩相界面的成熟的界面捕捉方法[5-9]。

本文基于守恒形式不可壓同位網(wǎng)格有限體積法離散N-S方程，采用THINC/QQ格式捕捉兩相間自由界面，對(duì)互不相溶兩相自由面問題進(jìn)行數(shù)值模擬。利用有限體積法框架下的方程離散和求解，通過典型算例驗(yàn)證數(shù)值方法在不同類型網(wǎng)格下的實(shí)施精度并得出結(jié)論。

1 兩相流數(shù)學(xué)模型

1.1 控制方程

假定兩相流體為不可壓、互不相溶的粘性流體，兩相之間的交界面被認(rèn)為是物性參數(shù)的間斷，界面單元中兩相流體共用速度和壓力。連續(xù)性方程、動(dòng)量方程以及體積分?jǐn)?shù)方程可表示為：

(1)

式中：u和p分別是速度和壓力；t為時(shí)間；g為重力加速度；Fσ表示兩相流體間表面張力；ρ和μ分別是密度和動(dòng)力粘性系數(shù)：

(2)

式中：下標(biāo)1和2分別指代2種流體；φ為單元中流體1所占的體積分?jǐn)?shù)。

根據(jù)連續(xù)表面力模型，兩相流體間表面張力為：

(3)

式中：σ為表面張力系數(shù)；κ為兩相界面曲率，可表示為：

(4)

采用守恒形式有限體積法在非結(jié)構(gòu)同位網(wǎng)格上對(duì)方程(1)進(jìn)行離散，采用一階歐拉隱格式處理時(shí)間項(xiàng)，采用二階迎風(fēng)格式離散對(duì)流項(xiàng)，對(duì)擴(kuò)散項(xiàng)采用中心差分格式離散，對(duì)連續(xù)性方程和動(dòng)量方程中的壓力速度耦合問題，采用SIMPLE算法進(jìn)行處理。采用THINC/QQ格式對(duì)體積分?jǐn)?shù)方程進(jìn)行求解。

1.2 THINC/QQ格式的數(shù)值實(shí)施方案

THINC格式借助于雙曲正切函數(shù)H(x)表征兩相分布，對(duì)于多維情況，界面單元中的分段雙曲正切函數(shù)可表示為：

(5)

式中：β是控制界面厚度的參數(shù)；Pi(x,y,z)+di=0表示網(wǎng)格單元中界面位置方程。網(wǎng)格單元的體積分?jǐn)?shù)是H函數(shù)在整個(gè)控制單元內(nèi)的積分平均：

(6)

為了計(jì)算簡化，將整個(gè)界面重構(gòu)過程映射到等參坐標(biāo)系(ξ,η,ζ)下進(jìn)行計(jì)算。表征單元中兩相界面位置的表達(dá)式Pi(ξ,η,ζ)采用二次多項(xiàng)式表示：

Pi(ξ,η,ζ)=a200ξ2+a020η2+a002ζ2+a110ξη+

a011ηζ+a101ξζ+a100ξ+a010η+a001ζ

(7)

式中：astr(s,t,r=0,1,2，s+t+r<=2)是確定界面形狀的參數(shù)，可以通過界面法線、曲率顯示表示。在確定astr后，可以通過求解以下方程計(jì)算得到di：

(8)

假定在第i個(gè)單元中設(shè)置G個(gè)高斯積分點(diǎn)，每個(gè)點(diǎn)對(duì)應(yīng)坐標(biāo)為ξig=(ξig,ηig,ζig)，則等式(8)可以轉(zhuǎn)換為：

(9)

式中ωig為第g個(gè)高斯積分點(diǎn)占的權(quán)重，且滿足：

(10)

根據(jù)正切函數(shù)恒等變形公式，得到：

tanh(β(Pi(ξig)+di))=

(11)

將式(10)和(11)代入到方程(9)中，可以得到：

(12)

方程(12)可以看作關(guān)于未知數(shù)tanh(βdi)的一元非線性方程，可以采用牛頓-拉夫遜法迭代求解，進(jìn)而得到di的值。

結(jié)合式(6)，體積分?jǐn)?shù)方程(1)可離散為：

(13)

式中：vnij=v·nij是面Γij的法向速度；下標(biāo)iup表示Γij的上游迎風(fēng)單元，考慮到流動(dòng)的方向性，采用上游迎風(fēng)單元中的體積分?jǐn)?shù)分布確定通過單元界面上的體積流量。一般建議采用三階龍格-庫塔方法(total variation diminishing，TVD)對(duì)方程(13)進(jìn)行顯示迭代更新。式(13)中右端的積分項(xiàng)用高斯積分進(jìn)行計(jì)算：

(14)

式中：G是高斯點(diǎn)數(shù)目；ωg是權(quán)重系數(shù)；(ξg,ηg,ζg)是局部坐標(biāo)系下Γij上的高斯點(diǎn)坐標(biāo)。

2 數(shù)值結(jié)果及分析

本節(jié)通過3個(gè)典型的自由面流動(dòng)算例測試有限體積法框架下THINC/QQ格式的性能。給定剪切流場，測試兩相界面在剪切流場中發(fā)生大變形時(shí)，當(dāng)前算法捕捉界面的精度。分別對(duì)重力、表面張力作用下的二維和三維氣泡上浮問題進(jìn)行數(shù)值模擬，說明當(dāng)前數(shù)值算法的適用性。

2.1 二維剪切流動(dòng)

考慮半徑0.15的圓在[0,1]×[0,1]計(jì)算域中發(fā)生剪切變形，初始圓心坐標(biāo)為(0.5,0.75)，速度場為：

(15)

計(jì)算周期T為 8.0 s。界面在剪切速度場變形為螺旋狀，在T/2時(shí)刻變形最大，在T時(shí)刻回到初始位置。采用非結(jié)構(gòu)化混合網(wǎng)格對(duì)當(dāng)前問題進(jìn)行測試，如圖1所示，網(wǎng)格數(shù)分別為2 564、10 025、40 217。

圖1 二維剪切流動(dòng)計(jì)算網(wǎng)格示意Fig.1 Computational mesh for the 2D shear flow test

采用當(dāng)前THINC/QQ格式對(duì)該問題進(jìn)行數(shù)值模擬，取界面厚度控制參數(shù)β=3.6，計(jì)算時(shí)間步長5.0×10-4s。得到不同網(wǎng)格數(shù)下數(shù)值誤差，如表1所示。其中數(shù)值誤差定義為：

(16)

從表1可以看出，隨著網(wǎng)格加密，計(jì)算誤差顯著降低，計(jì)算結(jié)果更加精確。

表1 不同網(wǎng)格下數(shù)值誤差Table 1 Numerical errors under different meshes

圖2給出了40 217網(wǎng)格下T/2時(shí)刻和T時(shí)刻的兩相界面形狀，可以看到，在T/2時(shí)刻，界面變成螺旋形，尾部有破碎的小液滴，這是因?yàn)榫W(wǎng)格尺度較尾部界面厚度尺寸過大。此外盡管經(jīng)過了很大變形，T時(shí)刻兩相界面仍能夠回復(fù)到最初圓形，數(shù)值解和解析解吻合良好。圖3為T時(shí)刻的體積分?jǐn)?shù)等值線放大云圖，圖中選取體積分?jǐn)?shù)為0.001、0.5和0.999的等值線，可以看到界面厚度均勻分布，長時(shí)間計(jì)算后THINC/QQ格式仍能夠控制界面厚度在3到4層網(wǎng)格之間。

圖2 T/2和T時(shí)刻界面形狀Fig.2 The interface profile plotted by at T/2 and T

2.2 二維氣泡上浮

本節(jié)考慮在重力以及表面張力作用下的二維氣泡上浮問題，驗(yàn)證當(dāng)前數(shù)值算法的計(jì)算精度和魯棒性。Hysing等[10]對(duì)2個(gè)典型二維氣泡上升算例采用3種不同數(shù)值算法進(jìn)行了定量計(jì)算分析，得到了典型時(shí)刻氣泡形狀、氣泡上浮位置以及氣泡上浮速度的時(shí)間變化曲線。其數(shù)值結(jié)果被廣泛用于程序驗(yàn)證。計(jì)算模型示意圖如圖4所示。

直徑d=0.5的氣泡位于[0,1]×[0,2]的矩形計(jì)算域中，上下邊界設(shè)置為不可滑移固壁，左右邊界為可滑移固壁。2個(gè)算例的物性參數(shù)設(shè)置如表2所示。在計(jì)算域內(nèi)劃分三角形網(wǎng)格，網(wǎng)格數(shù)為44 760。

圖3 T時(shí)刻體積分?jǐn)?shù)為0.001，0.5，0.999的等值線圖Fig.3 Counter maps for volume fractions values equaling 0.001，0.5 and 0.999 at T

圖4 二維氣泡上浮計(jì)算模型示意Fig.4 Computational setup for modeling the 2D rising bubble

圖5和圖6分別給出了t=0.0,1.0,2.0,3.0 s時(shí)刻2種物性參數(shù)設(shè)置下上浮氣泡的形狀，可以看到氣泡在重力和表面張力作用下發(fā)生變形，在算例2中兩相密度比和粘性系數(shù)比更大，而且兩相間表面張力系數(shù)更小，導(dǎo)致氣泡在上浮過程中逐漸變化為凹形，且逐步拉長拉細(xì)，并最終破碎成小氣泡。

表2 物性參數(shù)設(shè)置Table 2 Physical parameters

圖5 算例1不同時(shí)刻氣泡形狀Fig.5 Bubble shapes for case 1

圖6 算例2不同時(shí)刻氣泡形狀Fig.6 Bubble shapes for case 2

圖7和圖8分別給出了2個(gè)算例氣泡上浮位移和速度隨時(shí)間的變化曲線，并與Hysing等[10]的數(shù)值結(jié)果進(jìn)行比較。可以看到當(dāng)前數(shù)值算法在三角形網(wǎng)格下得到的氣泡位移和上浮速度時(shí)歷曲線與文獻(xiàn)[10]在結(jié)構(gòu)化網(wǎng)格得到的計(jì)算結(jié)果吻合良好。

圖7 算例1氣泡位移和速度隨時(shí)間變化曲線Fig.7 Curves of position and velocity against time for case 1

2.3 三維氣泡融合

本節(jié)考慮三維情況下2個(gè)氣泡在表面張力作用下發(fā)生融合。所用到的無量綱參數(shù)埃奧特沃斯數(shù)Eo、莫頓數(shù)Mo以及雷諾數(shù)Re分別定義為：

(17)

式中：U為氣泡上浮速度；考慮直徑為d、球心垂直方向相距1.5d的2個(gè)同軸氣泡，底部氣泡距離計(jì)算域底部為d，分別考慮2個(gè)同軸氣泡和2個(gè)水平方向相距0.8d的氣泡的上浮融合情況。計(jì)算模型示意圖如圖9所示。計(jì)算域大小為4d×4d×8d，采用規(guī)整的六面體網(wǎng)格，網(wǎng)格尺寸為h=1/20d。上下邊界(z=0，z=8d)采用無滑移固壁邊界，四周為滑移固壁邊界，初始時(shí)刻氣泡和周圍流體均處于靜止?fàn)顟B(tài)。兩相流體物性參數(shù)滿足Eo=16，Mo=2.0×10-4，兩相密度比和動(dòng)力粘度比均為100。計(jì)算過程中時(shí)間步長滿足庫朗數(shù)小于0.25。

圖8 算例2氣泡位移和速度隨時(shí)間變化曲線Fig.8 Curves of position and velocity against time for case 2

圖9 三維兩氣泡融合計(jì)算模型示意Fig.9 Computational setup for modeling the 3D merging bubbles

圖10給出了同軸2個(gè)氣泡和水平方向相距0.8d的兩傾斜氣泡上浮過程中的雷諾數(shù)時(shí)歷曲線，計(jì)算結(jié)果與Balcázar等[11-12]采用Level Set方法和耦合VOF/LS方法得到的數(shù)值結(jié)果進(jìn)行了對(duì)比，可以看到不管兩氣泡初始相對(duì)位置如何，本文采用的THINC/QQ格式能夠準(zhǔn)確捕捉兩氣泡上浮過程的速度特性。

圖10 氣泡上浮雷諾數(shù)時(shí)歷曲線Fig.10 Time evolution of the Reynold number

圖11和12給出了2同軸氣泡和水平方向相距0.8d的2傾斜氣泡的融合過程，結(jié)合圖10所示的雷諾數(shù)時(shí)歷曲線，可以看到氣泡上升初期，氣泡上升速度迅速增大，兩氣泡形狀變化相近，隨著時(shí)間推移，頂部氣泡逐漸趨向于扁平狀，底部氣泡進(jìn)入頂部氣泡的尾流區(qū)，導(dǎo)致底部氣泡逐漸變成子彈形狀，底部氣泡速度增大，追上頂部氣泡并與之融合，隨后融合氣泡繼續(xù)上浮，但上浮速度在減小。對(duì)初始時(shí)刻相對(duì)傾斜的兩上浮氣泡，其流場分布也不是對(duì)稱的，最終融合的大氣泡的主軸方向與z軸方向存在夾角。

圖11 同軸情況下不同時(shí)刻2氣泡的形狀Fig.11 Snapshots at different times of the coaxial coalescence

圖12 傾斜情況下不同時(shí)刻2氣泡的形狀Fig.12 Snapshots at different times of the oblique coalescence

3 結(jié)論

1)當(dāng)前算法能準(zhǔn)確模擬非結(jié)構(gòu)化網(wǎng)格下兩相界面的變形、分離以及融合，數(shù)值結(jié)果與相應(yīng)的解析解和文獻(xiàn)中數(shù)值解吻合良好，說明當(dāng)前兩相流模擬程序的準(zhǔn)確性。

2)氣相物性參數(shù)和表面張力系數(shù)直接影響上浮過程中氣泡形狀，減小氣相密度、粘性系數(shù)及表面張力系數(shù)更容易發(fā)生氣泡破碎。

在一定條件下，兩垂直方向存在指定間距的氣泡會(huì)發(fā)生融合，對(duì)影響氣泡融合的參數(shù)有待進(jìn)一步研究。