GEL方法在一維雙界面問題中的應(yīng)用*

姚 陽

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

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

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

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

1 基本方程

1.1 Euler方程

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

式中:ρ為密度，u為速度矢量，p為壓力，e為質(zhì)量內(nèi)能，t為時(shí)間，I是單位矩陣。狀態(tài)方程形式取為

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

1.2 Lagrange方程

連續(xù)介質(zhì)力學(xué)給出的基本方程有

式中:ρ、u和E分別表示質(zhì)點(diǎn)密度、速度和體積內(nèi)能，~b為體積體力，~σ為Cauchy應(yīng)力張量。當(dāng)采用流體彈塑性本構(gòu)模型時(shí)，應(yīng)力張量分解為

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

2 GEL方法

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

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

GEL耦合是根據(jù)E-L界面接觸間斷的性質(zhì)，界面兩側(cè)應(yīng)滿足的物理量條件，在兩個計(jì)算區(qū)域分別確定界面上連續(xù)或間斷物理量的過程中，內(nèi)邊界條件的確定。對E-L界面上的切向速度、密度(或熵)，在界面兩邊區(qū)域上的這些間斷量互不相關(guān)，因此兩個計(jì)算區(qū)域分別由各自區(qū)域計(jì)算結(jié)果確定E-L界面上的值。而由界面兩側(cè)法向速度的連續(xù)性，ghost點(diǎn)的法向速度由Lagrange計(jì)算的E-L界面上的法向速度確定，即相當(dāng)于對Euler域施加了速度邊界條件。利用壓力是連續(xù)量，以E-L界面附近Euler網(wǎng)格的壓力，在E-L界面可對Lagrange域施加壓力邊界條件。通過相互施加邊界條件，體現(xiàn)著Euler域和Lagrange域之間的相互作用。圖1是E-L界面連續(xù)量確定示意圖。耦合的過程需在每一時(shí)間步計(jì)算前進(jìn)行，然后以統(tǒng)一的時(shí)間步長用標(biāo)準(zhǔn)的Euler和Lagrange計(jì)算程序分別進(jìn)行獨(dú)立地計(jì)算。

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

2.1 Ghost域的定義

首先通過程函數(shù)

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

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

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

2.2 Isobaric修正技術(shù)

除壓力和速度外，還需要得到ghost網(wǎng)格上第3個獨(dú)立的物理量。因?yàn)樵谟屑げê途哂泄腆w不可穿透性質(zhì)的界面相互作用的數(shù)值模擬中，界面附近往往會出現(xiàn)所謂的overheating誤差，而壓力和速度在界面保持一致。為削弱此誤差，除了用isobaric修正變量，即熵、T'(ρ)和e'(ρ)作為第3個外插量外，還可以用isobaric修正技術(shù)將誤差最小化［7］。在GEL方法中，外插時(shí)凡是滿足

所有Euler域格點(diǎn)的物理量也隨之改變。式中:Δx為網(wǎng)格寬，k為無物理意義的非零常系數(shù)，可根據(jù)數(shù)值模擬結(jié)果調(diào)整具體值。采用isobaric修正技術(shù)后，對需要進(jìn)行速度合成的點(diǎn)也相應(yīng)地?cái)U(kuò)大范圍。

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

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

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

2.4 統(tǒng)一計(jì)算時(shí)間步長

時(shí)間步長［2］

式中:0.5為CFL系數(shù)，ΔtE和ΔtL分別是Euler和Lagrange計(jì)算得出的時(shí)間步長。最后，Euler域和Lagrange域均以此步長分別進(jìn)行獨(dú)立地計(jì)算。

2.5 多界面的處理

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

3 數(shù)值算例

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

圖2是不進(jìn)行isobaric修正時(shí)的密度、速度和壓力結(jié)果?？梢钥吹?當(dāng)t=0.16 ms時(shí)，左右兩方向的激波傳至Lagrange域，但未相遇;當(dāng)t=0.32 ms時(shí)，兩激波已經(jīng)過碰撞后進(jìn)行反射;當(dāng)t=0.43 ms時(shí)，兩激波穿過界面到達(dá)左右兩邊的Euler域。從3個計(jì)算時(shí)間的結(jié)果看，只有當(dāng)t=0.43 ms時(shí)左右界面處均出現(xiàn)明顯的overheating現(xiàn)象，尤其壓力較為明顯，引起的相對誤差為約8.6%。

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

圖2 無isobaric修正時(shí)的結(jié)果Fig.2 The results with no isobaric fix technique

圖3 Isobaric修正后與無isobaric修正時(shí)的結(jié)果比較(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方法從理論和計(jì)算上將整個計(jì)算流域分為Euler和Lagrange域，結(jié)合Euler和Lagrange兩種計(jì)算方法的優(yōu)越性;處理混合網(wǎng)格借助ghost-fluid的思想，程序編制相對簡單，計(jì)算量小;在耦合過程中，對每一界面的處理是獨(dú)立的，因此擴(kuò)展到雙界面甚至多界面問題是完全直接的。本文中GEL耦合程序中，對Lagrange域在界面上施加壓力邊界條件的方法和以往文獻(xiàn)給出的方法不同。模擬了雙界面一維理想氣體和水的Riemann問題，以及有機(jī)玻璃圓柱在沖擊波作用下變形的二維問題。

從一維黎曼問題算例計(jì)算結(jié)果看到，當(dāng)計(jì)算時(shí)間為0.43 ms時(shí)，即波陣面從Lagrange域穿過界面?zhèn)鞑サ紼uler域時(shí)，界面處overheating現(xiàn)象明顯比前兩時(shí)刻增大，而取其他一些isobaric修正量進(jìn)行計(jì)算時(shí)也出現(xiàn)了相同情況。這是否由于兩種計(jì)算格式對耦合的不同處理，以及外插方向與波行進(jìn)方向不一致而造成，仍然需要進(jìn)行進(jìn)一步的探討。

［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耦合計(jì)算的 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.