童石磊，鐘 敏，柏勁松，陳森華

(中國工程物理研究院流體物理研究所，四川綿陽 621999)

1 引 言

多介質(zhì)界面的數(shù)值模擬和界面運動追蹤是近年來計算流體力學(xué)研究的熱點問題，在材料科學(xué)、天體物理以及核物理等領(lǐng)域也有著廣泛應(yīng)用。目前，處理多介質(zhì)界面的數(shù)值模擬方法有Level-Set方法[1]、VOF方法[2]和Front Tracking (FT)方法[3]等。其中，F(xiàn)T方法是由James Glimm等人于20世紀(jì)80年代提出的，主要思想是通過參數(shù)表示的閉合曲線表征界面位置，界面隨著表征曲線的運動而運動。一般而言，這組閉合曲線由一些標(biāo)記點(示蹤點)組成，通過記錄示蹤點的運動即可捕捉到界面的運動情況。FT方法的優(yōu)點在于能準(zhǔn)確地捕捉界面的位置，使用了準(zhǔn)確的間斷條件，不引入數(shù)值擴(kuò)散，避免原本銳利的界面變得平滑；缺點在于算法比較復(fù)雜，難以處理拓?fù)浣Y(jié)構(gòu)變化，難以向高維推廣。Glimm及其團(tuán)隊對FT方法做了許多研究工作，取得了不少有價值的成果[4]，例如二維和三維的復(fù)雜數(shù)據(jù)存儲和處理方法、界面拓?fù)浠靵y的避免和糾正方法等。為避免和糾正拓?fù)浠靵y，Glimm等人[5]采用了Grid Free (GF)、Grid Based (GB)、GF和GB混雜以及Local Grid Based (LGB)等方法。GF方法的主要特點是不依賴輔助坐標(biāo)系即可刪除產(chǎn)生界面拓?fù)溴e誤的部分，重新調(diào)整網(wǎng)格分布和尺寸；GB方法則是引進(jìn)輔助坐標(biāo)系來完成這項任務(wù)。GF方法的數(shù)值擴(kuò)散較小，調(diào)整后的三角形元素形狀和尺寸好，計算精度高，但是不能處理拓?fù)湫再|(zhì)復(fù)雜的問題；而GB方法的數(shù)值擴(kuò)散大，生成的三角形元素質(zhì)量差，但是可以改善對復(fù)雜拓?fù)鋯栴}的處理效果?；祀s方法是周期性地交替使用GF和GB方法，平衡了兩種方法的優(yōu)缺點。LGB方法則是限制GB方法的使用范圍。

處理界面拓?fù)湫再|(zhì)復(fù)雜的問題時，F(xiàn)T方法出現(xiàn)困難的主要原因在于缺少對邊界處理的方法。本研究在經(jīng)典FT基礎(chǔ)上，將邊界上的拓?fù)浣Y(jié)構(gòu)變化通過示蹤點組成的線段進(jìn)行描述，邊界運動情況通過示蹤點線段的添加、合并、刪減自然地表達(dá)，產(chǎn)生的新界面則由新的示蹤點線段組成。通過這些改進(jìn)建立一種能夠有效處理界面拓?fù)渥兓姆椒?，以期為處理?fù)雜界面問題提供參考。

2 數(shù)值分析與數(shù)據(jù)結(jié)構(gòu)

對于FT方法，界面函數(shù)f上的示蹤點坐標(biāo)滿足

(1)

界面發(fā)展規(guī)律可以通過求解(2)式確定

(2)

式中：x為示蹤點位移，x0為示蹤點的初始位置，v為示蹤點的速度，t為時間。

2.1 傳統(tǒng)FT方法

傳統(tǒng)FT方法將示蹤點依次編號，并按編號順序?qū)⑹聚欬c相連組成封閉界面。對于二維問題，通常以二維數(shù)組Tx(N,L)和Ty(N,L)表示示蹤點坐標(biāo)，N為材料號，L為示蹤點編號。圖1為傳統(tǒng)FT方法示蹤點的分布，其中數(shù)字為不同的示蹤點編號，A物質(zhì)位于示蹤點走向的左側(cè)。

