一種強耦合預估-校正浸入邊界法*

張和濤，寧建國，許香照，馬天寶

（北京理工大學爆炸科學與技術國家重點實驗室，北京 100081）

計算力學領域中的流固耦合問題具有廣闊的應用背景，長期以來吸引著眾多研究者關注。流固耦合仿真的理論和數值模擬一直是計算力學領域中難度較大的課題之一。浸入邊界法由Peskin[1-2]于1972年首次提出，目前已發展成為流固耦合數值模擬中廣泛采用的一類方法[3]，被應用于研究聲場[4-6]、渦激振動[7-10]、納米流體[11]、湍流[12]等復雜問題。典型的浸入邊界法采用獨立的拉格朗日網格描述固體結構并將之嵌入歐拉網格描述的流場中。通過對邊界附近進行特殊處理，浸入邊界法確保了流固耦合邊界上的無滑移條件和力平衡關系得以滿足。

浸入邊界法的一個分支是銳利浸入邊界法[13]。銳利浸入邊界法嚴格保持邊界位置上的強間斷，根據邊界兩側的應力推斷相互作用能夠始終維持銳利邊界。借助于邊界單元的特殊處理規則，銳利浸入邊界法通常無需在流體Navier-Stokes方程中添加附加源項，因此對格式穩定性影響較小，有利于推廣到可壓縮流體情形。理論上，銳利浸入邊界法可以具有高階精度[14]，在模擬高雷諾數流動時相比浸入邊界法的另一分支模糊浸入邊界法[15-16]具有優勢。銳利浸入邊界法的核心問題是處理被邊界截斷的單元，簡稱截斷單元（cut-cell）。主流實現方式分為兩類：虛擬網格法（ghost-cell method）[17-28]和截斷單元法（cutcell method）[29-33]，其中應用相對較多的是虛擬網格法。對虛擬網格法的研究表明，為了處理截斷單元，虛擬網格法需要進行多維重構，常用方法包括形函數插值和移動最小二乘法等。

虛擬網格法存在的一個問題是當耦合邊界在流體的固定網格上移動時不能維持嚴格的質量守恒[23]。邊界運動時會改變固體的覆蓋域，在流體網格上產生許多失效單元和新增單元。被更新后的固體覆蓋的失效單元會損失其中包含的所有信息并造成失效單元質量丟失。為了獲得正確的壓力，虛擬網格法有時需要對新增單元進行填充，導致新增單元獲得額外質量。因此整體上虛擬網格法不能維持質量守恒。Seo等[34]的研究確認違反幾何守恒律及相應的質量損失是造成虛假振蕩的首要原因，并提出了一種針對不可壓縮流體的改善方案。目前對可壓縮流體仍無合適的解決方法。

截斷單元法考慮了截斷單元體積變化的過程，因此通過合適的單元處理可以實現精確的質量守恒。Monasse等[31]提出了一種針對可壓縮流體和剛體耦合的守恒型浸入邊界法，通過校正Godunov通量來還原損失的質量，但由于采用了有限體積法，涉及截斷單元的幾何計算會比較困難。與Monasse的方法類似，Sotiropoulos等[13]給出的結論表明，在截斷單元法中對截斷單元的拓撲管理是主要問題，妨礙了這一類方法推廣到三維情形。Brady等[33]指出，當截斷單元體積趨于零時會引起系統剛性問題。用有限差分法替代有限體積法可以顯著簡化截斷單元的計算。然而由于有限體積法在計算流體力學領域的廣泛應用，考慮有限體積法的耦合問題仍然是必要的。

