第九講：

VOF(Volume-Of-Fluid)

方法及其應用

兩相流數值模擬方法的分類分類方法很多氣(液)－液兩相流的相界面描述方法可分為兩大體系：即（１）界面捕捉法

(Front-CapturingMethods)（２）界面跟蹤法

(Front-TrackingMethods)。汽－液兩相流界面描述方法的分類

界面追蹤方法，直觀的理解，屬于Lagarange類方法。如果有足夠多雙眼睛盯住／“鎖定”界面上的各個點，通過對各個點的觀察、記錄就可以確定整個界面的運動過程和軌跡。這就是界面追蹤方法產生的最初的出發點和本質。

這類方法中最具代表性的是：

PIC(Particle-in-Cell)

MAC(Mark-and-Cell)界面追蹤法（FrontTrackingMethods）界面捕捉法中最具有代表性的是：

１LevelSet２VOF３高度函數法

界面捕捉法（FrontCapturingMethods）１．1

VOF方法的基本思想：

1981年，C.W.Hirt

和B.D.Nichols在InternationalJournalofComputationalPhysics雜志上首先正式發表了著名的VOF論文，開創性地提出了用VOF方法進行運動相界面追蹤的思想，并首先用VOF方法對潰壩和浪涌（BrokenDam,Breakingbore）及RayleighTaylor不穩定性現象進行了成功的數值模擬模擬。一、VOF方法的基本思想VOF方法的基本思想VOF方法在整個流場定義了一個相函數f，在每個網格中，這個函數的定義為一種流體（目標流體）的體積與網格體積的比值。相函數滿足一定的輸運方程，在得出相函數輸運方程之前先定義一個染色函數從而在控制體σ上，目標流體的體積函數表達式為

由拉格朗日流體體積的特性可知

將上述全導數展開可得輸運方程為相函數f的定義式為因此可得相函數滿足的輸運方程為當f=0時，網格內沒有目標流體，稱為空網格；當f=1時，網格被目標流體充滿，稱為滿網格；當0<f<1時，網格含有目標流體，但未被充滿，稱之為界面網格。（１）VOF方法用相函數(PhaseFunction)F取代了MAC方法中的虛擬無質量彩色粒子。從這個意義上說，VOF方法可看作是MAC方法的一個變種。（２）在一種流體相（比如說，液體）中，相函數F取值為1，而在另一種流體（比如，氣體或另一種液體）中取值為0；在相函數取0到1之間的數值的地方即為相界面位置。（３）相界面的取向可由界面附近各點上的F值來確定。一、VOF方法的基本思想

（4）VOF方法避開了采用工程浩大的Marker點的方案，比MAC方法對計算機內存的要求低，更容易實施。（5）Hirt和Nichols在最初的VOF模型中，設計了類似于MAC方法和PIC方法的Staggered型差分格式——即將壓力和相函數（體積分數）定義在格子的中心處，而將X方向的速度定義在格子的左、右格邊中點，將Y方向的速度定義在格子的上、下格邊中點。一、VOF方法的基本思想１．２關于相函數的概念：（１）相函數類似于氣－液兩相流中的截面含氣率（容積含氣率），表示某一相介質占據網格面積（二維）或體積（三維）的分數。（２）相函數是以一個網格為單元來定義的，與含氣率的概念不同。一、VOF方法的基本思想１．２關于相函數的概念：（３）相函數就是一個介質指針，指示著占據某一望個的介質種類。但這種指示不僅是定性的，而且是定量的。（４）對應于同一個相函數值，氣－液界面在網格內的形狀和方位是多值的。如圖所示。（５）F必須而且只能在0和1之間。一、VOF方法的基本思想VOF方法中的基本問題：1.

氣體的動量／運動控制方程2.

液體的動量／運動控制方程3.

相函數F的控制方程及其求解4.

如何由相函數F的分布獲得氣－液相界面——相界面的構造問題。二、VOF方法中的基本問題：三、氣－液兩相流的動量方程

