COMSOL Multiphysics? 軟件提供了四種可用于模擬自由液面的方法：水平集、相場、動網格和穩態自由表面。作為系列博客的第一節，我們將討論水平集和相場法，這兩種基于場的方法幾乎可以描述任何類型的自由液面。在第二節中，我們計劃將本文的求解結果與通過動網格接口獲得的結果進行比較。

何為水平集和相場法？

水平集和相場法都是基于場的方法，這類方法將自由液面表征為水平集或相場函數的等值面。自由液面對應固定網格框架下的液體和氣體之間的相分界面。

下圖為管道內兩顆液滴的表面，摘自附加產品“微流體模塊”的“案例庫”所提供的液滴破碎模型。從這張圖中可以看出，盡管液滴的表面非常明顯，但液滴周圍的單元并沒有貼合到液滴表面上。

不管采用水平集方法還是相場法，液滴表面與單元表面都不貼合。

水平集和相場函數都是由納維-斯托克斯方程計算的速度矢量進行對流傳輸的。在水平集和相場法中，對應公式為：

(1)

需要注意的是，水平集和相場函數都使用了 Φ。二者的不同在于方程右側的 F。在初始水平集方法中 F = 0，因此得到純對流傳輸方程。然而當 F = 0 時，數值解不僅不穩定，而且大部分情況下實用性很小。所以為了保持相界面清晰，我們在水平集方法的 F 中添加了高階導數項 Φ 。

在相場法中，F 代表設法將系統的自由能最小化的項。此項也引入了高階導數 Φ。實際上，相場方程中的源項中包含了四階項。這意味著，出于實用性考慮，方程經常被分解為兩個方程，與此同時，輔助因變量被定義為 Φ 的二階導數函數形式。COMSOL Multiphysics 中也采取了這種做法。

兩種方法均將自由液面的表面張力引入到納維-斯托克斯方程的源項中。水平集方法利用表征自由邊界的水平集等值面的曲率來描述表面張力。相場法根據化學勢計算出表面張力對納維-斯托克斯方程的源項貢獻，因為正是化學勢導致了界面附近產生的表面張力和相場函數梯度的產生。

在給定的水平集或相場函數值范圍內，即自由表面的數值范圍內，流體特性從液體平滑地過渡到氣體。

穿過自由表面，水平集函數 Φ 在 0 和 1 之間變化，在兩種流體內部分別為常數 0 或 1。比如在液相中為 0，氣相中為 1。至于自由液體面處，也就是液體和氣體之間的分界面上，對應的水平集函數值為 Φ = 0.5。根據下方公式，密度是水平集函數的函數：

(2)

動力粘度也是如此：

(3)

從這兩個方程中我們可以看出，當 Φ = 0 時，可得到流體 1（比如液體）的屬性，當 Φ = 1 時，可得到流體 2（比如氣體）的屬性。由此，流體屬性在整個自由表面上實現了平滑過渡，像水平集函數一樣平滑。

相場函數 Φ 在 -1 和 1 之間變化，當 Φ = 0 時的等值面為自由液面。粘度和密度的計算是通過與水平集法相似的方式完成的，只不過因為相場函數在 -1 和 1 之間變化，所以表達式不同。

水平集與相場接口的設置

水平集 和相場 接口中包含若干個參數，為了使計算取得最優效果，我們需要對這些參數進行調整。參數的值取決于被模擬的系統和方程的數值離散。幾何結構的尺寸（長寬高）、流體的特性、壁的潤濕條件、初始條件、工作條件以及網格的單元尺寸都會對模型設置選項產生影響。

水平集接口

我們從水平集方法的設置開始。用戶界面中的重新初始化參數 γ 可以確保水平集函數中的梯度隨時間推移逐漸集中到自由表面上，也就是下方示例中水和空氣之間的交界面。它沒有設定表面厚度，但能夠確保水平集函數的變化在厚度范圍之內。如果該參數值設置過低，水平集函數的變化可能被截留在一種流體域內，然后憑空生成液-氣界面。過高的值往往造成時間步過小、計算時間過長。

水平集接口的設置窗口。

ε 表示的界面厚度參數，它理解起來更加直觀，唯一的作用是控制水平集函數發生變化的區域的厚度。過厚會模糊自由表面的形狀。雖然這讓問題變得更容易求解，但是損失了自由表面的形狀精度。較小的 ε 值有利于生成清晰的表面，但這也要求網格具有相應的分辨率。設定小于網格單元尺寸的值是沒有意義的，因為這會使得界面變得不規則，而無法解析的表面張力可能會導致速度失真。ε 的設定值應與表面附近的最大單元大約相同。