為了克服虛擬網格法和截斷單元法存在的困難，本文中，將提出一種區別于這兩類方法的新思路，即適用于可壓縮流體和結構耦合問題的強耦合預估-校正銳利浸入邊界法。首先，分析虛擬網格浸入邊界法導致質量不守恒的原因，揭示網格拓撲變化的影響。其次，闡述一般流固耦合系統的矩陣表示，導出強耦合Gauss-Seidel迭代格式，進一步導出預估-校正格式，基于該格式，提出預估-校正算法，預估步驟使用無耦合邊界的流體模型求解，并對流體施加校正。然后，介紹獨立的介質求解器，并詳細闡述校正步驟，在校正步驟中，通過一組輸運規則來實現輸運過程，保證質量守恒。最后，使用該方法測試若干數值算例，檢驗該方法的準確性和收斂性。

1 虛擬網格法

有限體積法在重構流場的通量函數時需要使用相鄰多個單元組成的模板。當通量模板包含至少一個截斷單元時，一般采取構造虛擬網格的方式補全模板。一種常用的方法是使用基于移動最小二乘法的虛擬網格法補全模板。如圖1所示，下方的4個單元組成一個通量模板。單元G被浸入的固體網格覆蓋。斜實線代表浸入邊界。點G關于浸入邊界的鏡像為點I，兩點連線與浸入邊界相交于點H。距點I一定半徑范圍內的單元組成一個支撐域 ?I。在 ?I上使用徑向基函數移動最小二乘法擬合點I的物理量。

圖1 根據鏡像點I構造虛擬網格結點GFig.1 Construct the ghost pointGbased on its image pointI

根據點I的物理量和浸入邊界上的邊界條件可構造虛擬結點G。為了與固定邊界情形下的邊界條件一致，假設I和G關于H參考系滿足固壁邊界條件,可表達為：

式中： ρ 為密度，Ht為總能，u為速度。該邊界條件的兩種情形分別對應于Dirichlet邊界條件和Neumann邊界條件。

虛擬網格法會導致流體質量的改變。如圖2所示，當浸入邊界在固定的歐拉網格上移動時，會覆蓋一部分網格從而造成失效單元，同時一部分網格不再被覆蓋而成為新增單元。虛擬網格法在構造虛擬網格時沒有考慮到網格拓撲變化的影響。當出現失效單元時，網格拓撲發生突變，與連續的邊界條件產生矛盾。雖然施加了邊界條件，但是仍然會損失一部分質量。當出現新增單元時，為了擬合的光滑性，可能需要使用擬合方法對新增單元進行填充。同樣地，在新增單元上的擬合算法中加入邊界條件約束也不能保證質量守恒。

圖2 從tn 時刻到tn+1時刻邊界進行移動，造成紅色的失效單元和藍色的新增單元Fig.2 Boundary motion on a fixed grid from timetntotn+1.Dead(red)and fresh(blue)cells are generated by the motion

因此，設計質量守恒的銳利浸入邊界法時，應當充分考慮網格拓撲的變化對流場的影響。與虛擬網格法不同，截斷單元法可以實時計算網格變化，消除了質量損失，但也造成了算法實現上的困難。在第2節至第4節中，將推導一個預估-校正格式，基于該格式和網格拓撲變化構造一個質量守恒且易于實現的銳利浸入邊界法。

2 強耦合預估-校正格式

2.1 流固耦合系統

一般情況下的流固耦合系統可以用一組抽象的狀態量描述：式中：S=S(t)和F=F(t)分別代表在t時刻固體和流體的一組狀態量，V被稱為流固耦合系統的狀態向量。從tn時刻到tn+1時刻的系統演化過程可以視為對狀態向量的一次映射。假設流固耦合系統狀態向量的所有映射均屬于某個向量空間，且均存在唯一的逆映射。記x為向量， ? 為映射。為了描述方便，使用符號?表示將映射作用于向量，約定如下記法：（1） ? ?x=?(x) ；（2） ?2??1?x=?2(?1(x))。注意兩個映射之間不一定滿足交換律，即 ?2??1?x≠?1??2?x。

假定流固耦合系統從tn時刻到tn+1時 刻遵從某個映射。記該映射為 Φ ，那么系統演化可表示為：

