基于八叉樹自適應(yīng)網(wǎng)格技術(shù)的Level Set運(yùn)動(dòng)界面追蹤方法

殷亞軍，李陽(yáng)東，涂志新，李文，沈旭，周建新

?

基于八叉樹自適應(yīng)網(wǎng)格技術(shù)的Level Set運(yùn)動(dòng)界面追蹤方法

殷亞軍，李陽(yáng)東，涂志新，李文，沈旭，周建新

（華中科技大學(xué)材料成形與模具技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北武漢430074）

Level Set方法因能有效地處理界面處復(fù)雜的拓?fù)浣Y(jié)構(gòu)變化以及大變形問(wèn)題，廣泛應(yīng)用于界面追蹤領(lǐng)域。在Level Set方法追蹤運(yùn)動(dòng)界面時(shí)引入八叉樹網(wǎng)格技術(shù)，通過(guò)八叉樹網(wǎng)格的細(xì)化和粗化技術(shù)減少計(jì)算網(wǎng)格數(shù)量和計(jì)算內(nèi)存并提高計(jì)算效率和計(jì)算精度。因?yàn)榘瞬鏄渚W(wǎng)格為非均勻網(wǎng)格，其相鄰網(wǎng)格的層數(shù)值可能不相同，所以不能直接采用WENO格式離散Level Set函數(shù)得到網(wǎng)格處的函數(shù)值，進(jìn)而提出八叉樹網(wǎng)格離散模型解決這一問(wèn)題，并提出基于八叉樹網(wǎng)格距離場(chǎng)重新初始化方法減少Level Set方法的質(zhì)量損失，最后將基于八叉樹網(wǎng)格技術(shù)的Level Set方法應(yīng)用于兩個(gè)給定速度場(chǎng)的運(yùn)動(dòng)界面模擬算例以及基準(zhǔn)件方腔的鑄造充型過(guò)程的模擬。模擬結(jié)果表明該方法可以提高界面的精度，同時(shí)改善質(zhì)量守恒性。

界面追蹤；數(shù)值模擬；八叉樹網(wǎng)格技術(shù)；Level Set方法；網(wǎng)格自適應(yīng)；模型；優(yōu)化

引 言

許多物理問(wèn)題會(huì)涉及氣液自由界面。采用數(shù)值模擬計(jì)算自由界面時(shí)，界面附近流體的物性會(huì)發(fā)生明顯變化，需要準(zhǔn)確追蹤自由界面的位置。Level Set（水平集）方法是美國(guó)計(jì)算數(shù)學(xué)學(xué)者Osher等[1-4]提出的一種追蹤界面的方法，該方法能夠有效地解決復(fù)雜界面追蹤問(wèn)題，是目前界面追蹤問(wèn)題數(shù)值模擬的主流方法。

Level Set方法可以克服一般界面前沿追蹤方法必須構(gòu)造出具體的界面而難以處理界面前沿發(fā)生的拓?fù)浣Y(jié)構(gòu)變化的弱點(diǎn)，不需要顯式追蹤運(yùn)動(dòng)界面，從而可以較容易地處理復(fù)雜的物質(zhì)界面及其拓?fù)浣Y(jié)構(gòu)發(fā)生變化的情形，便于計(jì)算曲率及界面法向量，可以將間斷問(wèn)題轉(zhuǎn)化為光滑問(wèn)題加以解決。與VOF（volume of fluid）法相比，該方法易于從二維問(wèn)題推廣到三維問(wèn)題，追蹤界面精度較高。