——兩相流運動控制方程三、氣－液兩相流的動量方程

——兩相流運動控制方程三、氣－液兩相流的動量方程

——兩相流運動控制方程三、氣－液兩相流的動量方程

——兩相流運動控制方程四．氣－液相界面的控制方程四．氣－液相界面的控制方程五．氣－液相界面控制方程的求解方法五．氣－液相界面控制方程的求解方法五．氣－液相界面控制方程的求解方法與LevelSet方法不同的是：在LevelSet方法中，求得了LevelSet函數，由LevelSet函數的０等值面可立即得到相界面的位置和形狀；但是，在VOF方法中，得到了相函數F的分布之后，還有一個必須解決的關鍵問題？——即如何根據相函數F的分布準確地確定相界面在每一個時間層上的空間位置，也即如何實現所謂的“相界面重新構造”。六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法方程（４）中的最小函數“Min｛s｝”主要是為了防止從“施主”網格中過多的流出流體;——在一個時間單元內，流過（流出／流入）一個網格的距離不可能超過單元格的寬度！方程（５）中的最大函數“Max{s}”主要是為了考慮在計算中可能出現的、由于非目標流體的流動量超過了可能獲得的非目標流體量而導致的額外附加的流體流動。六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法

施主-受主的思想：Hirt和Nichols在提出VOF方法的同時，給出了相函數的求解方法，兩位學者設計了類似于MAC方法和PIC方法的Staggered型差分格式，即相函數定義在網格中心處，而速度定義在網格邊界上。在求解的過程中采用了施主-受主的思想，具體實施過程如下：

采用施主-受主型的差分對對流項進行逼近得：其中六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法CFL條件，以Courant,Friedrichs,Lewy三個人的名字命名界面重構的實施范圍六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法則單位法向為六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法6.3FLAIR型重構技術FLAIR型重構技術由Ashgriz和J.Y.Poo兩位教授首先于1990年提出。FLAIR一詞的英文原文是aFlux-Line-Segment

ModelforAdvectionandInterfaceReconstruction，即對流和界面重構的通量－線段模型。該方法運用直線段來近似兩個相鄰網格內的界面。基本的思路：由其名稱本身已可以看到其基本的思路：——即對任意與界面相關的兩個相鄰網格，通過構造一條帶有傾角的直線段作為跨過該網格邊界的近似界面。然后計算單位時間內流過該網格邊界的流體體積通量(用flux表示)，并作為修改流體體積函數的數值流通量。六．氣－液相界面的構造方法由于界面的構造要涉及到相鄰的兩個網格，因此需要分成多種情況進行討論計算。————根據相鄰網格的相含率的不同及施主與受主網格的判斷，即按通過網格邊界的流體速度確定施主單元和受主單元，分別用FD和FA表示，根據FD和FA的值可以分為以下五種情況。

值得首先說明的是，進行這種分類的目的，主要是為了計算通過網格邊界的流體體積通量。這種通量主要是由施主網格的狀態決定的，因此這里的分類也主要以施主網格的狀態為依據。這種分類主要為了厘清計算思路、方便計算過程的實施。六．氣－液相界面的構造方法六．氣－液相界面的構造方法（1）施主網格是滿網格，受主網格為滿網格或空網格，即fa=1.0，fb=0.0/1.0（2）施主網格是空網格，受主網格為滿網格或空網格，即fa=0.0，fb=0.0/1.0六．氣－液相界面的構造方法（1）施主網格是滿網格，受主網格為滿網格或空網格，即fa=1.0，fb=0.0/1.0（2）施主網格是空網格，受主網格為滿網格或空網格，即fa=0.0，fb=0.0/1.0六．氣－液相界面的構造方法六．氣－液相界面的構造方法針對施主和受主網格均為半網格時的上述四種情況，其判斷依據為：（1）若為第一種情況應滿足：（2）若為第二種情況應滿足：（3）若為第三種情況應滿足：（4）若為第四種情況應滿足：六．氣－液相界面的構造方法以下將簡單介紹第一種情況的輸運：（1）確定界面的斜率界面方程為：通過比較相鄰網格內的相含率fa、fb可得界面表達式的系數為：（2）確定界面的流量根據相鄰網格相鄰界面的速度以及界面方程可得界面的輸運為:式中：f,Fs,f+,f-六．氣－液相界面的構造方法六．氣－液相界面的構造方法上述四種情況的判斷依據為：（1）若為第一種情況應滿足：（2）若為第二種情況應滿足：