記逆映射Ψ =Φ?1，則：

寫成分塊矩陣形式：

可進一步對矩陣作LDU分解：

其中：

類比線性代數方程組，我們可以構造迭代格式。下面給出流固耦合系統的Gauss-Seidel迭代格式。

2.2 Gauss-Seidel格式

式(4)的Gauss-Seidel格式可表示為：

其中：

其中E為(D?L)的左逆，滿足：

其中I代表對角線元素均為恒等映射的單位矩陣。易得：

則：

于是有：

寫成如下形式：

其中：

根據S和F的物理含義，可認為G12對應于固體求解器，代表固體對自身的貢獻；G22對應于流體求解器，代表流體對自身的貢獻；G11代 表浸入邊界算法中的流體對固體的貢獻；G21代表浸入邊界算法中的固體對流體的貢獻。

將式(15)整理為如下形式：

（1）給定迭代初始值，

式(16)描述了一個強耦合的預估-校正算法流程。第3節中將依次闡述流體求解器和固體求解器，然后闡述浸入邊界的耦合力算和流體的校正算法。

3 控制方程和離散格式

本節簡要闡述流體和固體的控制方程與相應的離散格式。

3.1 流體控制方程和離散格式

二維黏性可壓縮流體方程組為：

其中：

熱通量qx和qy由 傅 里 葉定律給出：

式中： μ 為動力黏性系數；Pr為熱傳導系數；T=c2，c為聲速； γ 為多方指數。

對于非黏性項，我們使用帶有minmod限制器的二階MUSCL[35]有限體積格式，通量F和G采用基于Zha-Bilgen分裂[36]的AUSM+-up[37]方法。對于黏性項則采取中心差分格式。最后使用三階Runge-Kutta格式推進時間步。

3.2 固體控制方程和離散格式

3.2.1 剛體

剛體的Newton-Euler方程為：

式中：m為剛體質量；us為質心速度；fs為質心合力；fg為重力，在本文中均忽略重力影響；Is為質心慣性張量； ω 為角速度；Ts為質心合扭矩。

使用一階的顯式時間差分格式求解Newton-Euler方程，其中的合力和合扭矩來自流固界面上的耦合力。

3.2.2 線彈性體

線彈性體的控制方程及邊界條件為：

式中：σ為固體的應力；fV為體積力；ρs為固體密度；μd為阻尼系數，在本文中均假定μd=0；ds為位移；ε為應變；M為彈性矩陣；Sd和Sσ分別為位移邊界和應力邊界。

采取經典的有限元法對控制方程進行離散，并用Newmark-β隱式格式[38]推進時間步。

3.2.3 全局時間步長

計算流體的時間步長公式可以表示為：

式中：C為Courant數，l為網格最小尺寸，i為單元編號，ui和ci分別為單元i的速度和聲速。對于固體，除了固體求解器本身限制的最大時間步長 ?ts外，還必須考慮浸入邊界在網格上移動引起的變化。一般限制單步之內浸入邊界的移動不超過一個網格最小尺寸，故定義邊界移動時間步長：

式中： Γ 代表浸入邊界，常數Cb=0.95 。那么耦合場的全局時間步長為：

4 浸入邊界算法

本節將具體闡述流體對浸入邊界的耦合力算法和流體預估-校正算法中的邊界處理方法，并解釋算法保證質量守恒的原理。

4.1 耦合力

作用于固體表面的耦合力來自于流體的壓力。假設流體的壓力為p，浸入邊界的外法向量（指向流體一側）為n，在浸入邊界 Γ 上對表面的壓力積分可得到耦合力，即：

當 Γ 為直線段時，可采用高斯積分計算。關于高斯積分需要注意的是如果浸入邊界的長度遠大于流體網格尺寸，應先將浸入邊界分割為若干小段以改善精度。同時，在選取高斯積分的積分點時，為了減小邊界附近不光滑的壓力分布造成的壓力損失，可將積分點xΓ沿外法向量平移2h，即xΓ+2hn位置。