但是，經(jīng)典的Level Set方法守恒性差，初始時(shí)刻Level Set為距離函數(shù)，經(jīng)過(guò)有限的時(shí)間步長(zhǎng)后不再保持距離函數(shù)的性質(zhì)，需要重新初始化，重新初始化Level Set函數(shù)涉及的復(fù)雜數(shù)學(xué)過(guò)程極大地增加了計(jì)算量。此外，為了減少數(shù)值耗散，在求解Level Set方程時(shí)通常采用高階格式離散空間導(dǎo)數(shù)和時(shí)間導(dǎo)數(shù)，這些高階計(jì)算格式會(huì)影響其計(jì)算效率。為了提高Level Set方法的計(jì)算速率與精度，許多學(xué)者對(duì)此做了大量的研究。Osher等[5]對(duì)Level Set重新初始化方法進(jìn)行了研究。Sethian等提出了Fast Marching Method (FMM，快速步進(jìn)法)[6-7]對(duì)Level Set函數(shù)進(jìn)行重新初始化，提出窄帶方法[8]使Level Set方法更加快速有效。Peng等[9]在全局Level Set方法的基礎(chǔ)上提出了一種快速局部Level Set方法，大大節(jié)省了計(jì)算時(shí)間。Strain[10]首先引入樹型結(jié)構(gòu)的自適應(yīng)網(wǎng)格，在交界面附近采用高分辨率網(wǎng)格，而在遠(yuǎn)離交界面的區(qū)域保持較粗的網(wǎng)格。Sochnikov等[11]將Fast Marching（快速步進(jìn)）技術(shù)應(yīng)用到四叉樹型結(jié)構(gòu)的自適應(yīng)網(wǎng)格中，使得計(jì)算量進(jìn)一步減少。Costa等[12]采用了自適應(yīng)網(wǎng)格加密（adaptive mesh refinement, AMR）技術(shù)進(jìn)行結(jié)構(gòu)拓?fù)鋬?yōu)化。宮翔飛等[13]將Level Set方法與AMR技術(shù)融合用于數(shù)值模擬。黃筱云等提出了一種基于樹型結(jié)構(gòu)空間自適應(yīng)的快速粒子Level Set 方法[14]，并提出自適應(yīng)四叉樹網(wǎng)格下的N-S方程數(shù)值求解模型[15]。王生輝等[16]將四叉樹直角坐標(biāo)網(wǎng)格技術(shù)與Level Set方法相結(jié)合進(jìn)行運(yùn)動(dòng)界面的追蹤。Losasso等[17]利用八叉樹網(wǎng)格技術(shù)模擬了煙霧和水流的運(yùn)動(dòng)。朱雷等[18]采用八叉樹網(wǎng)格模擬焊接溫度場(chǎng)，對(duì)比了八叉樹網(wǎng)格與均勻網(wǎng)格的計(jì)算效率，采用八叉樹網(wǎng)格的計(jì)算時(shí)間只有均勻網(wǎng)格的40%。

為了進(jìn)一步提高Level Set方法的計(jì)算效率以及界面精度，本研究將八叉樹自適應(yīng)網(wǎng)格技術(shù)應(yīng)用于Level Set方法的求解，提出適用于Level Set方法追蹤界面的八叉樹自適應(yīng)網(wǎng)格模型，并建立了基于八叉樹自適應(yīng)網(wǎng)格的Level Set函數(shù)的離散和重新初始化模型。采用八叉樹網(wǎng)格容易進(jìn)行網(wǎng)格的局部加密或粗化，減少網(wǎng)格數(shù)量以及計(jì)算量。同時(shí)，局部網(wǎng)格細(xì)化能夠很好地提高界面追蹤的準(zhǔn)確性。本研究增加界面處網(wǎng)格的最大層數(shù)值，可以很大程度地提高界面追蹤的精度以及減少Level Set方法引起的數(shù)值耗散，提高界面追蹤的質(zhì)量守恒性，這對(duì)于改善流場(chǎng)失真具有重大的意義。

1 八叉樹網(wǎng)格技術(shù)

八叉樹是一種用于描述三維空間的樹狀數(shù)據(jù)結(jié)構(gòu)，其原理是將空間區(qū)域不斷分解為8個(gè)同樣大小的子區(qū)域，分解的次數(shù)越多，子區(qū)域就越小，一直到同一區(qū)域的屬性單一為止[19]。圖1所示即為一個(gè)簡(jiǎn)單的八叉樹節(jié)點(diǎn)產(chǎn)生過(guò)程。圖中左側(cè)為對(duì)應(yīng)的空間區(qū)域劃分，不同的區(qū)域用不同顏色的正方體表示；右側(cè)為八叉樹樹型結(jié)構(gòu)對(duì)應(yīng)的網(wǎng)格節(jié)點(diǎn)關(guān)系，圖中黑色節(jié)點(diǎn)表示子區(qū)域的父區(qū)域，葉子表示八叉樹網(wǎng)格的葉子單元網(wǎng)格。八叉樹網(wǎng)格以八叉樹作為基本理論基礎(chǔ)，網(wǎng)格組織方式與之類似，樹中指代的區(qū)域即為數(shù)值模擬中對(duì)應(yīng)的網(wǎng)格（網(wǎng)格單元統(tǒng)一為正方體單元）。

2 八叉樹網(wǎng)格自適應(yīng)

Level Set 方法追蹤運(yùn)動(dòng)界面時(shí)，界面精度和質(zhì)量損失都與網(wǎng)格的分辨率有關(guān)，如果將網(wǎng)格全部細(xì)化，則會(huì)很大程度地增加計(jì)算時(shí)間和存儲(chǔ)空間。需要自適應(yīng)網(wǎng)格技術(shù)對(duì)界面附近的區(qū)域進(jìn)行網(wǎng)格細(xì)化處理，對(duì)遠(yuǎn)離界面的區(qū)域進(jìn)行網(wǎng)格粗化處理，這樣既可以保證界面附近網(wǎng)格有較高的分辨率，提高界面區(qū)域的計(jì)算精度，又可以減少遠(yuǎn)離界面區(qū)域所占用的內(nèi)存，減少計(jì)算時(shí)間和計(jì)算內(nèi)存。

