基于層次-批發(fā)式maes方法的液體晃蕩數(shù)值模擬

基于層次-批發(fā)式maes方法的液體晃蕩數(shù)值模擬

1流場模擬研究堡壘方法是1988年提出的一種運動界面跟蹤技術。該方法具有隱式跟蹤運動界面、自動處理拓撲變形的優(yōu)點，克服了波前跟蹤方法難以處理復雜物質界面和拓撲變形的弱點。此外,該方法不涉及坐標變換,可使用的差分格式較多,容易設計出高精度的格式。它的缺點是守恒性不能很好的滿足,以及對大梯度界面的光滑化處理造成的界面失真。目前,該方法已在圖象處理、晶體生長、流場模擬以及計算機可視化等領域廣泛應用。在流場模擬領域,國內(nèi)外已有的研究主要有激波追蹤、潰壩模擬等。載液船舶液艙內(nèi)液體的劇烈晃蕩所引起的載荷效應對船舶穩(wěn)性及液艙結構存在很大的安全隱患,此類液體晃蕩現(xiàn)象引起了人們的深刻關注,這些載荷與效應已成為載液船舶安全性評估的重要內(nèi)容之一。因此,對于液艙內(nèi)液體晃蕩特性的研究具有現(xiàn)實的工程意義。目前,針對該領域的研究已經(jīng)發(fā)展了多種數(shù)值技術和方法,朱仁慶對其作了較全面的綜述。本文首次嘗試性地將Level-set方法這種高效的界面追蹤技術應用于晃蕩液體自由液面的追蹤,結合其他數(shù)值技術成功模擬了二維液體晃蕩現(xiàn)象。2方程和數(shù)值法的控制2.1動力學離散方程控制方程采用統(tǒng)一的形式予以描述,整個流場由氣體和液體兩種介質組成,存在一物質界面,即自由液面。密度和黏性系數(shù)在空間不均勻分布,因此連續(xù)方程和動量方程寫成如下守恒形式。連續(xù)方程:?ρ?t+??(ρV)=0(1)動量輸運方程:dVdt=Fb+fb-1ρ?p+μρ?2V(2)式中:V是流場速度,ρ是流體密度,μ為流體動力黏性系數(shù),Fb是流體質量力,fb為由于流體非勻速運動而引起的慣性力。采用有限差分法對流場控制方程進行求解,以壓力修正法作為求解內(nèi)核。在MAC交錯網(wǎng)格中采用LeapFrog差分格式和DuFort-Frankel差分格式(無條件穩(wěn)定)的組合對N-S方程進行離散。分別在i+1/2,j上建立速度分量u的動力學離散方程,在i,j+1/2上建立速度分量v的動力學離散方程,在i,j上建立離散連續(xù)方程。各物理量在網(wǎng)格單元上的定義如圖1所示,其中速度分量定義在各邊中點處,壓力p和Level-set函數(shù)φ定義在網(wǎng)格中心;兩個虛線框分別表示x向和y向的動量輸運體,質量輸運體為實線框。連續(xù)方程的離散采用中心差分格式。離散格式的動量輸運方程(x方向):un+1i+1/2,j-ui+1/2,jΔt=fx-pn+1i+1,j-pn+1i,jρΔx-(u?u?x)i+1/2,j-(v?u?y)i+1/2,j+ν[(?2u?x2)i+1/2,j+(?2u?y2)i+1/2,j](3)離散格式的連續(xù)方程:ρn+1i,j-ρni,jΔt+(ρu)n+1i+1/2,j-(ρu)n+1i-1/2,jΔx+(ρv)n+1i,j+1/2-(ρv)n+1i,j-1/2Δy=0(4)其中,廣義質量力fx=-¨X(t)-gsinθ+2v˙θ+xθ2+y¨θ?X(t)為平蕩位移,X(t)=X0cos(ωxt+εx),θ為橫搖角,θ(t)=θ0cos(ωθt+εθ)。(u?u?x)i+1/2,j=ui+1/2,jui+3/2,j-ui-1/2,j2Δx(5)(v?u?y)i+1/2,j=14(vi,j-1/2+vi,j+1/2+vi+1,j-1/2+vi+1,j+1/2)ui+1/2,j+1-ui+1/2,j-12Δy(6)(?2u?x2)i+1/2,j=ui+3/2,j-un+1i+1/2,j-un-1i+1/2,j+ui-1/2,jΔx2(7)(?2u?y2)i+1/2,j=ui+1/2,j+1-un+1i+1/2,j-un-1i+1/2,j+ui+1/2,j-1Δy2(8)2.2網(wǎng)格收斂條件s通過求解動量方程得到新時刻的速度值往往不能滿足連續(xù)條件,因此需要對速度進行修正。本文采用壓力修正法,結合超松弛迭代技術(SOR),對速度和壓力進行反復迭代調整,使之滿足連續(xù)方程。從動量方程可以得到壓力和速度之間的修正關系式:u(m)i+1/2,j=u(m-1)i+1/2,j+ΔtΔpi,jρi+1/2,j(xi+1-xi)(9)u(m)i-1/2,j=u(m-1)i-1/2,j-ΔtΔpi,jρi-1/2,j(xi-xi-1)(10)v(m)i,j+1/2=v(m-1)i,j+1/2+ΔtΔpi,jρi,j+1/2(yj+1-yj)(11)v(m)i,j-1/2=v(m-1)i,j-1/2-ΔtΔpi,jρi,j-1/2(yj-yj-1)(12)上標m、m-1分別表示第m、m-1次的迭代值。用S表示式(4)左邊的值,根據(jù)以上四式得:Sm=Sm-1+ΔtΔpi,j[1Δx2(ρn+1i+1/2,jρni+1/2,j+ρn+1i-1/2,jρni-1/2,j)+1Δy2(ρn+1i,j+1/2ρni,j+1/2+ρn+1i,j-1/2ρni,j-1/2)](13)假設Sm=0,則有:Δpi,j=-Sm-1/{Δt[1Δx2(ρn+1i+1/2,jρni+1/2,j+ρn+1i-1/2,jρni-1/2,j)+1Δy2(ρn+1i,j+1/2ρni,j+1/2+ρn+1i,j-1/2ρni,j-1/2)]}(14)pmi,j=pm-1i,j+Δpi,j(15)速度和壓力迭代步驟:首先計算壓力修正值Δpi,j,然后修正所有網(wǎng)格各邊上的速度值,完成一步迭代;重復上述迭代,直至所有網(wǎng)格均滿足收斂條件S≤ε,ε為一小量。為了加快收斂速度,采用超松弛迭代技術,超松弛因子ω取為1.9(在不同的取值下進行了計算比較,取1.9時迭代步數(shù)較少)。2.3模型的建立與求解Level-set方法的基本原理是把運動界面用Level-set函數(shù)φ(x,t)的零等值面來描述,φ(x,t)滿足相應的控制方程。通過求解計算域內(nèi)φ(x,t)的值確定其等值面,實現(xiàn)運動界面的追蹤。構造函數(shù)φ(x,t)滿足:φ(x,t)={d(x,Γ(t))x∈Ω1(t)0x∈Γ(t)-d(x,Γ(t))x∈Ω2(t)(16)Ω1,Ω2分別表示液體區(qū)域和氣體區(qū)域,Γ(t)表示t時刻的自由液面,d(x,Γ(t))表示流場中任意點到自由液面的距離。φ(x,t)滿足對應于流場速度場的輸運方程和重新初始化方程:dφdt=?φ?t+V??φ=0(17){φt=sign(φ0)(1-|?φ|)φ(x,0)=φ0(18)式中:sign(φ0)是符號函數(shù),φ0是上一時間步求得的值。為了方便求解,將sign(φ0)光滑化為:signε(φ0)=φ0φ02+ε2(19)采用積分平均型Superbee-TVD格式對輸運方程(17)進行求解。這種格式的基本思想是:將求解函數(shù)構造成分片的線性函數(shù),用特征線方法推進時間步,并利用Superbee限制器對斜率加以限制,以使格式保持單調,具有TVD(TotalVariationDiminishing)性質。針對本文中的二維問題,采用分步法進行計算,將方程(17)分裂成式(20)的形式。前半個時間步在x方向上計算得中間值φn+1/2,再以φn+1/2為初值在后半個時間步在y方向上計算得φn+1。前后兩個時間步的方程分別為:12?φ?t+u?φ?x=0,12?φ?t+v?φ?y=0(20)以x方向的輸運方程為例,在時間和空間控制體積元Ωin=[xi-1/2,xi+1/2]×[tn,tn+1]上積分得:∫xi-1/2xi+1/2[φn+1(x)-φn(x)]dx+∫tntn+1[uφi+1/2(t)-uφi-1/2(t)]dt=0(21)分別定義空間積分平均和時間積分平均:φˉin=1Δx∫xi-1/2xi+1/2φn(x)dx?uφˉi+1/2n?=1Δt∫tntn+1uφi+1/2(t)dt(22)通過沿特征線的輸運近似(Δt滿足CFL條件),并采用迎風思想,可將式(21)離散成如下格式:φˉin+1/2=φˉin-2rmax(0,u)(φˉin-φˉi-1n)+rΔxmax(0,u)(c-1)(Sin-Si-1n)-2rmin(0,u)(φˉi+1n-φˉin)+rΔxmin(0,u)(c+1)(Si+1n-Sin)(23)其中:c=uΔtΔx,Sin=?φˉin?x,r=ΔtΔx?Sin是空間控制元上的斜率,為了在界面處避免振蕩的發(fā)生,保證格式的單調,Sin應加以限制,取為:Sin=?(θin)(φˉi+1n-φˉin)/Δx(24)其中:?(θin)是Superbee限制器,取為:Sin=?(θin)=max(0,min(1.0,2θin),min(2.0,θin))(25)θin={φˉin-φˉi-1nφˉi+1n-φˉin,|φˉi+1n-φˉin|＞ε0,|φˉi+1n-φˉin|≤ε(26)ε為一小量,本文取10-7(視網(wǎng)格粗細而定),空間積分平均φˉin可近似取為φin。重新初始化方程(18)的求解采用修正的Godunov方法。將方程(18)寫成:φt+(S(φ0)φxφx2+φy2)φx+(S(φ0)φyφx2+φy2)φy=S(φ0)(27)其中,S(φ0)是符號函數(shù)sign(φ0)。定義s=S(φ0)(|φx+|-|φx-|)φx+-φx-(28)對φx做如下離散:φx={φx-當S(φ0)φx+≥0,S(φ0)φx-≥0φx+當S(φ0)φx+≤0,S(φ0)φx-≤00當S(φ0)φx+＞0,S(φ0)φx-＜0φx-當S(φ0)φx+＜0,S(φ0)φx-＞0,s＞0φx+當S(φ0)φx+＜0,S(φ0)φx-＞0,s≤0(29)其中,φx+,φx-為由5階WENO(weightedessentiallynon-oscillatory)格式構造出來的左側和右側的導數(shù)離散值。時間導數(shù)的離散采用三階TVD-Runge-Kutta方法。通過上述空間導數(shù)的離散,方程(18)等價于常微分方程:dφdt=L(φ)(30)其中:L(φ)是空間導數(shù)項的離散算子。三階TVD-Runge-Kutta法分為三步:{φ(1)=φ(0)+ΔtL(φ(0))φ(2)=34φ(0)+14φ(1)+14ΔtL(φ(1))φ(3)=13φ(0)+23φ(2)+23ΔtL(φ(2))(31)對方程(18)進行時間步的迭代直至解穩(wěn)定,即φt?0,在數(shù)值上要求:∑i,j|φi,jn+1-φi,jn|＜ε(32)其中,ε為一小量,本文取為0.005Δx。3計算值的示例3.1自由土地內(nèi)充放電溫度分布此算例模擬了二維矩形液箱內(nèi)具有初始液面形狀的駐波的自由運動,駐波的初始波形表示為:η=acos(2πx/L),其中a為波幅,x為長度方向的坐標,L為液箱長度。在理想流體和微幅波的假設前提下,液體在矩形液箱內(nèi)晃蕩的波面固有諧振頻率為:ωn2=2gmnπ·th(2mnπh/L)/L,其中:L是液箱長度,h是液深,g是重力加速度,n表示固有頻率的階數(shù),m表示L長度范圍內(nèi)液面所包含的波數(shù)。因此,線性勢流理論下液面波動的解析表達式可表示為η=acos(2πx/L)cosω1t,ω1為一階線性波面固有頻率。為了和線性勢流解析解比較,取黏性系數(shù)γ=0。分別在兩種不同的波幅a=0.05h和a=0.13h下進行計算,無因次化處理后的計算結果以液箱長度中點處x=L/2的自由液面起伏的時間歷經(jīng)曲線給出,如圖2-b、圖3-b所示,圖中實線為本文的數(shù)值結果,虛線為線性勢流理論下的解析解。此外,還將計算的晃蕩頻率與其余三種數(shù)值方法得到的結果進行了比較:VOF法(朱仁慶,2002)、有限元法(Greavesetal,1997)和有限水深駐波的傅里葉數(shù)值解法(TsaiandJeng,1994)。從圖2-b和圖3-b看出,本文采用的數(shù)值方法算得的液面起伏時間歷經(jīng)與解析解基本吻合,周期和波幅的微小偏差的產(chǎn)生是由于非線性的影響,類似的現(xiàn)象在文獻中也曾提及。計算得到的晃蕩頻率與其它數(shù)值方法得到的結果也保持了一致?？梢?Level-set法用于自由液面的追蹤是可行的。3.2壓電時間統(tǒng)計此算例模擬了Arai于1984年做的一個模型試驗。液箱裝有水,做受迫正弦橫搖運動,幅度θ=4°,頻率f=0.89Hz,與液體晃蕩的一階自振頻率0.8917Hz相近,橫搖中心位于箱底中部。液箱尺寸L=0.9m,H=0.6m,載液深度h=0.45m。根據(jù)所提供的實驗數(shù)據(jù),本文同樣給出了左壁面H/2和3H/4處的壓強時間歷經(jīng)曲線