劉 君,白曉征,王 巍

(1.國防科技大學航天與材料工程學院,湖南長沙 410073；2.空軍裝備研究院裝備總體論證研究所,北京 100076)

0 引 言

網(wǎng)格分區(qū)技術是結構網(wǎng)格求解復雜外形的基礎,廣泛應用的主要有重疊(overlapping)網(wǎng)格技術和對接(patched)網(wǎng)格技術兩類,計算過程中,首先采用流體力學解算器分別求解各子區(qū)域網(wǎng)格上物理量,然后通過邊界上信息傳遞,從而完成整個流場的計算。實際問題中相鄰子區(qū)域的邊界網(wǎng)格或重疊部分網(wǎng)格很少在位置上正好一一對應,因此需要采用一定的算法傳遞相鄰子區(qū)域的信息。對于包含有激波等間斷的流場模擬,目前常用線性插值或其他高精度插值方法,插值表達式通常為多項式,用網(wǎng)格自身及其周圍單元的信息來確定多項式的系數(shù)；線性插值優(yōu)點是滿足單調性、守恒性、形式簡單、效率高,但是只有一階精度,數(shù)值耗散大；高精度插值提高精度,但是涉及到網(wǎng)格內物理量導數(shù)計算,而且在應用非結構動網(wǎng)格技術和插值方法求解分離問題時發(fā)現(xiàn),網(wǎng)格重構以后流場出現(xiàn)的非物理振蕩會導致氣動力明顯的波動[1,2]。

在文獻[1]中給出一種新的信息傳遞算法,充分利用時間推進的特點,采用移動網(wǎng)格的方法,隨著時間推進完成舊網(wǎng)格信息到新網(wǎng)格的傳遞。這種新的流場信息傳遞算法與插值方法的原理完全不同,理論可以證明信息傳遞過程沒有引入誤差。文獻[1]中針對一維和二維問題,比較了新的信息傳遞算法和插值方法,給出了傳遞過程中時間推進步長的穩(wěn)定性條件。

本文進一步分析了插值方法引入誤差的原因,然后在介紹新的信息傳遞算法原理基礎上,推導了三維情況下傳遞過程時間推進步長的穩(wěn)定性條件,最后給出了三維算例。

1 插值方法引入誤差的原因分析

由于系數(shù)γl與網(wǎng)格空間分布間距Δx相關,計算得到的流場參數(shù)與網(wǎng)格分布相關,在其他條件相同的情況下,不同網(wǎng)格分布在同一物理位置的數(shù)值解也不相同。

設有兩套網(wǎng)格,在同一物理位置x P處,根據(jù)差分方程得到第二套網(wǎng)格數(shù)值解為uP,根據(jù)第一套網(wǎng)格結點變量插值得到的該點值為u T。插值方法是相同時刻不同空間位置物理量之間的關聯(lián),可以寫為二者之間距離的函數(shù):

數(shù)值解uP滿足修正方程,將上式帶入修正方程(1)中,略去高階小量,得到:

理論上uP和u T時間空間相等,如果假設uP和u T隨時間變化導數(shù)相等,根據(jù)上式可以推導出插值過程不影響差分方程數(shù)值解的條件:

由于插值方法和有限差分法相互獨立,不可能嚴格滿足以上條件,因此得出結論:插值過程必然引入誤差,導致流場計算過程精度降低。

根據(jù)CFD理論,對數(shù)值解影響較大的是修正方程第一項,假設第一項中沒有插值誤差就認為對有限差分法不影響。如果差分格式具有一階精度,為了保證修正方程第一項中系數(shù)γ2～o(Δx)不受影響,要求式(4)為高階小量,即r≥4；如果差分格式具有二階精度,保證第一項中系數(shù) γ3～o(Δx2),要求 r≥6,線性插值和ENO插值都達不到這樣的要求,因此插值以后出現(xiàn)明顯波動[1]是不可避免的。

2 動網(wǎng)格傳遞流場信息的基本原理

在兩套網(wǎng)格間進行流場插值或信息傳遞前,首先需要確定網(wǎng)格間的關系。以重構后新舊網(wǎng)格為例,需要確定新網(wǎng)格的網(wǎng)格單元中心位于舊網(wǎng)格區(qū)域的哪個單元內。在插值方法中,首先假設物理量的分布形式,如線性分布、拋物型分布,然后通過貢獻單元及其周圍單元的信息構造出分布系數(shù),最后得到新網(wǎng)格中心點的物理量,或者通過守恒型算法計算新網(wǎng)格單元內物理量的平均值。