2.1 細(xì)化模型

采用八叉樹自適應(yīng)網(wǎng)格技術(shù)可以提高界面附近區(qū)域的網(wǎng)格分辨率，保證界面附近區(qū)域的計(jì)算精度。圖2所示即為一個(gè)簡(jiǎn)單的網(wǎng)格細(xì)化實(shí)例，其中假設(shè)白色區(qū)域?yàn)檫\(yùn)動(dòng)界面，而如果網(wǎng)格如圖左側(cè)所示，則需將界面周圍網(wǎng)格細(xì)化（右圖中紅色網(wǎng)格），轉(zhuǎn)化為右側(cè)圖所示，才能滿足計(jì)算分析要求。

2.2 粗化模型

八叉樹自適應(yīng)網(wǎng)格與均勻網(wǎng)格的優(yōu)勢(shì)就在于效率高、計(jì)算網(wǎng)格數(shù)少，許多局部計(jì)算區(qū)域可以用粗化網(wǎng)格代替眾多細(xì)化網(wǎng)格。圖3即為某一計(jì)算時(shí)刻網(wǎng)格粗化的例子，白色環(huán)形區(qū)域?yàn)檫\(yùn)動(dòng)界面，左邊部分網(wǎng)格不太理想，因?yàn)橹行募t色顯示部分網(wǎng)格沒(méi)有必要采用細(xì)化網(wǎng)格，必須經(jīng)過(guò)一定的粗化模型粗化，如圖3右邊部分，這樣就無(wú)須在該紅色網(wǎng)格區(qū)域浪費(fèi)大量的計(jì)算資源。

3 Level Set方程的離散

Level Set方法的主要思想是將隨時(shí)間運(yùn)動(dòng)的物質(zhì)界面定義為一個(gè)函數(shù)的零等值面（線），即(,)=0。在任意某一時(shí)刻，只要求出值，并求出其零等值面，就能夠知道此時(shí)的活動(dòng)界面位置。

在兩相流或自由界面追蹤問(wèn)題中，速度等物理量的控制方程是N-S方程，所對(duì)應(yīng)的Level Set方程見式（1），其中表示Level Set函數(shù)。

本研究采用五階WENO格式離散Level Set函數(shù)對(duì)空間的偏導(dǎo)數(shù)[20]，采用三階TVD Runge-Kutta格式離散Level Set函數(shù)對(duì)時(shí)間的偏導(dǎo)數(shù)[21]。

4 八叉樹網(wǎng)格離散模型

對(duì)于均勻網(wǎng)格來(lái)說(shuō)，采用上述離散方式直接離散Level Set函數(shù)。但是八叉樹自適應(yīng)網(wǎng)格為非均勻網(wǎng)格，每個(gè)網(wǎng)格的鄰居網(wǎng)格層數(shù)值（層數(shù)值表征該網(wǎng)格的細(xì)化或粗化程度，其中樹根節(jié)點(diǎn)網(wǎng)格的層數(shù)值為0，層數(shù)值越大網(wǎng)格單元的劃分越細(xì)）可能不同，盡管單元數(shù)據(jù)仍然存儲(chǔ)在正方體網(wǎng)格中心，但相鄰單元網(wǎng)格的大小不同，不能直接采用上述的離散方式，需要對(duì)非均勻網(wǎng)格重新進(jìn)行處理，進(jìn)而獲得鄰居網(wǎng)格的Level Set函數(shù)值進(jìn)行導(dǎo)數(shù)離散。

八叉樹網(wǎng)格與其相鄰網(wǎng)格之間存在3種關(guān)系：當(dāng)前網(wǎng)格與相鄰網(wǎng)格均為葉子單元，如圖4（a）所示；當(dāng)前網(wǎng)格為葉子單元，相鄰網(wǎng)格為子單元網(wǎng)格，如圖4（b）所示；當(dāng)前網(wǎng)格為子單元網(wǎng)格，相鄰網(wǎng)格為葉子單元，如圖4（c）所示。其中，紅色網(wǎng)格為當(dāng)前網(wǎng)格，綠色網(wǎng)格為相鄰網(wǎng)格。