（3）若為第三種情況應滿足：（4）若為第四種情況應滿足：直線斜率?上述4種界面的構造最為復雜，其主要體現在界面斜率的計算。網格（i,j）界面的計算要在垂直方向上運用兩次類型3的界面構造技術，首先對網格（i,j）、（i,j-1）進行界面構造，得到斜率β1，然后再對網格（i,j）、（i,j+1）進行計算，得到斜率β2，最后取兩斜率的均值作為網格（i,j）界面的斜率β。

然后，再根據界面流速計算界面流量。

六．氣－液相界面的構造方法六．氣－液相界面的構造方法再進一步，可按類似于情形（３）中的方法構造直線，作為此網格內的界面近似，然后求出經過網格邊界流入受主網格的流體體積流量，同樣可以歸結為圖４所示的４種情形。六．氣－液相界面的構造方法

FLAIR算法的算例——剪切速度場取(x0,y0)=(0.5,0.5)，計算區域取[0,1]*[0,1]，初始界面為圓，圓心為(0.5,0.3)、半徑為0.2。計算網格為200*200時間步長為0.0001s。

圖3.3剪切速度場數值結果（上為旋轉1s，下為旋轉2s）六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法PLIC界面重構技術在二維直角坐標中的實現PLIC界面重構技術也是應用直線段對界面進行重構，但其與FLAIR不同，PLIC方法是在單個網格內對界面進行近似。PLIC界面存在20種情況，其中4種情況較為簡單，界面與坐標軸平行，其界面流量的計算非常簡單。其它16種情況的界面與兩個坐標軸均有交點，此界面的法向為：則單位法向為分母上系數，二維時為8；三維時為32？？通過坐標軸的旋轉以及坐標轉換可將各種界面轉換為的情況，的界面存在4種情況，分別為：圖3.4PLIC最終簡化類型上述4種情況是以圖3.5中的A、B兩點為分界的。圖3.5中界面的方程為：圖3.5PLIC界面示意圖當時，只存在圖3.4中的1、2、4情況，此時界面沿法線方向移動會先后經過A、B點。界面經過A點時，陰影面積為：界面經過B點時，陰影面積為：方程？若，則界面為情況1；若，則界面為情況2；若，則界面為情況4。當時，只存在圖3.4中的1、3、4情況，此時界面沿法線方向移動會先后經過B、A點。界面經過B點時，陰影面積為：界面經過A點時，陰影面積為：若，則界面為情況1；若，則界面為情況2；若，則界面為情況4。C=?相函數F，fC=?相函數F，f圖3.6中陰影面積為：式中，M、N為系數，若界面為情況1，M=1、N=1；若界面為情況2，M=1、N=0；若界面為情況3，M=0、N=1；若界面為情況4，M=0、N=0。

圖3.6PLIC界面示意圖通過上式可以求得各情況時，界面與x軸的交點為：第一種情況：第二種情況：第三種情況：第四種情況：

現以情況4為例介紹目標網格向周圍4個網格的流量輸運。圖3.6PLIC界面的第4種情況(1)網格i,j向i,j+1網格的流量輸運：如果Vi,j+1<0，不存在輸運；如果Vi,j+1>0，則存在輸運：

IFVi,j+1*dt<=dy-y1

else

（2）網格i,j向i,j-1網格的流量輸運：如果Vi,j>0，不存在輸運；如果Vi,j<0，則存在輸運：