與插值方法不同,信息傳遞算法并不顯式地指定物理量的分布形式。它的基本原理是,采用動網(wǎng)格技術,讓舊網(wǎng)格中心點在有限的時間內運動至新網(wǎng)格中心處,在移動網(wǎng)格的同時求解ALE形式的流動方程,從而得到新網(wǎng)格中心點的物理量。下面以二維為例介紹這種信息傳遞的原理。圖1(a)中,O點為n時刻舊網(wǎng)格的中心,n時刻物理量已知；P和Q是新網(wǎng)格中兩個單元的中心點。為了得到新網(wǎng)格P點的物理量,采用移動網(wǎng)格技術,保持O點附近網(wǎng)格不變形,在Δt時間內將舊網(wǎng)格單元從O點平移至P點,如圖1(b)所示,移動網(wǎng)格的同時求解流場,得到N+1時刻O′處物理量值,也即新網(wǎng)格單元中P點的值；同樣,在舊網(wǎng)格基礎上,還是從n時刻流場參數(shù)出發(fā),保證時間方向同樣推進Δt,同時平移網(wǎng)格使得O點在n+1時刻移到Q點,得到Q點n+1時刻的物理量。對新網(wǎng)格中每一單元重復此過程,即可得到新網(wǎng)格在n+1時刻的流場,在此基礎上放棄舊網(wǎng)格、采用新網(wǎng)格繼續(xù)計算,如圖1(c)所示。

圖1 信息傳遞原理演示Fig.1 Sketch map of the new method

從原理上可以看出,以上信息傳遞方法是一種全新的思想,與插值方法有明顯不同:

(1)插值方法在同一時間層內利用空間關聯(lián)完成信息傳遞,即用n時刻O點及其相鄰點參數(shù)得到n時刻P和Q點處流場參數(shù),然后采用新網(wǎng)格繼續(xù)計算；本文信息傳遞方法利用時間推進在兩個時間層內完成。

(2)上一節(jié)已經證明插值必然引入新的誤差,本文方法通過求解流動方程完成流場信息傳遞,得出的時刻流場參數(shù)與原有限體積方程完全相容,信息傳遞過程本身不會帶來額外誤差。

(3)插值方法與流動控制方程無關,其精度與假設的物理量的分布形式有關；本文方法通過求解流動方程實現(xiàn)信息傳遞,并不顯式假設物理量的分布形式,其精度即為流動求解器的精度。

采用本文方法進行信息傳遞時,為了保證n+1時刻的流場同步性,整個需要傳遞信息的區(qū)域應當使用根據(jù)穩(wěn)定性條件確定的統(tǒng)一的最小時間步長Δtmin。當新舊網(wǎng)格之間移動距離較大時,網(wǎng)格移動在一個Δt時間內無法實現(xiàn),需要采用多步進行。從CFD計算方法構造過程知道,采用時間顯式格式推進時,從n時刻推進到n+1時刻僅需要相鄰網(wǎng)格的n時刻流場,推進到n+2時刻則要相鄰網(wǎng)格n+1時刻流場,本質上涉及到相鄰網(wǎng)格以外更遠區(qū)域的n時刻流場。采用多步推進時,把移動區(qū)域擴大了,會使得計算過程更為復雜,計算量增大。因此,需要確定滿足Δtmin要求的最少時間推進步數(shù),這就涉及到移動網(wǎng)格過程中的穩(wěn)定性條件。

在文獻[1]中給出了一維和二維問題的穩(wěn)定性條件,下面推導三維問題的穩(wěn)定性條件。

3 三維穩(wěn)定性條件

如圖2示意的四面體,新網(wǎng)格中心點P位于中心點O的舊網(wǎng)格四面體ABCD中,假設O點在時間Δt內移動到P點,根據(jù)下式計算網(wǎng)格移動速度:

假設采用時間顯式格式的有限體積方法,穩(wěn)定性條件為:

進一步化為更嚴格的限制條件:

其中SX、SY和SZ為四面體ABCD在三個坐標平面的投影面積。聯(lián)立方程式(5)、(6),得到:

定義三維動網(wǎng)格Courant數(shù)NV為:

可以看出,NV與各點坐標絕對值大小及流動變量無關,僅取決于新舊網(wǎng)格的相對位置和舊網(wǎng)格的形狀參數(shù)。為便于分析,采用坐標系平移使O點位于原點,O點與舊網(wǎng)格的4個頂點相連,形成4個四面體,包含新網(wǎng)格P點的四面體記為ABCO,坐標系旋轉使得P點位于z軸負向,如圖2所示,設OP延長線與ΔABC平面交點為H。在 xoy平面內,ΔABC和D點的投影分別記為Δabc和d,由于P點位于四面體ABCO中,O點必然落在Δabc內,利用這一特性分析特征投影面積SZ的取值范圍。

圖2 NV取值范圍的證明過程Fig.2 Proving process of the varying range of NV