處理多物質(zhì)界面相交的問題時，傳統(tǒng)FT方法很容易在交界面處產(chǎn)生拓?fù)浠靵y。如圖2所示，曲線界面走向的左側(cè)為A物質(zhì)，直線界面走向的左側(cè)為B物質(zhì)，1、2、3、4、5為曲線界面上示蹤點。采用傳統(tǒng)FT方法計算時，由于無法保證界面相交之前恰好處于剛剛接觸狀態(tài)，因此，在下一個計算時間內(nèi)，很可能出現(xiàn)界面直接嵌入的情況(如圖2中虛線所示，示蹤點3、4、5走向的左側(cè)既含有A物質(zhì)，也含有B物質(zhì))，從而導(dǎo)致拓?fù)浠靵y，使計算無法繼續(xù)進(jìn)行。拓?fù)浠靵y的出現(xiàn)限制了傳統(tǒng)FT方法的應(yīng)用，針對該問題，Glimm等人對傳統(tǒng)FT方法進(jìn)行了改進(jìn)，使其能夠處理界面拓?fù)浣Y(jié)構(gòu)的變化，但是該方法仍然存在很大的局限性，并且要求問題中不能超過3種物質(zhì)。此外，為處理界面拓?fù)浣Y(jié)構(gòu)變化，在界面接觸之前，需要對界面示蹤點進(jìn)行一次重新分布，界面相交后，每運行一個時間步，就需要對示蹤點進(jìn)行一次重新分布，計算過程十分繁瑣，精度也不高。

圖1 傳統(tǒng)FT方法示蹤點分布Fig.1 Distribution of tracer points in classic FT method

圖2 拓?fù)浠靵y情況Fig.2 Topological chaos

2.2 改進(jìn)的FT方法

如圖3所示，改進(jìn)的FT方法用示蹤點線段描述界面，采用數(shù)組T(s,e,ml,mr)表示示蹤點線段，并對線段進(jìn)行編號，s表示線段的起點坐標(biāo)，e表示線段的終點坐標(biāo)，ml表示線段左側(cè)物質(zhì)的編號，mr表示線段右側(cè)物質(zhì)的編號。當(dāng)不同編號的線段具有相同的端點坐標(biāo)時，認(rèn)為線段相連，并由相連線段組成封閉界面。改進(jìn)的FT方法通過數(shù)組T即可清楚描述界面相關(guān)量，因此線段編號可以任意給定。

針對界面相交時拓?fù)浠靵y產(chǎn)生的問題，在已求得界面?zhèn)鞑ニ俣鹊臈l件下，改進(jìn)的FT方法通過設(shè)定界面Ⅰ上示蹤點到達(dá)界面Ⅱ的最短時間，來保證界面結(jié)構(gòu)發(fā)生拓?fù)渥兓扒『锰幱诮佑|狀態(tài)，具體方法如下：首先，通過數(shù)組T確定初始時刻的界面位置，再利用界面?zhèn)鞑ニ俣惹蟮媒缑姊裆系氖聚欬c到達(dá)界面Ⅱ位置所需的時間，選取其中最短的時間作為界面相交的判定條件，即可保證拓?fù)渥兓敖缑嫣幱谇『孟嘟佑|的狀態(tài)，避免拓?fù)浠靵y的產(chǎn)生。

圖3 改進(jìn)FT方法示蹤點線段的分布Fig.3 Distribution of segments using improved FT method

圖4 界面恰好接觸Fig.4 Interface just in contact

改進(jìn)的FT方法對界面拓?fù)浣Y(jié)構(gòu)變化的處理主要通過示蹤點線段實現(xiàn)。以兩界面相交問題為例，改進(jìn)的FT方法對界面拓?fù)浣Y(jié)構(gòu)變化問題的處理如圖5所示，其中共包含3種物質(zhì)，編號分別為0、1、2。選取部分界面進(jìn)行考察，設(shè)定最小距離d，當(dāng)示蹤點間距離小于d時，認(rèn)為示蹤點重合為一點；當(dāng)示蹤點與示蹤點線段間距離小于d時，認(rèn)為示蹤點與線段融合，并將線段分成兩部分。圖5(a)中，示蹤點線段A1B1與A2B2兩側(cè)的物質(zhì)不同，當(dāng)兩示蹤點線段相互接近至距離小于d時，線段融合產(chǎn)生新的界面A0B0；圖5(b)中，線段A1B1與A2B2發(fā)生融合后，得到的新線段兩側(cè)均為物質(zhì)1，因此新線段可以直接消去，完成界面拓?fù)浣Y(jié)構(gòu)的變化。