4.2 流體預估

式(17)中的預估步驟為：

式(28)對應的算法是將流體求解器應用于Fn。在式(28)中，固體的狀態量S不影響Fn的演化，因此將流體求解器應用于Fn時，將流固耦合邊界視為自由面，固體原本占據的空間初始化為零質量的單元，允許流體自由穿過耦合邊界。

4.3 流體校正

4.3.1 基本思路

考慮一個浸入流體域的固體，其占據空間為 ?s，邊 界 為Γs。流 體 被 固 體 擠 占 一 部 分 空 間后，該部分空間的網格單元上的質量為零相當于真空。本文中，將這部分真空域記為 ?v，邊界為Γv。如圖3所示，在tn時刻的真空域邊界與固體域邊界重合。在tn+1時刻二者分別演化為和，對應的內部空間為和。

圖3 流固耦合系統的子域變化（實線代表tn時刻的邊界，虛線代表tn+1時刻的邊界）Fig.3 Changes of the subdomains of the fluid-structure interaction system,where solid lines are boundaries attn ,and dashed lines are boundaries attn+1

4.3.2 輸運算法

如圖4所示，輸運算法的第1個步驟是使用感染-免疫法[39]標記單元。為直觀起見，用顏色代表不同標記。首先標記靠近浸入邊界外側的一層單元為藍色，接著標記與藍色單元相鄰的位于內側的一層單元為黃色，已標記的單元不再參與標記。以此類推，逐層向內染色。

圖4 流體標記和輸運方向，曲線為浸入邊界Fig.4 Fluid markers and direction of transportation,the curve is the immersed boundary

第2個步驟是根據染色順序移動流體。考慮一個紅色單元向橙色單元輸運的過程。記單元的守恒量w=(ρ,ρu,ρHt)T。假設輸運前紅色單元的守恒量為wr，橙色單元的守恒量為wo。僅考慮均勻正交網格上的有限體積法格式，則輸運導致單元的守恒量變為：

式中： λ 為輸運比例系數。假設該紅色單元周圍存在n個相鄰的橙色單元，其中n1個單元與該單元共享一條邊，稱為正相鄰單元，n?n1個單元與該單元共享一個頂點，稱為斜相鄰單元。進一步假設輸運比例與距離成反比，則在均勻正交網格上正相鄰單元和斜相鄰單元的輸運比例系數分別為：

輸運后紅色單元的守恒量變為零。以此類推，按照染色順序依次將流體從內側單元輸運到外側，直到內側單元流體完全轉移到藍色單元中，如圖5所示。圖中曲線為浸入邊界， λ1和 λ2為輸運比例系數。

圖5 依據染色順序逐層輸運流體Fig.5 Transport the fluid in sequence of colors

第3個步驟是校正藍色單元的速度。單元的守恒量w與物理量( ρ,u,e)的域之間存在雙射關系。假設當前藍色單元的守恒量為，則唯一導出當前密度、速度和內能為：

式中：w為的第k個分量。假設藍色單元的中心向浸入邊界的投影位置為點B，記點B的速度為uB。根據無滑移條件，將藍色單元的速度替換為點B的速度，保持密度和內能不變，重新導出守恒量：

至此，所有單元的校正結束。

上述的推導表明，該方法基于網格拓撲變化，無需構造虛擬網格。式(28)～(31)描述了校正步驟中的流體輸運規則，本質上是守恒律和無滑移條件的應用，因此是一種嚴格保證質量守恒的方法。與同樣保證質量守恒的截斷單元法相比，該方法的第一個步驟中的標記操作僅需要獲取網格單元與浸入邊界的相對距離，無需計算截斷單元的體積變化，簡化了算法實現。

5 數值算例

5.1 一維活塞問題