IF|Vi,j|*dt>y1else（3）網格i,j向i+1,j網格的流量輸運：如果Ui+1,j<0，不存在輸運；如果Ui+1,j>0，則存在輸運：

IFUi+1,j*dt<dx-x1else（4）網格i,j向i-1,j網格的流量輸運：如果Ui,j>0，不存在輸運；如果Ui,j<0，則存在輸運：

IF|Ui,j|*dt>x1else其它幾種情況的界面輸運與情況4相近，求得各界面的輸運以后就可以得到新時刻的相函數分布。

PLIC算法的驗證算例如下：（1）剪切流場取(x0,y0)=(0.5,0.5)，計算區域取[0,1]*[0,1]，初始界面為圓，圓心為(0.5,0.3)、半徑為0.2。計算網格為200*200時間步長為0.0001s。

圖3.7剪切速度場數值結果（上為旋轉1s，下為旋轉2s）（2）旋轉流場：取(x0,y0)=(0.5,0.5)，計算區域取[0,1]*[0,1]，初始界面為開口圓，圓心為(0.5,0.5)、半徑為0.25、開口寬度為0.08、開口底部距圓心的距離為0.04、開口向下。計算網格為200*200時間步長為0.0001s。圖3.8旋轉速度場數值結果（依次旋轉0.5s、1.0s、1.5s、2s）

PLIC方法在多相流中的應用：（1）氣泡在水中的上升過程計算區域為：XL=0.2m，YL=0.2m；計算所采用的網格為102*102；計算時間步長為0.0001s；初始氣泡的半徑為0.02m；圓心為（0.1,0.1）。氣水兩相的物性分別為：水的密度為1000.0kg/m3、水的粘度為0.001Pa.s；空氣的密度為1.266kg/m3、空氣的粘度為17.8E-6Pa.s；氣水表面張力為0.0728N/m。界面重構方法圖3.9氣泡上升的數值結果（依次0ms、50ms、100ms、130ms、170ms、200ms）圖3.10氣泡上升的流場（依次100ms200ms）（2）PLIC方法在油水兩相流動中的應用r界面錐

油水混合區不滲透固壁油水相界面計算區域(a)油洞抽象示意圖thLd(b)物理模型及尺寸示意圖圖3.11溶洞油藏的物理模型與計算模型六．氣－液相界面的構造方法

3、PLIC界面重構技術在三維直角坐標中的實現

3、PLIC界面重構技術在三維直角坐標中的實現在三維網格中同樣存在滿網格、空網格、相網格，空網格不需要計算界面的流量輸出，而滿網格的界面輸運計算非常簡單，工作的重點為相網格。根據相網格中界面的法向可將界面分為以下三種：（1），即界面與網格的一個面平行，此界面的輸運比較好計算；（2），即界面與網格的一個坐標軸平行，此界面的輸運與二維PLIC的界面輸運相近，因此求解并不是很困難；（3），界面與所在的平面與三個坐標軸均有交點，此界面存在情況較多，求解也較為困難。第一種類型界面的輸運求解非常簡單；第二種類型界面與二維PLIC相近；因此下面將不對這兩種類型的界面進行分析。圖3.16三種界面類型的示意圖

的界面情況非常多，因此需要簡化。首先可通過坐標軸旋轉、坐標轉換及目標流體與非目標流體的轉化將各種情況轉為，此時可將界面分為5種情況，如下：情況1情況2情況3情況4情況5圖3.173D-PLIC的5種界面情況

界面方程為

式中，

上述5種情況的判斷準則為：（1）當界面為第一種情況時，應滿足下式：

（2）當界面為第二種情況時，應滿足下式：（3）當界面為第三種情況時，應滿足下式：（4）當界面為第四種情況時，應滿足下式：（5）當界面為第五種情況時，應滿足下式：以情況1為例來分析目標網格向周圍6個的流量輸運目標流體所占的體積為：因此：

