適用二維切割體網格的THINC/QQ算法擴展研究

李超, 段文洋, 趙彬彬

(哈爾濱工程大學 船舶工程學院,黑龍江 哈爾濱 150001)

當采用計算流體力學(computational fluid dynamics,CFD)方法對多相流的自由界面進行模擬時,通常要對界面附近的網格進行局部加密,如波浪的二維傳播、變形和破碎等問題。存在海底突起或斜坡等不規則邊界時,采用非結構網格進行空間離散更加靈活。其中,切割體網格[1-2]既能支持懸掛點式的高效率加密方法,又能保證網格的整體正交性,因此具有很大優勢。

在自由界面捕捉方面,流體體積法(volume of fluid,VOF)由于守恒性能好、能處理自由面大變形流動,故被廣泛應用。而雙曲正切函數界面捕捉(tangent of hyperbola for interface capturing,THINC)[3]類型VOF算法,通過引入雙曲正切函數進行界面重構,有效地抑制了數值耗散。其中,THINC/QQ[4]作為一種混合幾何/代數型VOF算法,在沒有進行顯式幾何界面重構的前提下,實現了與幾何型VOF算法相當的精度,具有很好的應用前景。對于二維網格,THINC/QQ算法只能處理常規四邊形和三角形單元。然而,二維切割體網格還包含更為復雜的多邊形單元,例如在粗密網格過渡處會存在懸掛點式四邊形/三角形單元,導致THINC/QQ算法不能直接適用。

Xie等[5]提出將THINC/QQ與FVMS3(finite volume method based on merged stencil with 3rd-order reconstruction)[6]方法進行整合(本文將其記作THINC/QQ-FVMS3),能夠在任意二維多邊形單元中進行界面重構(也可用于三維問題);但是需要將目標單元分割為較多的子單元才能進行單元高斯積分,且三維的基本目標單元應是線性的[6]。Li等[7]提出了適用切割表面式非結構單元的擴展THINC/QQ算法,即THINC/QQ-SF,可用于常見形式的三維切割體網格,以及復雜三維切割體網格中的主要區域。其核心思想來自文獻[8],根據復雜單元的主要特征并借助THINC/QQ進行界面重構,當計算VOF流量時再考慮真實單元的拓撲結構。這樣就無需采用FVMS3方法,且單元高斯積分也更簡便,尤其適合采用固定網格(或彈簧網格)的多相流數值模擬。

本文進一步將THINC/QQ-SF的原理應用到二維切割體網格中,同時設計了一系列測試網格用于典型二維VOF對流和潰壩兩相流算例的驗證。 為了更好地對比,本文還采用主流CFD軟件常用的HRIC[9]算法同步進行模擬。

1 界面捕捉算法的基本原理

1.1 THINC/QQ算法基本原理

考慮物理上互不混溶的流體1和流體2,則流體1的體積分數φ遵循對流方程:

(1)

式中u是流體速度矢量。對于單元Ωi,φi為:

(2)

式中:|Ωi|表示單元Ωi的面積;x=x(x,y)是全局坐標系下的位置矢量;Hi(x)是流體1的指示函數,若x位于流體1中其值為1,否則為0。

(3)

式中:β是控制界面跳躍厚度的斜率因子(在本文中取為6);Pi(ξ)是一個二次多項式,其系數要根據分界面的法向和曲率進行計算,使得THINC/QQ隱式地具有幾何型VOF算法的特點。

(4)

di根據式(4)構造非線性方程求解。對方程(1)應用有限體積離散。在Ωi內可以得到如下形式的半離散控制方程:

(5)

THINC/QQ算法默認對方程(5)采用三階TVD龍格-庫塔(RK3)法進行瞬態求解;在計算di時對常規四邊形和三角形單元分別采用9點和6點高斯面積分;在計算φij時,在Eij上采用4點高斯線積分。

1.2 二維THINC/QQ-SF算法基本原理

圖1 將單元節點分為懸掛節點和重構節點的流程Fig.1 The procedure of separating the cell nodes into hanging nodes and reconstructing nodes

THINC/QQ-FVMS3在計算di時,需要將目標多邊形單元分解為J個三角形子單元,共需要3J個高斯積分點[6]。因此,THINC/QQ-SF在每一時間步內要更為簡便;尤其是當網格拓撲關系不發生變化時(如固定網格或彈簧網格),單元虛擬簡化步驟只需在網格前處理環節進行一次。對于一般的數據結構,線性四邊形與含有一個懸掛點的三角形具有相同的存儲特征[7];但THINC/QQ-FVMS3自身無法區分這2種單元,而THINC/QQ-SF則能區別對待。

1.3 HRIC及多相流求解算法