如果d點落在Δabc內部,那么 SZ=Sabc。如果d點不在Δabc內,分為d點在線段ab外側和d點在c點外側兩種情況,如圖2中示意,ab中點記為e,cd中點記為f。O點為四面體重心,有:xa+xb+xc+xd=0和y a+yb+y c+y d=0,推出 e點和 f點滿足的關系:

(1)對于d點在線段ab外側情況,兩個四邊形面積相等Safbc=Safbd=SZ,由于Sabc≥Safbc,得到:Sabc≥Z。

由于 NV≤3,對穩(wěn)定CFL數(shù)為1的顯式格式,根據(jù)式(7),網(wǎng)格移動可以在3步內完成。

按照以下關系可以得到滿足全流場穩(wěn)定性和動網(wǎng)格Courant數(shù)的時間推進步長:

即新網(wǎng)格中點P位于任意四面體ABCD中,三步內可以將O點移動到P實現(xiàn)流場信息傳遞。

4 算例驗證

設計如下圖3所示意的激波誘導物體運動問題,兩圓柱位于長∶寬∶高=2∶0.25∶1的矩形計算區(qū)域內,其中質心為(-0.4,0.78,0.125)的圓柱固定不動,質心為(-0.4,0.22,0.125)的圓柱的無量綱質量為m=1,不考慮重力和摩擦力等因素,圓柱運動特性完全由氣動力控制,可以自由移動。在環(huán)境大氣條件下,激波馬赫數(shù) M∞=2.96進入計算區(qū)域與物體相互作用,產生繞射、反射等非定常流動現(xiàn)象的同時,采用動網(wǎng)格技術開展流動模擬,具體計算方法、幾何守恒律等問題參見文獻[3]。

圖3中給出無量綱時間為t=0.2和t=0.4的壓力云圖,計算中在t=0.332時網(wǎng)格變形嚴重,需要進行局部網(wǎng)格重構,重構前后物體附近網(wǎng)格見圖4。分別采用線性插值、Weighted Least-Squares[4-5]插值方法和動網(wǎng)格信息傳遞方法得到重構后網(wǎng)格點的流場參數(shù),然后繼續(xù)計算。

圖3 計算區(qū)域、網(wǎng)格和不同時刻流場壓力云圖Fig.3 Computational model,grids and pressure contours

圖4 重構前后局部網(wǎng)格Fig.4 Grids beforeand after remeshing

圖5給出不同方法得到的圓柱受力情況,不同的信息傳值方法有明顯差異,動網(wǎng)格信息傳遞沒有引入新的誤差,作為比較基礎,分析如下:

(1)對于靜止圓柱,盡管物體附近網(wǎng)格沒有進行重構,由于處于運動圓柱插值區(qū)域馬赫錐內,還是受到插值誤差影響,大致t=0.36開始軸向力FX出現(xiàn)變化,插值誤差隨流動向下游傳播,到t=0.39前FX逐漸趨近于本文方法的計算曲線。法向力FY到了t=0.42以后三條曲線才開始明顯區(qū)別,說明插值誤差影響機理較為復雜。

(2)對于運動圓柱,線性插值以后立刻引起軸向力F X劇烈變化,通過流場動畫比對,在t=0.39以后計算曲線的兩次波動是插值誤差經過壁面和另一圓柱反射所致。法向力FY在t=0.35前兩種插值計算曲線和本文方法的比較一致,到了t=0.41后兩種插值計算曲線之間相差較大,本文方法計算曲線落在它們之間,變化規(guī)律不一致。

通過算例,可以看出插值方法引入誤差在流場內擴散規(guī)律較為復雜,本文方法計算曲線光滑,驗證了前面的分析結論。

圖5 上:靜止圓柱受力 下:運動圓柱受力Fig.5 Aerodynamic forces for the steady column(above)and the moving column(below)

5 結束語

理論和實際算例表明,新的信息傳遞算法比目前常用的插值方法精度高,通過本文工作把這一方法推廣到三維有限體積方法中,解決了其應用到復雜流場時遇到的難題。

[1]LIU J,BAI X Z,GUO Z,et al.A new method for transferring flow information among meshes[J].Computational Fluid Dynamics Journal,2007,15(4):509-514.

[2]白曉征,郭正,劉君,等.非結構動網(wǎng)格方法中的若干問題研究[A].中國第一屆近代空氣動力學與氣動熱力學會議論文集[C].綿陽,2006:301-306.

[3]劉君,白曉征,郭正,等.有相對運動多體動力學系統(tǒng)流動計算方法若干問題討論[J].空氣動力學學報,2008,26(1):14-18.

[4]ZEEUW D L D.A quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations[D].the University of Michigan,1993.

[5]HUNT JD.An adaptive 3d Cartesian approach for the parallel computation of inviscid flow about static and dynamic configurations[D].The University of Michigan,2004.