相場接口

接下來是相場法，界面厚度參數 ε 與水平集方法基本一致，前文中的論證同樣適用于此設定值。

相場接口的設置。

相場法的遷移調整參數 χ 在某種程度上決定了相場的擴散系數。它的值必須足夠大，才能保證相場方程穩定，但是為了保證表面尖銳，同時也要足夠小。合理值與表面的速度成正比，與表面張力系數成反比。

(4)

這意味著一旦確定了某組工作條件的 χ 值，就可以利用上述關系為一組新的條件設定 χ 的值。

自由液面建模

通過下圖中的示例問題，我們將演示如何在 COMSOL Multiphysics 中利用水平集和相場法進行自由液面建模。

實心條一半浸沒在小型管道內的水中。沿與水面相切的方向來回移動實心條，使水面產生柔和的波浪。為了保持層流狀態，管道和實心條的尺寸必須足夠小。

示例問題的幾何結構和定義。

在此例中，我們利用動網格功能指定小方形條在水面上反復運動。不過，在水平集和相場法中，網格不會隨液體表面移動。

下面是關于在 COMSOL Multiphysics 中創建模型的幾點備注：

• 添加重力作為動量方程的源

• 添加參考壓力和壓力約束點

• 為了方便比較結果，對壁應用 Navier 滑移條件，并使滑移長度等于單元長度

水平集與相場法的仿真結果

下方的結果圖顯示了 0.07 秒、0.57 秒和 1.0秒的流場和水面形狀。兩種方法在各個時間點上的流線和水面形狀均大致相同。最大速度和水面高度稍有差異，原因可能在于兩種方法處理表面張力的不同方式。

再經過一段時間，流線也出現了差別，例如 t = 1.0 s 秒時的回流區。通過相場法得到的液面相對更加平靜。水平集方法利用基于水平集函數的梯度獲得的表面曲率來計算表面張力。所以，與相場法中更加平滑的力相比，水平集方法得到的表面力更“尖銳”。

從上至下依次為：0.07 秒、0.57 秒和 1.0 秒后水平集（左）和相場（右）方法得到的結果。

水平集和相場法還能計算自由液面上的空氣域流場。從圖片中可以看出，長方形條的運動在整個相界面上方引起了持續顯著變化的流場。

利用相場法計算 0.57 秒后水域和空氣域內的流場。

如果增強攪動能使液面運動更劇烈，那么液面可能先破裂再合并，如下方動畫所示。這也是水平集與相場法的一個優勢：我們能夠更簡單地處理自由表面的拓撲變化。

方形條的攪動頻率和幅度增大，致使較小的波浪破裂，并生成了滯留在水相中的氣泡。

盡管水平集和相場法類似，但表面張力的處理對二者的穩定性產生很大的影響，至少在 COMSOL Multiphysics 中存在區別。在處理對表面張力非常敏感的問題時，相場法在求解時間方面比水平集方法表現得更好，原因在于利用水平集方法計算表面曲率時，瞬態求解器需要采用比相場法小得多的時間步長。

在此例中，水平集方法的平均時間步長是相場法的六倍，所以水平集方法需要六倍的計算時間。因此，對于較小尺度自由液面問題和對表面張力極其敏感的層流（例如微流體）問題，相場法通常是更好的選擇。

自適應網格細化與自由表面

針對本文探討的二維示例，在采用了水平集和相場法的整個自由表面區域中，網格足夠密集。不過，對于三維模型而言，這種水平的分辨率的計算成本太高。一種替代方案是使用自適應網格細化功能，它能夠自動根據我們選擇的任何函數創建更密集的網格。

舉例來說，下方動畫展示了基于最大速度梯度的位置（上圖）和相場函數的最大梯度（下圖）生成的網格自適應結果。繪圖展示了氣相和水相。需要注意的是，與剪切率相關的自適應網格在氣相中生成了細化網格，而水相的細化程度相對較低。在此例中，我們對水相更感興趣，因此需要通過縮放相場函數來修改自適應網格的誤差指示器。

利用速度梯度（上）和相場梯度（下）作為誤差指示器的自適應網格細化效果。

動網格與自由表面

本文討論了兩種基于場的自由表面建模方法：相場法和水平集方法。在此系列的第二篇文章中，我們將使用動網格為自由表面建模，并比較兩篇文章介紹的方法。敬請關注！

來源：COMSOL