不同于THINC類型算法,HRIC作為代數型VOF算法,是基于歸一化變量圖法 (normalized variable diagram,NVD)保證對流格式的有界性[9],同時實現自由面的高分辨率求解。本文根據STAR-CCM+的默認參數來設置HRIC的相關公式[2]。并且與THINC/QQ一致,本文中HRIC也采用RK3方法對方程(5)進行瞬態求解。

值得說明的是,HRIC的公式中有低階對流格式(一階迎風)的貢獻。當流動方向與網格的單元面不正交時,一階迎風格式會使流場變量發生假擴散[10]。所以,在流體分界面發生大變形時,HRIC可能將界面由“清晰”變為“模糊”。

此外,本文還采用壓力隱式算子分裂算法 (pressure implicit with splitting of operators,PISO)[11]求解不可壓縮真實流動,具體按照文獻[12]的方式將其應用于非結構同位網格兩相流。對于動量方程,采用一階隱式歐拉法進行瞬態求解。

在時間步長控制方面,本文的庫朗數Co按照OpenFOAM的方式進行定義,對于Ωi有[7]:

(6)

并且同時采用2個最大庫朗數,其中Coφmax適用于主體界面單元(0.01<φ<0.99),Comax適用于其他單元。

2 典型對流算例的模擬與分析

VOF對流算例采用定量形狀誤差[4]:

(7)

式中:Ncel是網格單元總數;φei代表 VOF值的初始解(或精確解);φni代表數值解。

2.1 開槽圓盤剛體旋轉

首先以廣泛使用的Zalesak開槽圓盤剛體旋轉算例[13]進行測試。該圓盤的圓心位于(0.5,0.75),半徑為0.15,開槽區域為:

{(x,y)|(|x-0.5|≤0.025)且

(0.6≤y≤0.85)}

(8)

旋轉速度場為:

(u,v)=(ω(y-0.5),ω(0.5-x))

(9)

其中x、y∈[0,1]。在本文,ω=2π,且最大求解時間T=1, 相當于繞中心(0.5,0.5)旋轉一周;而庫朗數設置為:Coφmax=0.25且Comax=1.0×1010。

本例所用的測試網格包括A和B這2類二維切割體網格,每類網格又分為粗、細2個級別,分別記作A1、A2和B1、B2。每個網格都有2個內部加密區域,分別定義為:

{(x,y)|(|x-0.5|≤0.17)且

(|y-0.25|≤0.17)}

(10)

{(x,y)|(|x-0.5|≤0.17)且

(|y-0.75|≤0.17)}

(11)

圖2對網格A1和B1及其初始流體分界面進行了展示。在粗、細網格的每個邊上分別均勻分布50和100個單元。

圖2 開槽圓盤剛體旋轉的兩類粗網格及其初始流體分界面Fig.2 Two types of coarse meshes for the solid-body rotation of a slotted disk with the initial interface

表1列出了每個網格的單元總數,以及在t=T時刻的E(L1)數值誤差。在當前的計算條件下,HRIC的誤差是THINC/QQ-SF的3倍以上。

表1 開槽圓盤剛體旋轉算例在t=T時刻的E(L1)數值誤差Table 1 The numerical errors of E(L1) for the solid-body rotation of a slotted disk at t=T

在定性方面,圖3展示了在網格A1和B1上的代表性VOF等值線(0.05、0.5和0.95)。其中,THINC/QQ-SF算法能夠在二維切割體網格上保持流體分界面的緊致性;而HRIC算法的數值耗散性則比較明顯,尤其是對于較粗的網格。

注:在本文所有彩圖中值為0.05、0.5和0.95的VOF等值線分別用藍色、綠色和紅色線表示圖3 開槽圓盤剛體旋轉算例在t=T時刻網格A1和B1上0.05、0.5和0.95的VOF等值線Fig.3 The VOF contour lines of 0.05, 0.5 and 0.95 on Mesh A1 and B1 for the solid-body rotation of a slotted disk at t=T

2.2 圓盤剪切變形

在二維剛體旋轉算例的基礎上,采用經典的single vortex圓盤剪切變形算例[14]做進一步測試。同樣地,圓盤的圓心位于(0.5,0.75),半徑為0.15;其速度場通過流函數進行定義:

(12)

其中:x、y∈[0,1];庫朗數Coφmax為0.3且Comax為1.0×1010;T=8。

本例的網格類型也與開槽圓盤剛體旋轉一致,只是2個內部加密區域的定義稍有差異,分別為:

{(x,y)|(|x-0.5|≤0.187 5)且(|y-0.25|≤

0.187 5)}、{(x,y)|(|x-0.5|≤0.187 5)且

(|y-0.75|≤0.187 5)}

(13)

