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

姚 陽(yáng)

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

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

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

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

1 基本方程

1.1 Euler方程

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

式中:ρ為密度，u為速度矢量，p為壓力，e為質(zhì)量?jī)?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二維動(dòng)力有限元程序［10］。

2 GEL方法

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

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

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

首先通過(guò)程函數(shù)

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

Ghost點(diǎn)的速度的切向分量來(lái)自外插速度的切向分量，法向分量還需由Lagrange計(jì)算的E-L界面上點(diǎn)的法向速度確定。以二維情況為例，此E-L界面上的點(diǎn)應(yīng)是距ghost點(diǎn)最近的E-L界面線段上的點(diǎn)，速度值通過(guò)此點(diǎn)所在的界面線段兩節(jié)點(diǎn)的速度值線性?xú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個(gè)獨(dú)立的物理量。因?yàn)樵谟屑げê途哂泄腆w不可穿透性質(zhì)的界面相互作用的數(shù)值模擬中，界面附近往往會(huì)出現(xiàn)所謂的overheating誤差，而壓力和速度在界面保持一致。為削弱此誤差，除了用isobaric修正變量，即熵、T'(ρ)和e'(ρ)作為第3個(gè)外插量外，還可以用isobaric修正技術(shù)將誤差最小化［7］。在GEL方法中，外插時(shí)凡是滿(mǎn)足

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

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

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

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

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

時(shí)間步長(zhǎng)［2］

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

2.5 多界面的處理

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

3 數(shù)值算例

考慮一維理想氣體和水的黎曼問(wèn)題，計(jì)算區(qū)域?yàn)椋? m，3 m］。流場(chǎng)左右邊界均為流動(dòng)邊界，左邊流場(chǎng)和右邊均為初始是靜止?fàn)顟B(tài)的高壓雙原子理想氣體，壓力p=7.81 GPa，密度ρ=0.6 g/cm3，多方指數(shù)γ=1.4，這兩部分的計(jì)算采用以SCB格式編制的計(jì)算程序。流場(chǎng)中間為水，初始也是靜止?fàn)顟B(tài)，密度ρ=1 g/cm3，壓力 p=0.1 MPa，此部分流場(chǎng)的計(jì)算采用 DEFEL二維動(dòng)力有限元程序，狀態(tài)方程為Grüneisen狀態(tài)方程。初始左右兩界面分別在x=0.5，2.5 m處。左右兩邊流場(chǎng)的網(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)過(guò)碰撞后進(jìn)行反射;當(dāng)t=0.43 ms時(shí)，兩激波穿過(guò)界面到達(dá)左右兩邊的Euler域。從3個(gè)計(jì)算時(shí)間的結(jié)果看，只有當(dāng)t=0.43 ms時(shí)左右界面處均出現(xiàn)明顯的overheating現(xiàn)象，尤其壓力較為明顯，引起的相對(duì)誤差為約8.6%。

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

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

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

從一維黎曼問(wèn)題算例計(jì)算結(jié)果看到，當(dāng)計(jì)算時(shí)間為0.43 ms時(shí)，即波陣面從Lagrange域穿過(guò)界面?zhèn)鞑サ紼uler域時(shí)，界面處overheating現(xiàn)象明顯比前兩時(shí)刻增大，而取其他一些isobaric修正量進(jìn)行計(jì)算時(shí)也出現(xiàn)了相同情況。這是否由于兩種計(jì)算格式對(duì)耦合的不同處理，以及外插方向與波行進(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］姚陽(yáng)，李平，柏勁松，等.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.