GEL方法在一維雙界面問題中的應用*

姚 陽

(中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理國防科技重點實驗室，四川綿陽621900)

將兩種不同坐標系下的計算格式耦合起來的Euler-Lagrange耦合計算方法，在許多計算性學科領域中應用前景非常廣闊。例如在爆炸力學問題的數值計算中，往往希望在整體計算區域中流體或大變形的部分，采用Euler方法計算，而對小變形以及固體材料的部分，采用Lagrange方法計算，以能夠對整體計算區域的不同區域部分獲得更細致的描述，如水下爆炸對船體的影響、低速侵徹問題和一些爆轟驅動問題等。對于這兩種不同計算格式，怎樣用合理有效的耦合方式、簡單的邏輯過程以及不對計算存貯方面有較高限制的Euler-Lagrange耦合方法，已越來越成為關注的熱點。

本文中，對一種較新的Euler-Lagrange耦合方法進行研究，它能夠自由地對標準Euler和Lagrange計算程序包進行連接，并且不必對他們進行任何實質上的改動，連接過程簡單易行［1］。這種方法由R.P.Fedkiw［2］提出，并用于多介質流問題的計算。對于Euler方法面對因運動界面而出現混合網格的難題，利用ghost fluid method(GFM)［3］的思想進行處理。E.Morano等［1］用這種方法進行了在流固耦合方面的研究，并稱為 ghost-fluid Euler-Lagrange(GEL)方法。相對 CEL［4］、PISCES［5］、PELE［6］等 Euler-Lagrange方法，GEL對自由邊界的處理過程和程序的編制簡單、計算量小。

在處理混合網格時，為得到ghost網格上的質量、動量和能量，需要將外插量從Euler真實域穿過界面延拓到ghost域。鑒于在有激波和具有固體不可穿透性質的界面相互作用的數值模擬中，如果選擇外插溫度或密度，在界面附近往往會出現所謂的overheating誤差。為削弱此誤差，用isobaric修正變量，即熵、T'(ρ)、e'(ρ)作為第3個外插量將減小這種誤差［7］。進一步用isobaric修正技術不僅將外插量向ghost區域延拓，而且還向緊鄰界面的Euler真實點延拓，可將誤差進一步減小。本文中，擴展了原自行編制的二維GEL程序［8］，將程序用于處理雙界面問題以及材料在激波作用下的變形。通過一維黎曼問題，計算當初始間斷形成的左右行進激波相遇后，反射波分別穿過界面傳播的情況，并通過對isobaric修正方法的具體應用，有效地減小overheating現象。

1 基本方程

1.1 Euler方程

理想可壓縮流體Euler控制方程是

式中:ρ為密度，u為速度矢量，p為壓力，e為質量內能，t為時間，I是單位矩陣。狀態方程形式取為

Euler計算采用二階SCB格式［9］編制的計算程序。

1.2 Lagrange方程

連續介質力學給出的基本方程有

式中:ρ、u和E分別表示質點密度、速度和體積內能，~b為體積體力，~σ為Cauchy應力張量。當采用流體彈塑性本構模型時，應力張量分解為

Lagrange計算采用DEFEL二維動力有限元程序［10］。

2 GEL方法

GEL方法中，計算區域中的Euler域和Lagrange域有各自生成的Euler固定網格和Lagrange隨體網格。兩個區域的交界面稱為Euler-Lagrange界面(E-L界面)。E-L界面將整個流場分為Euler和Lagrange域，而且被兩個區域共同享有。因為Lagrange節點以隨體的速度運動，因此，界面的幾何特性可以由Lagrange計算結果確定［2］。由于運動界面的關系，部分的Euler網格被Lagrange網格覆蓋，這些被覆蓋的Euler網格就是ghost網格，也稱為ghost域。因此界面將Euler域分為兩個部分:一部分是未被覆蓋的真實Euler域，另一部分是ghost域。

與GFM方法相同的是，ghost域上物理量定義點稱為ghost點，對物理量的定義可看作是對Euler域內邊界物理量的確定。