改進(jìn)的FT方法處理界面拓?fù)浣Y(jié)構(gòu)變化問題的計算流程如圖6所示。通過以上分析可以發(fā)現(xiàn)，改進(jìn)FT方法采用示蹤點線段對界面進(jìn)行處理十分簡便，由此我們可以建立新的界面拓?fù)浣Y(jié)構(gòu)變化追蹤和模擬機(jī)制，并進(jìn)行計算程序編寫。

圖5 界面拓?fù)浣Y(jié)構(gòu)變化Fig.5 Topological changes of interface

圖6 改進(jìn)的FT方法處理界面拓?fù)浣Y(jié)構(gòu)問題流程圖Fig.6 Flowchart of the improved FT method

3 數(shù)值算例

采取二維算例檢驗改進(jìn)FT方法的可靠性，以下物理量均為無量綱量。算例采用一階格式進(jìn)行離散并設(shè)定誤差允許區(qū)間，計算結(jié)果滿足精度要求，在以后的應(yīng)用中可以采用更高階離散格式得到更高精度的結(jié)果。由于兩算例在界面變化過程中均沒有進(jìn)行質(zhì)量輸運，因此滿足質(zhì)量守恒。

3.1 圓柱貫穿靶

本算例主要用于驗證改進(jìn)的FT方法能否有效處理界面拓?fù)涞淖兓D7為圓柱物質(zhì)A貫穿矩形靶板物質(zhì)B的過程，其中物質(zhì)A為剛體，物質(zhì)B為無強度、質(zhì)量為零的理想可變形體。該模型的特點是排除了狀態(tài)方程、本構(gòu)關(guān)系等真實物理參數(shù)的測量誤差對計算結(jié)果的影響。C物質(zhì)為真空，計算區(qū)域為[-200,200]×[-150,300]，計算時間為72，速度條件為

(3)

對貫穿過程進(jìn)行分析可以發(fā)現(xiàn)，界面拓?fù)渥兓譃橐韵?個過程。初始狀態(tài)下，界面如圖7(a)所示，A物質(zhì)與B物質(zhì)相互獨立，有各自的邊界；當(dāng)t=24時，界面位置如圖7(b)所示，A物質(zhì)穿入B物質(zhì)，兩物質(zhì)之間有一條公共邊界，拓?fù)湫再|(zhì)發(fā)生第1次變化；當(dāng)t=48時，界面位置如圖7(c)所示，A物質(zhì)開始部分穿出B物質(zhì)，B物質(zhì)被分為兩部分，原來的公共邊界消失并產(chǎn)生了兩條新的公共邊界，拓?fù)湫再|(zhì)發(fā)生第2次變化；當(dāng)t=72時，界面位置如圖7(d)所示，A物質(zhì)完全穿出B物質(zhì)，公共邊界消失，B物質(zhì)被分解為兩個完全獨立的部分，拓?fù)湫再|(zhì)發(fā)生第3次變化，計算終止。

由圖7圓柱貫穿靶過程中界面的變形情況可以看出，改進(jìn)的FT方法在處理界面拓?fù)浣Y(jié)構(gòu)變化時繼承了傳統(tǒng)方法的優(yōu)點，物質(zhì)界面運動情況能夠得到清晰的追蹤，界面形狀得到了較好的保持，同時保證了界面尖銳部分不會變得圓滑。

圖7 圓柱貫穿靶過程中的界面變形情況Fig.7 Deformation of cylinder-target interface in penetration process

3.2 正向和反向速度場中界面還原測試