所以可求得界面在x軸上的截距為：圖3.18界面的第一種情況

相網格向周圍6個網格的輸運如下：(1)網格i,j,k通過PQGF面向i+1,j,k網格的流量輸運：如果Ui+1,j,k<0，不存在輸運；如果Ui+1,j,k>0，則存在輸運：

IFUi+1,j,k*dt<=AF

else

(2)網格i,j,k通過MNKE面向i-1,j,k網格的流量輸運：如果Ui,j,k>0，不存在輸運；如果Ui,j,k<0，則存在輸運：

IF|Ui,j,k|*dt>=AE

else

(3)網格i,j,k通過BMKG面向i,j+1,k網格的流量輸運：如果Vi,j+1,k<0，不存在輸運；如果Vi,j+1,k>0，則存在輸運：

IFVi,j+1,k*dt<=BK

else

(4)網格i,j,k通過PNEF面向i,j-1,k網格的流量輸運：如果Vi,j,k>0，不存在輸運；如果Vi,j,k<0，則存在輸運：

IF|Vi,j,k|*dt>=BE

else

(5)網格i,j,k通過PNMQ面向i,j,k+1網格的流量輸運：如果Wi,j,k+1<0，不存在輸運；如果Wi,j,k+1>0，則存在輸運：

IFWi,j,k+1*dt<=CN

else

(6)網格i,j,k通過FGKE面向i,j,k-1網格的流量輸運：如果Wi,j,k>0，不存在輸運；如果Wi,j,k<0，則存在輸運：

IF|Wi,j,k|*dt>=CE

else

其它幾種情況的界面輸運與情況1相近，求得各界面的輸運以后就可以得到新時刻的相函數分布。

PLIC界面重構技術的驗證算例：

（1）常數速度場算例為立方體在上述流場中的變化，三維計算區域為，所計算立方體的下角點為、邊長為0.4，計算網格分別為、，時間步長為0.001s。圖3.19常數速度場的驗證（網格為52*52*52、時間分別為0.0s、0.5s、1.0s）

圖3.20常數速度場的驗證（網格為102*102*102、時間分別為0.0s、0.5s、1.0s）

（2）旋轉速度場考察了Zalesak’s球在流場中的變化，此計算區域為，Zalesak’s球的球心為，半徑為0.15，

計算網格分別為、，計算時間步長為0.001s。

圖3.21旋轉速度場的驗證（網格為102*102*102、時間分別為0.0s、1.0s、2.0s）

圖3.22旋轉速度場的驗證（網格為152*152*152、時間分別為0.0s、1.0s、2.0s）

（3）剪切速度場算例為圓球在上述流場中的變化，三維計算區域為，所計算球體的求心為、半徑為0.15，計算周期為T=3s，計算網格為、，計算時間步長為0.001s。圖3.23剪切速度場的驗證（網格為127*127*127、時間分別為0s、0.5s、1.0s、1.5s、3.0s、2.5s、2.0s

）

圖3.24剪切速度場的驗證（網格為152*152*152、時間分別為0s、0.5s、1.0s、1.5s、3.0s、2.5s、2.0s

）

3D-PLIC技術在氣泡動力學中的應用：（1）單個氣泡的上升過程t=0mst=10mst=20mst=30mst=35mst=40ms圖3.25單個氣泡的上升過程（2）兩個氣泡的融合過程t=0mst=20mst=40mst=60mst=80ms圖3.26兩個氣泡的融合過程六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法六．氣－液相界面的構造方法七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較七．數值模擬算例與比較八．關于相界面的構造問題的幾點說明八．關于相界面的構造問題的幾點說明九．VOF方法中氣、液兩相流場中的物性表示方法十．VOF方法的優點十一．VOF方法的缺點十一．VOF方法的缺點十二．VOF方法的應用VOF方法在氣泡動力學中的應用：圖1.1VOF方法在氣泡動力學中的應用

（M.van

Sint

Annaland2004