GEL耦合是根據E-L界面接觸間斷的性質，界面兩側應滿足的物理量條件，在兩個計算區域分別確定界面上連續或間斷物理量的過程中，內邊界條件的確定。對E-L界面上的切向速度、密度(或熵)，在界面兩邊區域上的這些間斷量互不相關，因此兩個計算區域分別由各自區域計算結果確定E-L界面上的值。而由界面兩側法向速度的連續性，ghost點的法向速度由Lagrange計算的E-L界面上的法向速度確定，即相當于對Euler域施加了速度邊界條件。利用壓力是連續量，以E-L界面附近Euler網格的壓力，在E-L界面可對Lagrange域施加壓力邊界條件。通過相互施加邊界條件，體現著Euler域和Lagrange域之間的相互作用。圖1是E-L界面連續量確定示意圖。耦合的過程需在每一時間步計算前進行，然后以統一的時間步長用標準的Euler和Lagrange計算程序分別進行獨立地計算。

圖1 E-L界面連續量的確定Fig.1 The definition of the continuous quantity in E-L interface

2.1 Ghost域的定義

首先通過程函數

從未被Lagrange域覆蓋的真實Euler域，外插壓力、熵或密度、以及速度到ghost域。式中:I表示外插量，n是Level set距離函數Φ梯度方向上的單位矢量。Level set函數是點到界面的距離函數。每一時間步，由Lagrange計算的E-L界面，都需重新建立Φ。采用文獻［2］中介紹的方法判別Φ的符號。在E-L界面上，Φ=0，n的方向從Euler域指向Lagrange域。

Ghost點的速度的切向分量來自外插速度的切向分量，法向分量還需由Lagrange計算的E-L界面上點的法向速度確定。以二維情況為例，此E-L界面上的點應是距ghost點最近的E-L界面線段上的點，速度值通過此點所在的界面線段兩節點的速度值線性內插得到［1-2］。最后，ghost點的速度vG可由一定外插形式的速度合成公式計算。采用反射外插公式［1-2］

式中:vJ表示距ghost點最近的E-L界面上點的速度，vext表示外插的速度，t為垂直于n的單位矢量。式(10)反映了界面的不可穿透性［1］。

2.2 Isobaric修正技術

除壓力和速度外，還需要得到ghost網格上第3個獨立的物理量。因為在有激波和具有固體不可穿透性質的界面相互作用的數值模擬中，界面附近往往會出現所謂的overheating誤差，而壓力和速度在界面保持一致。為削弱此誤差，除了用isobaric修正變量，即熵、T'(ρ)和e'(ρ)作為第3個外插量外，還可以用isobaric修正技術將誤差最小化［7］。在GEL方法中，外插時凡是滿足

所有Euler域格點的物理量也隨之改變。式中:Δx為網格寬，k為無物理意義的非零常系數，可根據數值模擬結果調整具體值。采用isobaric修正技術后，對需要進行速度合成的點也相應地擴大范圍。

2.3 Lagrange域壓力邊界條件的確定

對Lagrange域的計算來說，在E-L界面上施加壓力邊界條件，體現來自Euler域的作用。GEL中，這實際上是對于組成E-L界面上每條線段(二維情況下)的兩個節點，受到的來自Euler域壓力值的確定。文獻［1-2］中由內插得到線段中點受到的壓力值，作為整條線段受到的平均壓力值。而本文中分別計算了界面線段跨越的所有Euler網格對線段各部分施加的壓力值，將這些壓力值累加后再平均分配給線段的兩個端點，即［8］

式中:m為每條E-L界面線段跨越的Euler網格數總和，pj為Euler網格的壓力，lj是線段在Euler網格中的長度。

2.4 統一計算時間步長

時間步長［2］

式中:0.5為CFL系數，ΔtE和ΔtL分別是Euler和Lagrange計算得出的時間步長。最后，Euler域和Lagrange域均以此步長分別進行獨立地計算。

2.5 多界面的處理

GEL對兩種計算格式在界面的處理，對Euler計算來說，實際上就是對ghost域上各物理量的定義;而對Lagrange計算，是E-L界面上節點壓力邊界條件的計算。這些耦合的過程相對于他們各自的計算，是完全獨立的，因此在多界面問題中，界面與界面的處理過程，或者說每一個在界面上發生的Euler和Lagrange耦合過程是相互獨立的。因此，GEL擴展到多界面問題是完全直接的。

3 數值算例