為了考察方法的準確性和收斂性，對Monasse等[31]提出的一維活塞問題進行了數值模擬。假設沿x軸方向在[0,7]區間放置一個一維管道，并放置一個密度為2，長度為0.5的剛性活塞，其中心位于x=2處。在(2.25,5)內填充理想氣體，密度為1，壓力為1 05，速度為零，多方指數為1.4。在(0,1.75)∪(5,7)內填充相同的理想氣體，密度和壓力分別改為10和1 06。計算域兩側設定為周期邊界條件。網格數N=1 440。

初始時刻氣體和活塞速度均為零。兩側不平衡的氣體壓力推動活塞向右運動，產生向右的壓縮波。利用移動最小二乘虛擬網格法和預估-校正方法分別計算，并與Monasse使用截斷單元法的結果進行對比。圖6展示了t=0.003時刻的壓力分布曲線，表明3種方法計算結果較接近。虛擬網格法在左側靠近浸入邊界的壓力和右側波峰壓力均偏高。預估-校正方法在左側與截斷單元法非常吻合，在右側波峰則偏低一些。

圖7展示了t=0.003之前的流體相對質量的變化。虛擬網格法計算的流場質量經歷了先減少后增加的過程，并且由于網格拓撲的突變而產生了振蕩。預估-校正方法計算的流場質量始終不變。通過對比，證明了預估-校正方法具有精確維持質量守恒的優勢。

圖7 流體相對質量的歷史曲線Fig.7 History of the relative mass of the fluid

為檢驗收斂性，使用如下公式計算壓力和速度分布的無量綱L2范數誤差：

圖6t=0.003時刻的壓力分布（網格數為1 440，虛線表示浸入邊界）Fig.6 Pressure distribution att=0.003(The number of cells is 1 440,and the dashed lines stand for the immersed boundaries.)

式中： ?i為壓力或速度，i表示分布曲線的數據點編號，為參考分布， ?0為 特征量。此處采用網格數N=7 290的計算結果作為參考分布，選取壓力和速度的特征量分別為p0=106和u0=300 。圖8展示了預估-校正方法計算的流體壓力和速度的誤差收斂曲線，證明壓力和速度分布曲線關于流體網格數N呈一階收斂。

圖8 壓力和速度的無量綱L2范數誤差Fig.8 Dimensionless L2 norms of error of the pressure and velocity

圖9展示了預估-校正方法計算的流場密度關于時空坐標(x,t)的云圖。圖9中深藍色條形區域代表活塞的歷史軌跡，反映了活塞加速和減速的過程，與Monasse使用截斷單元浸入邊界法的計算結果也基本符合。

圖9 兩種方法獲得的流場密度 ρ (x,t)云圖Fig.9 Mass density ρ (x,t)contours of the fluid by two methods

總體而言，該算例證明了本方法在一維剛體的流固耦合問題中能取得較準確的結果，并證明了該方法是一階收斂的。

5.2 激波沖擊平板

為了檢驗該方法在彈性結構流固耦合問題中的準確性，開展了相應的實驗和數值模擬。考慮一個激波與彈性平板耦合問題。如圖10所示，假設存在靜止的二維均勻流場，計算域為(?120,60)×(0,65) ，單位為mm，下同。在( 0,4)×(0,50)區域內放置一塊有機玻璃板，下方邊緣固定。流場介質為多方指數 γ =1.4 的理想氣體，初始流場的密度為1.29 kg/ m3，壓力為101.3 kPa。有機玻璃板的密度為1 180 kg/ m3，彈性模量為3.1 GPa，泊松比為0.4。在流場右邊界輸入一個馬赫數為1.2的激波。

圖10 激波管初始條件（藍色部分為有機玻璃板，灰色部分為靜止流場，右側紅色部分為輸入邊界）Fig.10 Initial conditions of the shock tube(The PMMA panel is blue,the static fluid is grey,the right boundary in red color is the inflow.)