對(duì)于圖4（a）所示情況，計(jì)算當(dāng)前網(wǎng)格Level Set函數(shù)的導(dǎo)數(shù)可以采用均勻網(wǎng)格計(jì)算方法得到，直接選取相鄰網(wǎng)格的值進(jìn)行計(jì)算。對(duì)于圖4（b）所示情況，在計(jì)算當(dāng)前網(wǎng)格Level Set函數(shù)的導(dǎo)數(shù)時(shí)首先計(jì)算相鄰網(wǎng)格值，在其相鄰網(wǎng)格的8個(gè)子單元中找出與當(dāng)前計(jì)算網(wǎng)格面相鄰的4個(gè)子單元，因?yàn)楫?dāng)前網(wǎng)格在該面的物理場(chǎng)變量的變化只與該面的流通量有關(guān)，故需要這4個(gè)單元才能獲得相鄰網(wǎng)格值。本研究采用一種較為簡(jiǎn)單的方式獲取相鄰網(wǎng)格值，將4個(gè)子單元對(duì)應(yīng)的中心網(wǎng)格值取平均值計(jì)算相鄰網(wǎng)格的Level Set函數(shù)值，然后用相鄰網(wǎng)格值計(jì)算當(dāng)前網(wǎng)格Level Set函數(shù)的空間導(dǎo)數(shù)。相鄰網(wǎng)格值的計(jì)算如式（2）所示

式中，+1表示導(dǎo)數(shù)計(jì)算時(shí)鄰居網(wǎng)格值，、、和分別表示與當(dāng)前計(jì)算網(wǎng)格面接觸的4個(gè)鄰居子單元網(wǎng)格值。

對(duì)于圖4（c）所示情況，計(jì)算子單元網(wǎng)格的Level Set函數(shù)值，需要借助一定的插值方法獲得鄰居網(wǎng)格值，將相鄰網(wǎng)格構(gòu)造出8個(gè)相同的虛單元網(wǎng)格（圖5中綠色網(wǎng)格），利用插值公式計(jì)算相鄰虛單元網(wǎng)格的函數(shù)值并作為鄰居網(wǎng)格值，然后用該虛單元的值計(jì)算當(dāng)前子單元網(wǎng)格的Level Set函數(shù)導(dǎo)數(shù)值。鄰居網(wǎng)格的插值公式如式（3）所示

式中，0+1表示相鄰虛單元網(wǎng)格值，+1表示相鄰葉子單元網(wǎng)格值，D、D、D分別表示虛子單元網(wǎng)格中心到相鄰葉子單元網(wǎng)格的距離。

圖5 相鄰虛單元網(wǎng)格

Fig.5 Illustration of neighbor virtural subunit grid

5 Level Set函數(shù)重新初始化

求解Level Set函數(shù)時(shí)，保持(,)為距離函數(shù)非常重要，但是Level Set方法守恒性差，經(jīng)過(guò)幾個(gè)計(jì)算時(shí)間步長(zhǎng)后(,)不再保持距離函數(shù)性質(zhì)，需要對(duì)時(shí)刻求解得到的(,)函數(shù)進(jìn)行修正，使函數(shù)的近似值滿足距離場(chǎng)的條件。Sussman等[3]提出了一種新的重新初始化方法

(5)

其中，(0)是函數(shù)符號(hào)。Peng等[9]采用式（5）所示格式可以得到更好的效果。本研究采用式（5）形式對(duì)Level Set函數(shù)進(jìn)行重新初始化。

重新初始化Level Set函數(shù)值會(huì)發(fā)生改變，所以八叉樹網(wǎng)格中的Level Set函數(shù)值也需要更新，八叉樹網(wǎng)格中距離場(chǎng)函數(shù)更新步驟如下：①將界面處網(wǎng)格單元（圖6中紫色網(wǎng)格）標(biāo)記為1，其余單元標(biāo)記為0，若邊界位置網(wǎng)格不是界面單元?jiǎng)t將其Level Set函數(shù)值設(shè)置為一個(gè)極大值；②從邊界位置網(wǎng)格[圖6（a）中綠色網(wǎng)格]掃描遍歷網(wǎng)格單元，當(dāng)掃描到細(xì)化網(wǎng)格單元如圖6（b）紅色圓圈中的網(wǎng)格[圖6（b）為圖6（a）中紅色圈出網(wǎng)格的局部放大圖]，本研究規(guī)定先掃描與當(dāng)前掃描網(wǎng)格相鄰的4個(gè)子單元網(wǎng)格，再掃描與這4個(gè)子單元網(wǎng)格相鄰的子單元網(wǎng)格，在未掃描到界面單元之前這些網(wǎng)格單元的距離場(chǎng)函數(shù)值不更新；③當(dāng)掃描到標(biāo)記為1的界面網(wǎng)格單元時(shí)，直接跳過(guò)該單元，不更新該網(wǎng)格單元的距離場(chǎng)函數(shù)值，繼續(xù)沿該方向掃描并更新后面標(biāo)記為0的網(wǎng)格單元的Level Set函數(shù)值，直至掃描到下一個(gè)邊界位置；④從已更新函數(shù)值的邊界網(wǎng)格[圖6（a）中綠色圓圈中的網(wǎng)格]的相反方向開始掃描，在掃描到界面網(wǎng)格之前不更新已掃描網(wǎng)格的函數(shù)值，當(dāng)掃描到界面網(wǎng)格單元時(shí)直接跳過(guò)該網(wǎng)格，繼續(xù)掃描網(wǎng)格單元并更新之后掃描網(wǎng)格單元的函數(shù)值，直至掃描到界面網(wǎng)格為止，如果當(dāng)前掃描網(wǎng)格的鄰居網(wǎng)格有多個(gè)界面網(wǎng)格則該網(wǎng)格的距離場(chǎng)函數(shù)值為3個(gè)方向鄰居網(wǎng)格的最小值；⑤重復(fù)步驟②～④對(duì)其他網(wǎng)格進(jìn)行重新初始化，直到所有八叉樹網(wǎng)格單元距離場(chǎng)函數(shù)值全部重新初始化為止。

