大密度比兩流體系統(tǒng)計(jì)算的一種改進(jìn)模式

汪迎春,唐洪武,汪德,唐 亮

(1.蚌埠學(xué)院食品與生物工程系,安徽蚌埠 233030;2.河海大學(xué)水利水電學(xué)院,江蘇南京 210098;3.中國(guó)水電建設(shè)集團(tuán)國(guó)際工程有限公司,北京 100048)

對(duì)于密度相近的兩流體,Δt時(shí)間內(nèi)網(wǎng)格單元密度可以近似認(rèn)為不變,但如果兩流體密度差異很大或者時(shí)間步長(zhǎng)Δt較大,網(wǎng)格單元密度在Δt時(shí)間內(nèi)相對(duì)變化需要考慮,這時(shí)密度的變化將對(duì)計(jì)算結(jié)果產(chǎn)生影響.大密度比兩流體問(wèn)題的求解需增加密度方程

式中:ρ——流體密度;uj——流速標(biāo)量;t——時(shí)間;xj——空間坐標(biāo).

目前常見(jiàn)的解決此類問(wèn)題的方法有2種:一種是令時(shí)間步長(zhǎng)Δt網(wǎng)格單元密度不變的方法,這種方法適用于弱密度比的兩流體問(wèn)題;另一種是通過(guò)增加密度方程并進(jìn)行耦合求解的方法,這種方法適用于大密度比的兩流體問(wèn)題.直接求解密度方程,會(huì)使計(jì)算過(guò)程中的非線性問(wèn)題更加突出,還會(huì)產(chǎn)生密度越界現(xiàn)象.本文通過(guò)對(duì)動(dòng)量方程的一些變換,把密度方程和動(dòng)量方程統(tǒng)一為一個(gè)積分形式,從而有效地解決了大密度比問(wèn)題.

1 控制方程

對(duì)于滿足連續(xù)條件不可壓縮的兩相流,其連續(xù)方程、密度方程和動(dòng)量方程為

式中:fi——質(zhì)量力;p——壓強(qiáng);μ——?jiǎng)恿︷ざ?

式(3)左邊第1項(xiàng)可變換為

代入式(3)得

方程(4)的離散采用非結(jié)構(gòu)同位網(wǎng)格,物理變量布置在單元中心,如圖1所示,用有限體積法離散.將式(4)寫(xiě)成微分方程的積分形式

式中:φ——物理變量;φp——控制體中心的物理變量;Sφ——源項(xiàng);n——單元邊界A的外法向單位向量;d A——面積分微元;d V——體積分微元;u——流速向量.

從方程(5)可直觀地看出,對(duì)流項(xiàng)的穩(wěn)定性得到了改善.經(jīng)以上變換,方程(2)和(3)組合成了方程(5),只需求解方程(1)和(5)即可.

對(duì)于方程(5),對(duì)流項(xiàng)離散可用二階迎風(fēng)格式,單元邊界均值 φe由迎風(fēng)單元直接線性插值得到,

式中:φu——迎風(fēng)單元中心物理量;(gradφ)u——迎風(fēng)單元物理量梯度;d r——從迎風(fēng)單元中心至單元面的距離矢量.

其中的單元中心梯度由高斯公式計(jì)算:

圖1 網(wǎng)格單元Fig.1 Distribution of cell variables

式中:Aji——單元j面的i方向面積分量;φj——界面處物理變量均值,可直接通過(guò)線性插值求得,如果是高度扭曲的畸變網(wǎng)格,則會(huì)引入網(wǎng)格誤差;α——保證重構(gòu)不導(dǎo)致局部極值產(chǎn)生的限制性因子,0≤α≤1.擴(kuò)散項(xiàng)的離散方法可參考文獻(xiàn)[1].

2 速度壓力耦合與相界面的處理

不可壓縮流體N-S速度壓力耦合方程的求解方法有很多,由于人工壓縮法[2-6]形式簡(jiǎn)單,且可以獲得很高的精度,本文采用人工壓縮法.其中,偽時(shí)間導(dǎo)數(shù)項(xiàng)的離散采用一階向后差分格式,物理時(shí)間項(xiàng)的離散采用二階向后差分格式,對(duì)流項(xiàng)的離散采用二階迎風(fēng)格式,擴(kuò)散項(xiàng)的離散采用中心差分格式.兩相流界面的追蹤采用Level Set方法[6-10].

3 算 例

本文對(duì)Rayleigh-Taylor不穩(wěn)定分層問(wèn)題[11]進(jìn)行了數(shù)值模擬.