選取矩形在速度場中先正向旋轉(zhuǎn)一段時間后，再反向旋轉(zhuǎn)相同時間回到初始位置作為算例，同時采用改進(jìn)的FT方法和Level-Set方法對該算例進(jìn)行計算，通過比較兩種方法計算所得旋轉(zhuǎn)前、后矩形邊界的變化情況，驗證改進(jìn)的FT方法的精確度。

矩形界面初始位置如圖8所示，計算區(qū)域為[0,π]×[0,π]，矩形界面所圍的區(qū)域為[0.3π,0.7π]×[0.05π,0.45π]。令u、v分別為x和y方向速度分量，當(dāng)t

(4)

當(dāng)t≥t*時，給定反向速度場為

(5)

式中：t*為速度場反向時刻。

圖8 初始時刻界面位置Fig.8 Position of interface at initial time

通過理論分析可知：在t=t*時刻，界面變形最大；在t=2t*時刻，界面應(yīng)當(dāng)恢復(fù)到初始時刻的形狀。在t*=8、t*=16兩種條件下，采用改進(jìn)的FT方法對矩形界面的演化情況進(jìn)行數(shù)值計算，得到0、t*、2t*時刻界面的形狀，如圖9所示。

由于t*時刻速度場反向，再旋轉(zhuǎn)t*后停止，此時界面恰好還原。因此，通過對計算停止時的界面形狀與初始時刻矩形的比較，便能夠反映出算法的準(zhǔn)確性。從圖9可以看出，旋轉(zhuǎn)后界面位置與初始時刻位置基本重合，旋轉(zhuǎn)過程中，界面形狀得到了很好的保持，尤其是界面尖銳部分不會被抹平或者丟失。由此可見，改進(jìn)的FT方法能夠得到精度較高的模擬結(jié)果。

作為對比，采用Level-Set方法對相同條件下矩形界面的演化情況進(jìn)行數(shù)值計算。計算采用二階離散格式，在界面尖銳處采用自適應(yīng)網(wǎng)格加密，每種情況分別給出0、t*、2t*時刻的界面形狀，結(jié)果如圖10所示。通過對圖9和圖10中結(jié)果進(jìn)行比較可以看出，隨著計算時間增加，改進(jìn)FT方法對界面保持更好，精確度更高。此外，由于改進(jìn)的FT方法只采用了一階離散格式，并未采用自適應(yīng)網(wǎng)格，因此，改進(jìn)的FT方法耗費的資源更少，并且有進(jìn)一步提高精度的空間。

圖9 改進(jìn)FT方法得到的不同時刻界面形狀Fig.9 Interface shapes at different time using improved FT method

圖10 Level-Set方法得到的不同時刻界面形狀Fig.10 Interface shapes at different time using Level-Set method

4 結(jié) 論

提出了一種新的數(shù)據(jù)結(jié)構(gòu)處理方法，建立了界面拓?fù)浣Y(jié)構(gòu)變化的追蹤和模擬機(jī)制。該方法不但繼承了傳統(tǒng)FT方法能夠準(zhǔn)確模擬界面運動的優(yōu)點，而且能夠處理復(fù)雜界面拓?fù)浣Y(jié)構(gòu)的變化。在二維情況下，對2個算例數(shù)值模擬的結(jié)果顯示，改進(jìn)的FT方法能夠準(zhǔn)確地處理多物質(zhì)邊界拓?fù)浣Y(jié)構(gòu)的變化。

[1] OSHER S,SETHIAN J A.Fronts propagating with curvature-dependent speed:algorithms based on Hamilton-Jacobi formulations [J].J Comput Phys,1988,79(1):12-49.

[2] HIRT C W,NICHOLS B D.Volume of fluid (VOF) method for the dynamics of free boundaries [J].J Comput Phys,1981,39(1):201-225.

[3] CHERN I L,GLIMM J,MCBRYAN O,et al.Front tracking for gas dynamics [J].J Comput Phys,1986,62(1):83-110.

[4] GLIMM J,GROVE J W,LI X L,et al.Three-dimensional front tracking [J].SIAM J Sci Comput,1998,19(3):703-727.

[5] GLIMM J,GROVE J W,LI X L,et al.Robust computational algorithms for dynamic interface tracking in three dimensions [J].SIAM J Sci Comput,2000,21(6):2240-2256.