6 應(yīng) 用

以基于八叉樹網(wǎng)格技術(shù)的Level Set方法模擬不同網(wǎng)格細(xì)化程度的Zalesak圓球在旋轉(zhuǎn)場(chǎng)中的運(yùn)動(dòng)情況以及圓球在剪切場(chǎng)中的運(yùn)動(dòng)情況為例，討論基于八叉樹網(wǎng)格技術(shù)的Level Set方法在三維空間的應(yīng)用，并將該方法應(yīng)用于追蹤方腔鑄件充型過(guò)程的界面前沿，模擬充型過(guò)程中液體的流動(dòng)。

6.1 旋轉(zhuǎn)流場(chǎng)

本算例模擬的是一個(gè)半徑為0.25的Zalesak圓球界面在旋轉(zhuǎn)流場(chǎng)作用下的變化情況，初始位置為[0.5,0.5,0.5]，切口寬度為0.125，切口深度為0.3125，計(jì)算區(qū)域?yàn)閇0,1]×[0,1]×[0,1]，網(wǎng)格尺寸為128×128×128。網(wǎng)格的最大層數(shù)分別為8層和9層，最小網(wǎng)格層數(shù)為4層。旋轉(zhuǎn)流場(chǎng)的速度場(chǎng)為

流場(chǎng)周期為2p，時(shí)間步長(zhǎng)為0.05。

圖7為Zalesak圓球在三維旋轉(zhuǎn)流場(chǎng)中八叉樹自適應(yīng)網(wǎng)格示意圖,圖8和圖9分別是網(wǎng)格最大層數(shù)為8層和9層時(shí)Zalesak圓球在旋轉(zhuǎn)流場(chǎng)中的運(yùn)動(dòng)時(shí)間分別為0、/4、/2、3/4和時(shí)刻的數(shù)值模擬結(jié)果。圖10為劉長(zhǎng)春[22]采用Youngs-VOF方法模擬缺口圓在旋轉(zhuǎn)流場(chǎng)中的運(yùn)動(dòng)的結(jié)果，圖11為張明遠(yuǎn)[23]采用均勻網(wǎng)格的Level Set方法追蹤Zalesak圓球運(yùn)動(dòng)界面的結(jié)果。

通過(guò)圖8和圖9的數(shù)值模擬結(jié)果可以得出：Zalesak圓球經(jīng)過(guò)一個(gè)計(jì)算周期后回到了初始位置，八叉樹網(wǎng)格最大層數(shù)為9層時(shí)Zalesak圓球表面的精度遠(yuǎn)高于網(wǎng)格最大層數(shù)為8層時(shí)Zalesak圓球表面的精度，而且在Zalesak圓球缺口處仍然可以保持尖銳。比較圖9、圖10和圖11的模擬結(jié)果可以發(fā)現(xiàn)：圖10 Youngs-VOF方法中紅色圓圈中的尖銳部分有一些失真，本研究提出的方法在尖銳界面處的模擬精度略好一些；圖11中綠色圓圈中界面部分質(zhì)量損失嚴(yán)重且界面精度也較低，其中的平面部分已經(jīng)平滑成為了曲面，本研究圖9的模擬結(jié)果則很好地保持了原來(lái)的界面形狀，提高了界面追蹤的精度和質(zhì)量守恒性。

6.2 剪切流場(chǎng)

本算例模擬的是一個(gè)半徑為0.125的圓球在剪切流場(chǎng)作用下的變化情況，初始位置為[0.5,0.5,0.25]，計(jì)算區(qū)域?yàn)閇0,1]×[0,1]×[0,1]，網(wǎng)格尺寸為128×128×128。網(wǎng)格的最大層數(shù)分別為8層和9層，最小層數(shù)為4層。剪切流場(chǎng)的速度場(chǎng)為