Rayleigh-Taylor不穩(wěn)定分層系指上下兩層不可壓縮的黏性流體,如圖2所示.上層為重質(zhì)流體,其密度和運(yùn)動(dòng)黏度分別為 ρ1,ν1;下層為輕質(zhì)流體,其密度和運(yùn)動(dòng)黏度分別為 ρ2,ν2.ρ1/ρ2=998/499=2,ν1/ν2=1.初始界面擾動(dòng)為

圖2 Rayleigh-Taylor不穩(wěn)定分層問(wèn)題定義Fig.2 Definition of Rayleigh-Taylor unstable layer problem

式中:u,v——x,y方向的流量;Uref——參考點(diǎn)流速.

本研究取α=0.25,Re=283,左右兩邊取對(duì)稱邊界,上下兩邊取固壁邊界.

量綱一的初始?jí)毫θ?/p>

其中

式中:pref——參考點(diǎn)壓力;g——重力加速度;p(y)——基準(zhǔn)壓力.因

此,上層流體

下層流體

其中

圖3為采用大密度比常用處理方法和本文處理方法計(jì)算得到的分層界面形狀.由圖3可以看出,2種大密度比處理方法計(jì)算得到的分層界面形狀很接近,由于離散格式未采用高精度的格式,故不穩(wěn)定界面向下翻轉(zhuǎn)的卷舌處模擬還不夠精細(xì).圖3(a)由于采用常用的處理方法,界面坦化嚴(yán)重,結(jié)果和理論分析成果[12]出現(xiàn)較大的差別.

圖3 計(jì)算的分層界面形狀Fig.3 Shapes of calculated layer interfaces

4 結(jié) 語(yǔ)

本文把密度方程和動(dòng)量方程組合成一個(gè)方程,變換為一個(gè)積分形式,提出了一種新的處理大密度比兩流體系統(tǒng)的方法.該方法有效地避免了兩方程的耦合問(wèn)題以及單獨(dú)求解密度方程帶來(lái)的一些不合理現(xiàn)象.數(shù)值試驗(yàn)表明,本文方法是有效的,可用于大密度比的兩流體問(wèn)題的數(shù)值模擬和分析.

[2]薛具奎,趙金保.人工壓縮法數(shù)值模擬非定常不可壓縮熱對(duì)流流場(chǎng)[J].西北師范大學(xué)學(xué)報(bào):自然科學(xué)版,1998,34(2):32-37.(XUE Ju-kui,ZHAO Jin-bao.Numerical simulation for unsteady incompressible thermoconvective flow with artificial compressibility method[J].Journal of Northwest Normal University:Natural Science Edition,1998,34(2):32-37.(in Chinese))

[3]陳傳淡.守恒雙曲型方程式激波計(jì)算區(qū)域分解法中的追蹤及人工壓縮方法[J].廈門(mén)大學(xué)學(xué)報(bào):自然科學(xué)版,1999,38(3):329-333.(CHEN Chuan-dan.Tracing and artificial compression method of theshock′s calculation for the conservative hyperbolic equation in domain decomposition[J].Journal of Xiamen University:Natural Science,1999,38(3):329-333.(in Chinese))

[4]朱自強(qiáng),賈劍波.三角翼大迎角不可壓粘流的數(shù)值模擬[J].力學(xué)學(xué)報(bào),1996,28(6):736-740.(ZHU Zi-qiang,JIA Jian-bo.Numerical Simulation of incompressible viscous flow over a deltawing at high angle of attack[J].Acta Mechanica Sinica,1996,28(6):736-740.(in Chinese))

[6]汪迎春.運(yùn)動(dòng)界面模擬技術(shù)及在環(huán)境分層流中的應(yīng)用[D].南京:河海大學(xué),2007.

[7]OSHERS,SETHAIN J A.Fronts progagating with curvature-dependent speed:algorithms based on Hamilton-Jacobi formulations[J].J Comput Phys,1988,79(1):12-49.

[8]劉儒勛,王志峰.數(shù)值模擬方法和運(yùn)動(dòng)界面追蹤[M].合肥:中國(guó)科技大學(xué)出版社,2001.

[9]TSAI Y H R.Rapid and accurate computation of the distance function using grids[J].JComput Phys,2002,178(1):175-195.

[11]FELECY F J,PLETCHER RH.Numerical simulation of free surface flows in closed containers using free surface capture approach,in advances in computational methods in fluid dynamics[J].Am Soc Mech Eng,1994,196:33-51.

[12]ZHAO Yong,TAN Hsiang-hui,ZHANG Bai-li.A high-resolution characteristics-based implicit dual time-stepping VOF method for free surface flow simulation on unstructured grids[J].Journal of Computational Physics,2002,183(1):233-273.