處理多介質界面的改進型Front Tracking方法*

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

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

1 引 言

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

處理界面拓撲性質復雜的問題時，FT方法出現困難的主要原因在于缺少對邊界處理的方法。本研究在經典FT基礎上，將邊界上的拓撲結構變化通過示蹤點組成的線段進行描述，邊界運動情況通過示蹤點線段的添加、合并、刪減自然地表達，產生的新界面則由新的示蹤點線段組成。通過這些改進建立一種能夠有效處理界面拓撲變化的方法，以期為處理復雜界面問題提供參考。

2 數值分析與數據結構

對于FT方法，界面函數f上的示蹤點坐標滿足

(1)

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

(2)

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

2.1 傳統FT方法

傳統FT方法將示蹤點依次編號，并按編號順序將示蹤點相連組成封閉界面。對于二維問題，通常以二維數組Tx(N,L)和Ty(N,L)表示示蹤點坐標，N為材料號，L為示蹤點編號。圖1為傳統FT方法示蹤點的分布，其中數字為不同的示蹤點編號，A物質位于示蹤點走向的左側。

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

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

圖2 拓撲混亂情況Fig.2 Topological chaos

2.2 改進的FT方法

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

針對界面相交時拓撲混亂產生的問題，在已求得界面傳播速度的條件下，改進的FT方法通過設定界面Ⅰ上示蹤點到達界面Ⅱ的最短時間，來保證界面結構發生拓撲變化前恰好處于接觸狀態，具體方法如下：首先，通過數組T確定初始時刻的界面位置，再利用界面傳播速度求得界面Ⅰ上的示蹤點到達界面Ⅱ位置所需的時間，選取其中最短的時間作為界面相交的判定條件，即可保證拓撲變化前界面處于恰好相接觸的狀態，避免拓撲混亂的產生。

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

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

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

改進的FT方法處理界面拓撲結構變化問題的計算流程如圖6所示。通過以上分析可以發現，改進FT方法采用示蹤點線段對界面進行處理十分簡便，由此我們可以建立新的界面拓撲結構變化追蹤和模擬機制，并進行計算程序編寫。

圖5 界面拓撲結構變化Fig.5 Topological changes of interface

圖6 改進的FT方法處理界面拓撲結構問題流程圖Fig.6 Flowchart of the improved FT method

3 數值算例

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

3.1 圓柱貫穿靶

本算例主要用于驗證改進的FT方法能否有效處理界面拓撲的變化。圖7為圓柱物質A貫穿矩形靶板物質B的過程，其中物質A為剛體，物質B為無強度、質量為零的理想可變形體。該模型的特點是排除了狀態方程、本構關系等真實物理參數的測量誤差對計算結果的影響。C物質為真空，計算區域為[-200,200]×[-150,300]，計算時間為72，速度條件為

(3)

對貫穿過程進行分析可以發現，界面拓撲變化分為以下4個過程。初始狀態下，界面如圖7(a)所示，A物質與B物質相互獨立，有各自的邊界；當t=24時，界面位置如圖7(b)所示，A物質穿入B物質，兩物質之間有一條公共邊界，拓撲性質發生第1次變化；當t=48時，界面位置如圖7(c)所示，A物質開始部分穿出B物質，B物質被分為兩部分，原來的公共邊界消失并產生了兩條新的公共邊界，拓撲性質發生第2次變化；當t=72時，界面位置如圖7(d)所示，A物質完全穿出B物質，公共邊界消失，B物質被分解為兩個完全獨立的部分，拓撲性質發生第3次變化，計算終止。

由圖7圓柱貫穿靶過程中界面的變形情況可以看出，改進的FT方法在處理界面拓撲結構變化時繼承了傳統方法的優點，物質界面運動情況能夠得到清晰的追蹤，界面形狀得到了較好的保持，同時保證了界面尖銳部分不會變得圓滑。

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

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

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

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

(4)

當t≥t*時，給定反向速度場為

(5)

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

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

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

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

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

圖9 改進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 結 論

提出了一種新的數據結構處理方法，建立了界面拓撲結構變化的追蹤和模擬機制。該方法不但繼承了傳統FT方法能夠準確模擬界面運動的優點，而且能夠處理復雜界面拓撲結構的變化。在二維情況下，對2個算例數值模擬的結果顯示，改進的FT方法能夠準確地處理多物質邊界拓撲結構的變化。

[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.