流場(chǎng)周期為2p，時(shí)間步長(zhǎng)為0.05。

圖12為圓球在三維剪切流場(chǎng)中八叉樹自適應(yīng)網(wǎng)格模型,圖13和圖14分別是網(wǎng)格最大層數(shù)為8層和9層時(shí)運(yùn)動(dòng)時(shí)間分別為0、/4、/2、3/4和時(shí)刻的數(shù)值模擬結(jié)果。圖15為劉長(zhǎng)春[22]采用Youngs-VOF方法模擬圓在二維剪切流場(chǎng)中的運(yùn)動(dòng)的結(jié)果，圖16為張明遠(yuǎn)[23]采用均勻網(wǎng)格的Level Set方法追蹤圓球在三維剪切流場(chǎng)中的運(yùn)動(dòng)界面的結(jié)果。

通過(guò)圖13和圖14的數(shù)值模擬結(jié)果可以得到：圓球在剪切流場(chǎng)中經(jīng)過(guò)一個(gè)時(shí)間周期后又恢復(fù)到原來(lái)的形狀，但是由于Level Set方法的質(zhì)量守恒性差，圓球發(fā)生一些變形。比較圖13和圖14，增加網(wǎng)格最大層數(shù)可以改善Level Set方法在剪切流場(chǎng)中的數(shù)值模擬結(jié)果，提高Level Set方法的質(zhì)量守恒性以及界面的精度。比較圖14、圖15和圖16可以發(fā)現(xiàn)，圖15中的界面追蹤結(jié)果出現(xiàn)了破碎現(xiàn)象，其結(jié)果很不理想，圖16中的均勻網(wǎng)格Level Set方法模擬結(jié)果的質(zhì)量守恒性差，本研究提出的方法則可以較好地保持質(zhì)量守恒性且界面形狀保持較好。

6.3 方腔鑄件充型過(guò)程數(shù)值模擬算例

為了進(jìn)一步證明基于八叉樹自適應(yīng)網(wǎng)格技術(shù)的Level Set方法追蹤界面運(yùn)動(dòng)的可行性，本算例模擬方腔鑄件充型流動(dòng)過(guò)程中金屬液體界面的運(yùn)動(dòng)過(guò)程，該過(guò)程中金屬液體通過(guò)澆注系統(tǒng)流入方腔鑄件并進(jìn)行充型，其計(jì)算區(qū)域?yàn)閇0,1]×[0,1]×[0,1]，網(wǎng)格尺寸為32×32×32。

數(shù)值模擬結(jié)果如圖17所示，圖中紅色區(qū)域表示鑄型，藍(lán)色區(qū)域表示充型金屬液體，紫色區(qū)域表示鑄型中的空氣即待充型區(qū)域。由圖17可以看出采用八叉樹技術(shù)的Level Set方法可以很好地模擬鑄件充型過(guò)程中液體界面的運(yùn)動(dòng)位置，能夠較為準(zhǔn)確地描述鑄造過(guò)程中金屬液體的充型流動(dòng)。

7 結(jié) 論

研究了基于八叉樹網(wǎng)格技術(shù)的Level Set方法追蹤運(yùn)動(dòng)界面，引入八叉樹技術(shù)可以改善Level Set方法模擬運(yùn)動(dòng)界面的一些不足，主要的研究工作總結(jié)如下。

（1）提出了適用于Level Set方法追蹤界面的八叉樹網(wǎng)格自適應(yīng)模型（八叉樹網(wǎng)格技術(shù)的細(xì)化與粗化模型），通過(guò)對(duì)局部運(yùn)動(dòng)網(wǎng)格的細(xì)化與粗化可以減少數(shù)值模擬的計(jì)算網(wǎng)格與內(nèi)存，提高計(jì)算的效率與精度。

（2）建立了一套基于八叉樹自適應(yīng)網(wǎng)格的Level Set方法追蹤運(yùn)動(dòng)界面求解方法，并提出了基于八叉樹自適應(yīng)網(wǎng)格Level Set函數(shù)值離散模型以及求解八叉樹相鄰網(wǎng)格距離場(chǎng)函數(shù)值的插值方法。

（3）提出了八叉樹網(wǎng)格距離場(chǎng)重新初始化方法，通過(guò)快速掃描八叉樹網(wǎng)格并判斷網(wǎng)格是否為界面網(wǎng)格更新非界面網(wǎng)格單元的Level Set函數(shù)值，重新初始化八叉樹網(wǎng)格中的距離場(chǎng)函數(shù)值。