考慮一維理想氣體和水的黎曼問題，計算區域為［0 m，3 m］。流場左右邊界均為流動邊界，左邊流場和右邊均為初始是靜止狀態的高壓雙原子理想氣體，壓力p=7.81 GPa，密度ρ=0.6 g/cm3，多方指數γ=1.4，這兩部分的計算采用以SCB格式編制的計算程序。流場中間為水，初始也是靜止狀態，密度ρ=1 g/cm3，壓力 p=0.1 MPa，此部分流場的計算采用 DEFEL二維動力有限元程序，狀態方程為Grüneisen狀態方程。初始左右兩界面分別在x=0.5，2.5 m處。左右兩邊流場的網格數為100，而中間部分為400，以本文中GEL耦合程序進行耦合。

圖2是不進行isobaric修正時的密度、速度和壓力結果。可以看到:當t=0.16 ms時，左右兩方向的激波傳至Lagrange域，但未相遇;當t=0.32 ms時，兩激波已經過碰撞后進行反射;當t=0.43 ms時，兩激波穿過界面到達左右兩邊的Euler域。從3個計算時間的結果看，只有當t=0.43 ms時左右界面處均出現明顯的overheating現象，尤其壓力較為明顯，引起的相對誤差為約8.6%。

取isobaric修正值Φ≥－0.36Δx，在t=0.43 ms與不進行修正計算結果的比較見圖3，最后為壓力圖像在左邊界面處的放大圖。此時overheating現象有效減小，引起的相對誤差為約2.6%。

圖2 無isobaric修正時的結果Fig.2 The results with no isobaric fix technique

圖3 Isobaric修正后與無isobaric修正時的結果比較(t=0.43 ms)Fig.3 The resultswith isobaric fix technique are compared against the one with no fix in t=0.43 ms

4 討論

GEL方法從理論和計算上將整個計算流域分為Euler和Lagrange域，結合Euler和Lagrange兩種計算方法的優越性;處理混合網格借助ghost-fluid的思想，程序編制相對簡單，計算量小;在耦合過程中，對每一界面的處理是獨立的，因此擴展到雙界面甚至多界面問題是完全直接的。本文中GEL耦合程序中，對Lagrange域在界面上施加壓力邊界條件的方法和以往文獻給出的方法不同。模擬了雙界面一維理想氣體和水的Riemann問題，以及有機玻璃圓柱在沖擊波作用下變形的二維問題。

從一維黎曼問題算例計算結果看到，當計算時間為0.43 ms時，即波陣面從Lagrange域穿過界面傳播到Euler域時，界面處overheating現象明顯比前兩時刻增大，而取其他一些isobaric修正量進行計算時也出現了相同情況。這是否由于兩種計算格式對耦合的不同處理，以及外插方向與波行進方向不一致而造成，仍然需要進行進一步的探討。

［1］Morano E，ArientiM，Hung P，etal.A level setapproach to Eulerian-Lagrangian coupling［J］.Journal of Computational Physics，2003，185:213-251.

［2］Fedkiw R P.Coupling an Eulerian fluid calculation to a Lagrangian solid calculation with the Ghost fluid method［J］.Journal of Computational Physics，2002，175:200-224.

［3］Fedkiw R，Aslam T，Merriman B，et al.A non-oscillatory Eulerian approach to interfaces in multimaterial flows(The Ghost fluid method)［J］.Journal of Computational Physics，1999，152:457-492.

［4］Noh W.CEL:A time-dependent，two space dimensional，coupled Eulerian-Lagrangian code［C］∥Methods in computational physics.New York:Academic Press，1964:117-179.

［5］Hancock S.PISCES 2DELK theoretical Manual［Z］.Physics International，1985.

［6］McMasterW H.Computercodes for fluid-structure interactions［C］∥Proc 1984 Pressure Vessel and Piping Conference.San Antonio，1984.

［7］Fedkiw R，Marquina A，Merriman B.An osobaric fix for the overheating problem inmultimaterial compressible flows［J］.Journal of Computational Physics，1999，48:545-578.

［8］姚陽，李平，柏勁松，等.Euler-Lagrange耦合計算的 GEL方法［J］.爆炸與沖擊，2007，27(5):420-425.YAO Yang，LIPing，BAIJin-song，et al.Euler-Lagrange coupling with GEL［J］.Explosion and Shock Waves，2007，27(5):420-425.

［9］Zhao N.High resolution SCB scheme for hyperbolic systems of2-D conservation laws［J］.Journal of Computational Physics，1998，16(2):179-192.

［10］DEFEL Users Manual［Z］.Dyna East Corporration，1988.