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

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

1 引 言

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

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

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

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

(1)

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

(2)

式中：x為示蹤點(diǎn)位移，x0為示蹤點(diǎn)的初始位置，v為示蹤點(diǎn)的速度，t為時(shí)間。

2.1 傳統(tǒng)FT方法

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

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

圖1 傳統(tǒng)FT方法示蹤點(diǎn)分布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方法用示蹤點(diǎn)線段描述界面，采用數(shù)組T(s,e,ml,mr)表示示蹤點(diǎn)線段，并對(duì)線段進(jìn)行編號(hào)，s表示線段的起點(diǎn)坐標(biāo)，e表示線段的終點(diǎn)坐標(biāo)，ml表示線段左側(cè)物質(zhì)的編號(hào)，mr表示線段右側(cè)物質(zhì)的編號(hào)。當(dāng)不同編號(hào)的線段具有相同的端點(diǎn)坐標(biāo)時(shí)，認(rèn)為線段相連，并由相連線段組成封閉界面。改進(jìn)的FT方法通過(guò)數(shù)組T即可清楚描述界面相關(guān)量，因此線段編號(hào)可以任意給定。

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

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

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

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

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

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

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

3 數(shù)值算例

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

3.1 圓柱貫穿靶

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

(3)

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

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

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

3.2 正向和反向速度場(chǎng)中界面還原測(cè)試

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

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

(4)

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

(5)

式中：t*為速度場(chǎng)反向時(shí)刻。

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

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

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

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

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

圖10 Level-Set方法得到的不同時(shí)刻界面形狀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ùn)動(dòng)的優(yōu)點(diǎn)，而且能夠處理復(fù)雜界面拓?fù)浣Y(jié)構(gòu)的變化。在二維情況下，對(duì)2個(gè)算例數(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.