（4）兩個(gè)運(yùn)動(dòng)界面算例模擬結(jié)果表明：將八叉樹技術(shù)引入Level Set方法可以成功地追蹤運(yùn)動(dòng)界面，局部細(xì)化八叉樹網(wǎng)格以及增加網(wǎng)格的最大層數(shù)可以提高Level Set方法追蹤界面的精度以及質(zhì)量守恒性，減少Level Set方法產(chǎn)生的數(shù)值耗散。并用該方法模擬鑄件充型過(guò)程的金屬液體流動(dòng)，結(jié)果表明該方法可以較為準(zhǔn)確地模擬鑄件充型過(guò)程。

符 號(hào) 說(shuō) 明

S(j0)——距離函數(shù)符號(hào) T——運(yùn)動(dòng)周期，s t——時(shí)間，s Dt——時(shí)間步長(zhǎng)，s u, v,w——分別為直角坐標(biāo)3個(gè)方向的速度，s-1 x,y,z——三維直角坐標(biāo)系的3個(gè)方向 Dx,Dy,Dz——分別為3個(gè)方向網(wǎng)格單元長(zhǎng)度 j——Level Set 函數(shù)

References

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

[2] OSHER S, FEDKIW R P. Level set methods: an overview and some recent results[J]. Journal of Computational Physics, 2001, 169(2): 463-502.

[3] SUSSMAN M, SMEREKA P, OSHER S. A level set approach for computing solutions to incompressible two-phase flow[J]. Journal of Computational Physics, 1994, 114(1): 146-159.

[4] WU K T, HAO L, WANG C,. Level set interface treatment and its application in Euler method[J]. Science China Physics Mechanics & Astronomy, 2010, 53(2):227-236.

[5] CHEN S, MERRIMAN B, OSHER S,. A simple level set method for solving Stefan problems[J]. Journal of Computational Physics, 1997, 135(1): 8-29.

[6] SETHIAN J A. A fast marching level set method for monotonically advancing fronts[J]. Proceedings of the National Academy of Sciences, 1996, 93(4): 1591-1595.

[7] SETHIAN J A. Evolution, implementation, and application of level set and fast marching methods for advancing fronts[J]. Journal of Computational Physics, 2001, 169(2): 503-555.

[8] ADALSTEINSSON D, SETHIAN J A. A fast level set method for propagating interfaces[J]. Journal of Computational Physics, 1995, 118(2): 269-277.

[9] PENG D, MERRIMAN B, OSHER S,. A PDE-based fast local level set method[J]. Journal of Computational Physics, 1999, 155(2): 410-438.

[10] STRAIN J. Tree methods for moving interfaces[J]. Journal of Computational Physics, 1999, 151(2): 616-648.

[11] SOCHNIKOV V, EFRIMA S. Level set calculations of the evolution of boundaries on a dynamically adaptive grid[J]. International Journal for Numerical Methods in Engineering, 2003, 56(13): 1913-1929.

[12] COSTA J C A, ALVES M K. Layout optimization with h-adaptivity of structures[J]. International Journal for Numerical Methods in Engineering, 2003, 58(1): 83-102.

[13] 宮翔飛, 張樹道, 江松. 界面捕捉Level Set方法的（AMR）數(shù)值模擬[J].計(jì)算物理, 2006, 23(4): 391-395.GONG X F, ZHANG S D, JIANG S. Numerical simulation of fluid interfaces with the adaptive mesh refinement method, the ghost fluid method and the Level Set method[J]. Chinese Journal of Computational Physics, 2006, 23(4): 391-395.

[14] 黃筱云, 李紹武.空間自適應(yīng)的快速粒子Level Set方法[J]. 天津大學(xué)學(xué)報(bào), 2010, 43(11): 981-987. HUANG X Y, LI S W. Spatially adaptive fast particle Level Set method[J]. Journal of Tianjin University, 2010, 43(11): 981-987.

[15] 黃筱云, 李紹武. 自適應(yīng)四叉樹網(wǎng)格下的N-S方程數(shù)值求解模型[J]. 天津大學(xué)學(xué)報(bào)(自然科學(xué)與工程技術(shù)版), 2013, 46(1): 58-66. HUANG X Y, LI S W. Numerical N-S equation solver based on adaptive quadtree mesh[J]. Journal of Tianjin University(Science and Technology), 2013, 46(1):58-66.

[16] 王生輝, 尤偉, 李增耀. 基于四叉樹直角坐標(biāo)網(wǎng)格的運(yùn)動(dòng)界面追蹤方法[J]. 化工學(xué)報(bào), 2014, 65(S1): 44-50. WANG S H, YOU W, LI Z Y. Interface-capturing method on quadtree-based adaptive Cartesian grid[J]. CIESC Journal, 2014, 65(S1): 44-50.

[17] LOSASSO F, GIBOU F, FEDKIW R. Simulating water and smoke with an octree data structure[J]. ACM Transactions on Graphics, 2004, 23(3):457-462.