實驗裝置的主體為一個激波管，實驗段如圖11所示。在激波管實驗段的下壁面安裝一個剛性基座，并在基座靠近前端位置插入平板。實驗段內初始氣體為空氣，激波的驅動氣體為氮氣。實驗期間使用高速相機拍攝光學紋影。

圖11 激波管實驗段（底部的金屬方塊為基座，豎立薄片為實驗測量的平板）Fig.11 Experimental section of the shock tube(The metal block on the bottom is the base,and the vertical sheet is the tested panel.)

在數值模擬中，對模型進行了簡化，略去了基座部分。流場使用均勻正交網格，單元密度為3 mm?1。平板的有限元網格包含253個Q4單元。局部網格如圖12所示。

圖12 局部網格（灰色部分為流場，藍色部分為平板）Fig.12 Local grids(The fluid is grey,and the panel is blue.)

圖13展示了實驗光學紋影與數值模擬紋影的對比。各時刻對比圖中上方為實驗紋影，下方為模擬紋影，即密度梯度云圖。模擬結果表明，激波向左運動一段時間后接觸平板，產生反射和繞射。繞流在平板上端右側頂點形成一個較大的渦。繞射后的激波改變方向，在平板后方的下壁面發生反射。在190～456 μs期間，繞流產生的渦逐漸擴大半徑并顯示出脫落的趨勢。在456 μs時刻下壁面反射波與逐漸擴大的渦接觸。通過對比實驗紋影，證明了模擬結果基本正確。圖14展示了在0～456 μs時間內實驗測量的平板的最大撓度與數值模擬結果的對比，證明使用的數值模擬方法能夠較準確地預測平板在激波作用下的變形。

圖13 不同時刻實驗紋影與模擬紋影的對比Fig.13 Comparison of experimental and simulated

圖14 平板的最大撓度Fig.14 Maximum deflections of the panel

值得注意的是，圖13中實驗拍攝和模擬的激波速度存在小幅度的偏差，推測形成這一偏差的主要原因可能來自于實驗的誤差。理想條件下，在實驗過程中可以通過調節壓縮參數改變激波管產生的激波，達到預設的馬赫數。但實際操作時會有一些不穩定因素，比如當天的氣溫較低、導致空氣介質密度和壓力變化、管道輕微漏氣等，這些因素使得真實的馬赫數偏離預設的馬赫數。此外，由于實驗紋影照片分辨率的限制，根據實驗紋影照片來確定時間零點也會造成測量誤差，因此該實驗方案仍有改進的空間。

6 結 論

討論了銳利浸入邊界法中的質量不守恒問題，指出守恒性與網格拓撲相關，提出了一種基于網格拓撲的預估-校正浸入邊界法。首先建立了流固耦合系統的演化模型，推導了強耦合的預估-校正格式。接著闡述了求解器和邊界處理方法，分析了流固耦合系統各介質的子域及其邊界的變化過程，給出了校正算法的具體步驟。對校正算法的推導過程表明，該算法不需要構造虛擬網格，因此避免了虛擬網格法引起的質量誤差。同時，該算法不計算截斷單元的體積變化，顯著簡化了算法的實現。

用該方法分別計算了一維剛體和二維線彈性平板的流固耦合問題，并與虛擬網格法和截斷單元法的結果進行對比。計算結果表明，該方法在一維剛體和二維線彈性平板的流固耦合問題中均能取得較準確的結果。在一維問題中，該方法嚴格地維持了全流場質量守恒，與理論預測一致；在二維問題中，該方法準確預測了激波與平板的相互作用過程，與實驗測量結果基本吻合。總體而言，該預估-校正浸入邊界法相比于傳統的虛擬網格法和截斷單元法具有精確質量守恒和易于實現的優勢，可作為一種替代方法。其主要缺陷在于一階精度相對較低，有待后續研究完善。