在粗、細網格的每個邊上分別均勻分布64和128個單元。表2列出了每個網格的單元總數,以及在t=T時刻的E(L1)數值誤差。在本文計算條件下,HRIC的誤差同樣都大于THINC/QQ-SF的誤差,且前者處于后者的1～2倍。

表2 圓盤剪切變形算例在t=T時刻的E(L1)數值誤差Table 2 The numerical errors of E(L1) for the single-vortex shearing flow of a disk at t=T

圖4定性地展示了在中間時刻t=T/2流體分界面(即VOF的0.5等值線)結果,可以看到,即使在分界面被嚴重拉伸變形的情況下,THINC/QQ-SF算法依然能在粗密網格之間實現分界面的光順過渡,而且比HRIC的數值解耗散性更低。

圖4 圓盤剪切變形算例在t=T/2時刻網格A2和B2上的流體分界面及其附近單元Fig.4 The interfaces along with the neighboring cells on mesh A2 and B2 for the single-vortex shearing flow of a disk at t=T/2

3 潰壩兩相流算例模擬與分析

為了進一步驗證二維THINC/QQ-SF算法的性能,采用典型二維潰壩算例[15]進行研究,包括有障礙物和無障礙物2種條件。計算域在水平x方向長0.584 m,在垂直y方向高0.34 m,初始水域為:

{(x,y)|(0≤x≤a)且(0≤y≤2a)},

a=0.146 m

(14)

在物理屬性方面,水和空氣的密度分別設置為1 000 kg/m3和1.0 kg/m3,而運動粘度分別為1.0×10-6m2/s和1.48×10-5m2/s;重力加速度值g=9.81 m/s2;且不考慮表面張力。庫朗數設置為:Coφmax=0.25且Comax=0.5。

3.1 有障礙物潰壩流動

網格的單元總數為9 841,并進行局部加密。在邊界條件方面,頂部為壓力入口(壓力等于0 Pa),其余均為不可滑移壁面。圖5為典型時刻的VOF等值線對比結果,其中,THINC/QQ-SF和HRIC都能捕捉到流動的主要特征、所得的自由面在粗密網格之間可以光順過渡;但THINC/QQ-SF的自由面更緊致,尤其是在界面發生大變形時。

圖5 有障礙物潰壩的瞬時0.05、0.5和0.95 VOF等值線Fig.5 The transient VOF contour lines of 0.05, 0.5 and 0.95 for the dambreak flow with an obstacle

3.2 無障礙物潰壩流動

計算域外尺寸和邊界條件與有障礙物潰壩一致,網格單元總數為7 424。圖6展示了在t=1.0 s的VOF等值線,雖然HRIC同樣能捕捉到典型位置的氣泡,但是THINC/QQ-SF在控制界面厚度方面更有優勢。

圖6 無障礙物潰壩在1.0 s的0.05、0.5和0.95VOF等值線Fig.6 The VOF contour lines of 0.05, 0.5 and 0.95 for the dambreak flow without any obstacles at 1.0 s

圖7 無障礙物潰壩流的自由面無量綱位移時歷曲線Fig.7 Time histories of the non-dimensional interface displacements for the dambreak flow without any obstacles

在垂向位移方面,本文所用的2種VOF算法數值解與文獻[12]的CICSAM數值解(采用8 400個網格單元)幾乎重合,并都與文獻[16]的實驗值非常一致。

在水平位移方面,不同VOF算法的數值解與實驗值[16-17]相比都存在一定的偏差,但整體趨勢比較吻合。這與水密閘門的運動、實驗測量方法等都有一定關系,有系統誤差的影響。

整體來講,在界面形變相對較小時,HRIC的界面質量與THINC/QQ-SF無顯著差別。當界面出現較大翻卷(尤其是破碎)時,THINC/QQ-SF更有助于抑制VOF的數值耗散,并能很好地適用于二維切割體網格。

4 結論

1)本文將THINC/QQ-SF算法的原理應用到二維切割體網格中,使得原始THIINC/QQ算法能夠用于懸掛點式四邊形/三角形單元的界面重構。而且,與THINC/QQ-FVMS3相比,THINC/QQ-SF需要的面高斯積分點更少、更簡便。

2)針對典型VOF對流算例(開槽圓盤剛體旋轉和圓盤剪切變形)以及潰壩兩相流算例,本文設計了一系列二維切割體網格,驗證了所用擴展方法的可行性。同時,與代數型的HRIC算法相比,THIINC/QQ-SF的精度更高、得到的自由面更緊致,能夠很好地繼承原始THIINC/QQ的優點。

3)對于需要局部加密或計算域不規則的二維數值模擬,如波浪傳播、經過海底凸起后的波浪變形、在斜坡上的波浪破碎等自由表面問題,THIINC/QQ-SF算法通過與二維切割體網格進行配合,具有很好的應用前景。