[18] 朱雷, 周建新, 龐盛永, 等. 基于八叉樹網(wǎng)格的激光焊接溫度場(chǎng)數(shù)值模擬方法[J]. 熱加工工藝, 2014, 43(5): 164-167. ZHU L, ZHOU J X, PANG S Y,. Numerical simulation method of laser welding temperature field based on octree grid[J]. Hot Working Technology, 2014, 43(5): 164-167.

[19] 殷亞軍. 基于八叉樹網(wǎng)格技術(shù)的相場(chǎng)法金屬凝固過(guò)程組織模擬的研究[D]. 武漢: 華中科技大學(xué), 2013. YIN Y J. Research on the microstructure simulation of the alloy solidification based on phase field method by the octree mesh technology[D]. Wuhan: Huazhong University of Science and Technology, 2013.

[20] 劉元元. WENO方法的研究及其在對(duì)流擴(kuò)散方程中的應(yīng)用[D]. 合肥: 中國(guó)科學(xué)技術(shù)大學(xué), 2012. LIU Y Y. WENO methods and the applications in convection diffusion equations[D]. Hefei: University of Science and Technology of China, 2001.

[21] 鄭素佩, 歐陽(yáng)潔, 張紅平, 等. 帶有矩形嵌件薄壁型腔內(nèi)熔接過(guò)程動(dòng)態(tài)模擬[J]. 化工學(xué)報(bào), 2008, 59(1): 232-238. ZHENG S P, OUYANG J, ZHANG H P,. Dynamic simulation for weld lines in thin mold with rectangle cylinder[J]. Journal of Chemical Industry and Engineering (China), 2008, 59(1): 232-238.

[22] 劉長(zhǎng)春. Level Set方法模擬界面運(yùn)動(dòng)[D]. 北京: 中國(guó)科學(xué)院工程熱物理研究所, 2007. LIU C C. Level Set methods for moving interfaces[D].Beijing: Institute of Engineering Thermophysics, Chinese Academy of Sciences, 2007.

[23] 張明遠(yuǎn). 基于Projection方法的鑄造充型過(guò)程數(shù)值模擬[D].武漢: 華中科技大學(xué), 2011. ZHANG M Y. Numerical simulations of casting’s filling process based on a Projection method[D]. Wuhan: Huazhong University of Science and Technology, 2011.

Interface capturing method based on Level Set with octree adaptive mesh technology

YIN Yajun, LI Yangdong, TU Zhixin, LI Wen, SHEN Xu, ZHOU Jianxin

(State Key Laboratory of Materials Processing and Die & Mould Technology, Huazhong University of Science and Technology, Wuhan 430074, Hubei, China)

Level Set method is one of the dominant methods to calculate and track evolution of interfaces because it can efficiently handle complex topology change of interfaces and large deformation. This method often suffers from large amount of numerical dissipation and high order accurate numerical discretization in time and space is used to improve mass conservation properties, which leads to an increase of computational overhead. When octree adaptive mesh technology is incorporated into Level Set to capture moving interfaces, the amount of computational meshes and memory usage is reduced and the calculation precision and efficiency are enhanced by adaptive octree mesh refining or coarsening according to interface distance. Because octree mesh is non-uniform and layer numbers of neighbor grids may be different, the WENO formula could not be directly applied to set values of octree grids. A discretization model and a technique of re-initializing distance field were proposed to calculate value of neighbor octree grids as well as to reduce mass loss in Level Set function. Numerical simulation on moving interface in rotation flow field of Zalesak sphere, in shear flow field of sphere and in cast filling process of square cavity benchmark indicated that combination of Level Set method and octree adaptive mesh technology could successfully capture moving interface. The interface-tracking precision was enhanced by refinement of octree meshes and increase of maximum octree mesh layers. When the maximum octree mesh layer was increased from 8 to 9 in the shear flow field, mass conservation of the sphere was obviously improved and numerical dissipation caused by Level Set method was significantly reduced.

interface capture; numerical simulation; octree mesh technology; Level Set method; adaptive mesh; model; optimization

2016-04-05.

Prof. ZHOU Jianxin,zhoujianxin@hust.edu.cn

10.11949/j.issn.0438-1157.20160421

O 351.2

A

0438—1157（2016）11—4732—10

殷亞軍（1985—），男，博士，講師。

國(guó)家自然科學(xué)基金項(xiàng)目（51475181，51305149）；中國(guó)博士后科學(xué)基金特別資助項(xiàng)目（2015T80795）。

2016-04-05收到初稿，2016-08-03收到修改稿。

聯(lián)系人：周建新。

supported. by the National Natural Science Foundation of China(51475181, 51305149) and the China Postdoctoral Science Foundation Special Funding(2015T80795).