張和濤,寧建國,許香照,馬天寶
(北京理工大學爆炸科學與技術(shù)國家重點實驗室,北京 100081)
計算力學領(lǐng)域中的流固耦合問題具有廣闊的應用背景,長期以來吸引著眾多研究者關(guān)注。流固耦合仿真的理論和數(shù)值模擬一直是計算力學領(lǐng)域中難度較大的課題之一。浸入邊界法由Peskin[1-2]于1972年首次提出,目前已發(fā)展成為流固耦合數(shù)值模擬中廣泛采用的一類方法[3],被應用于研究聲場[4-6]、渦激振動[7-10]、納米流體[11]、湍流[12]等復雜問題。典型的浸入邊界法采用獨立的拉格朗日網(wǎng)格描述固體結(jié)構(gòu)并將之嵌入歐拉網(wǎng)格描述的流場中。通過對邊界附近進行特殊處理,浸入邊界法確保了流固耦合邊界上的無滑移條件和力平衡關(guān)系得以滿足。
浸入邊界法的一個分支是銳利浸入邊界法[13]。銳利浸入邊界法嚴格保持邊界位置上的強間斷,根據(jù)邊界兩側(cè)的應力推斷相互作用能夠始終維持銳利邊界。借助于邊界單元的特殊處理規(guī)則,銳利浸入邊界法通常無需在流體Navier-Stokes方程中添加附加源項,因此對格式穩(wěn)定性影響較小,有利于推廣到可壓縮流體情形。理論上,銳利浸入邊界法可以具有高階精度[14],在模擬高雷諾數(shù)流動時相比浸入邊界法的另一分支模糊浸入邊界法[15-16]具有優(yōu)勢。銳利浸入邊界法的核心問題是處理被邊界截斷的單元,簡稱截斷單元(cut-cell)。主流實現(xiàn)方式分為兩類:虛擬網(wǎng)格法(ghost-cell method)[17-28]和截斷單元法(cutcell method)[29-33],其中應用相對較多的是虛擬網(wǎng)格法。對虛擬網(wǎng)格法的研究表明,為了處理截斷單元,虛擬網(wǎng)格法需要進行多維重構(gòu),常用方法包括形函數(shù)插值和移動最小二乘法等。
虛擬網(wǎng)格法存在的一個問題是當耦合邊界在流體的固定網(wǎng)格上移動時不能維持嚴格的質(zhì)量守恒[23]。邊界運動時會改變固體的覆蓋域,在流體網(wǎng)格上產(chǎn)生許多失效單元和新增單元。被更新后的固體覆蓋的失效單元會損失其中包含的所有信息并造成失效單元質(zhì)量丟失。為了獲得正確的壓力,虛擬網(wǎng)格法有時需要對新增單元進行填充,導致新增單元獲得額外質(zhì)量。因此整體上虛擬網(wǎng)格法不能維持質(zhì)量守恒。Seo等[34]的研究確認違反幾何守恒律及相應的質(zhì)量損失是造成虛假振蕩的首要原因,并提出了一種針對不可壓縮流體的改善方案。目前對可壓縮流體仍無合適的解決方法。
截斷單元法考慮了截斷單元體積變化的過程,因此通過合適的單元處理可以實現(xiàn)精確的質(zhì)量守恒。Monasse等[31]提出了一種針對可壓縮流體和剛體耦合的守恒型浸入邊界法,通過校正Godunov通量來還原損失的質(zhì)量,但由于采用了有限體積法,涉及截斷單元的幾何計算會比較困難。與Monasse的方法類似,Sotiropoulos等[13]給出的結(jié)論表明,在截斷單元法中對截斷單元的拓撲管理是主要問題,妨礙了這一類方法推廣到三維情形。Brady等[33]指出,當截斷單元體積趨于零時會引起系統(tǒng)剛性問題。用有限差分法替代有限體積法可以顯著簡化截斷單元的計算。然而由于有限體積法在計算流體力學領(lǐng)域的廣泛應用,考慮有限體積法的耦合問題仍然是必要的。
為了克服虛擬網(wǎng)格法和截斷單元法存在的困難,本文中,將提出一種區(qū)別于這兩類方法的新思路,即適用于可壓縮流體和結(jié)構(gòu)耦合問題的強耦合預估-校正銳利浸入邊界法。首先,分析虛擬網(wǎng)格浸入邊界法導致質(zhì)量不守恒的原因,揭示網(wǎng)格拓撲變化的影響。其次,闡述一般流固耦合系統(tǒng)的矩陣表示,導出強耦合Gauss-Seidel迭代格式,進一步導出預估-校正格式,基于該格式,提出預估-校正算法,預估步驟使用無耦合邊界的流體模型求解,并對流體施加校正。然后,介紹獨立的介質(zhì)求解器,并詳細闡述校正步驟,在校正步驟中,通過一組輸運規(guī)則來實現(xiàn)輸運過程,保證質(zhì)量守恒。最后,使用該方法測試若干數(shù)值算例,檢驗該方法的準確性和收斂性。
有限體積法在重構(gòu)流場的通量函數(shù)時需要使用相鄰多個單元組成的模板。當通量模板包含至少一個截斷單元時,一般采取構(gòu)造虛擬網(wǎng)格的方式補全模板。一種常用的方法是使用基于移動最小二乘法的虛擬網(wǎng)格法補全模板。如圖1所示,下方的4個單元組成一個通量模板。單元G被浸入的固體網(wǎng)格覆蓋。斜實線代表浸入邊界。點G關(guān)于浸入邊界的鏡像為點I,兩點連線與浸入邊界相交于點H。距點I一定半徑范圍內(nèi)的單元組成一個支撐域 ?I。在 ?I上使用徑向基函數(shù)移動最小二乘法擬合點I的物理量。

圖1 根據(jù)鏡像點I構(gòu)造虛擬網(wǎng)格結(jié)點GFig.1 Construct the ghost pointGbased on its image pointI
根據(jù)點I的物理量和浸入邊界上的邊界條件可構(gòu)造虛擬結(jié)點G。為了與固定邊界情形下的邊界條件一致,假設(shè)I和G關(guān)于H參考系滿足固壁邊界條件,可表達為:

式中: ρ 為密度,Ht為總能,u為速度。該邊界條件的兩種情形分別對應于Dirichlet邊界條件和Neumann邊界條件。
虛擬網(wǎng)格法會導致流體質(zhì)量的改變。如圖2所示,當浸入邊界在固定的歐拉網(wǎng)格上移動時,會覆蓋一部分網(wǎng)格從而造成失效單元,同時一部分網(wǎng)格不再被覆蓋而成為新增單元。虛擬網(wǎng)格法在構(gòu)造虛擬網(wǎng)格時沒有考慮到網(wǎng)格拓撲變化的影響。當出現(xiàn)失效單元時,網(wǎng)格拓撲發(fā)生突變,與連續(xù)的邊界條件產(chǎn)生矛盾。雖然施加了邊界條件,但是仍然會損失一部分質(zhì)量。當出現(xiàn)新增單元時,為了擬合的光滑性,可能需要使用擬合方法對新增單元進行填充。同樣地,在新增單元上的擬合算法中加入邊界條件約束也不能保證質(zhì)量守恒。

圖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
因此,設(shè)計質(zhì)量守恒的銳利浸入邊界法時,應當充分考慮網(wǎng)格拓撲的變化對流場的影響。與虛擬網(wǎng)格法不同,截斷單元法可以實時計算網(wǎng)格變化,消除了質(zhì)量損失,但也造成了算法實現(xiàn)上的困難。在第2節(jié)至第4節(jié)中,將推導一個預估-校正格式,基于該格式和網(wǎng)格拓撲變化構(gòu)造一個質(zhì)量守恒且易于實現(xiàn)的銳利浸入邊界法。

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

記逆映射Ψ =Φ?1,則:

寫成分塊矩陣形式:

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

其中:

類比線性代數(shù)方程組,我們可以構(gòu)造迭代格式。下面給出流固耦合系統(tǒng)的Gauss-Seidel迭代格式。
式(4)的Gauss-Seidel格式可表示為:

其中:

其中E為(D?L)的左逆,滿足:

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

則:

于是有:

寫成如下形式:

其中:

根據(jù)S和F的物理含義,可認為G12對應于固體求解器,代表固體對自身的貢獻;G22對應于流體求解器,代表流體對自身的貢獻;G11代 表浸入邊界算法中的流體對固體的貢獻;G21代表浸入邊界算法中的固體對流體的貢獻。
將式(15)整理為如下形式:

(1)給定迭代初始值,

式(16)描述了一個強耦合的預估-校正算法流程。第3節(jié)中將依次闡述流體求解器和固體求解器,然后闡述浸入邊界的耦合力算和流體的校正算法。
本節(jié)簡要闡述流體和固體的控制方程與相應的離散格式。
二維黏性可壓縮流體方程組為:

其中:



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

式中: μ 為動力黏性系數(shù);Pr為熱傳導系數(shù);T=c2,c為聲速; γ 為多方指數(shù)。
對于非黏性項,我們使用帶有minmod限制器的二階MUSCL[35]有限體積格式,通量F和G采用基于Zha-Bilgen分裂[36]的AUSM+-up[37]方法。對于黏性項則采取中心差分格式。最后使用三階Runge-Kutta格式推進時間步。
3.2.1 剛體
剛體的Newton-Euler方程為:

式中:m為剛體質(zhì)量;us為質(zhì)心速度;fs為質(zhì)心合力;fg為重力,在本文中均忽略重力影響;Is為質(zhì)心慣性張量; ω 為角速度;Ts為質(zhì)心合扭矩。
使用一階的顯式時間差分格式求解Newton-Euler方程,其中的合力和合扭矩來自流固界面上的耦合力。
3.2.2 線彈性體
線彈性體的控制方程及邊界條件為:

式中:σ為固體的應力;fV為體積力;ρs為固體密度;μd為阻尼系數(shù),在本文中均假定μd=0;ds為位移;ε為應變;M為彈性矩陣;Sd和Sσ分別為位移邊界和應力邊界。
采取經(jīng)典的有限元法對控制方程進行離散,并用Newmark-β隱式格式[38]推進時間步。
3.2.3 全局時間步長
計算流體的時間步長公式可以表示為:

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

式中: Γ 代表浸入邊界,常數(shù)Cb=0.95 。那么耦合場的全局時間步長為:

本節(jié)將具體闡述流體對浸入邊界的耦合力算法和流體預估-校正算法中的邊界處理方法,并解釋算法保證質(zhì)量守恒的原理。
作用于固體表面的耦合力來自于流體的壓力。假設(shè)流體的壓力為p,浸入邊界的外法向量(指向流體一側(cè))為n,在浸入邊界 Γ 上對表面的壓力積分可得到耦合力,即:

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

式(28)對應的算法是將流體求解器應用于Fn。在式(28)中,固體的狀態(tài)量S不影響Fn的演化,因此將流體求解器應用于Fn時,將流固耦合邊界視為自由面,固體原本占據(jù)的空間初始化為零質(zhì)量的單元,允許流體自由穿過耦合邊界。
4.3.1 基本思路
考慮一個浸入流體域的固體,其占據(jù)空間為 ?s,邊 界 為Γs。流 體 被 固 體 擠 占 一 部 分 空 間后,該部分空間的網(wǎng)格單元上的質(zhì)量為零相當于真空。本文中,將這部分真空域記為 ?v,邊界為Γv。如圖3所示,在tn時刻的真空域邊界與固體域邊界重合。在tn+1時刻二者分別演化為和,對應的內(nèi)部空間為和。

圖3 流固耦合系統(tǒng)的子域變化(實線代表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]標記單元。為直觀起見,用顏色代表不同標記。首先標記靠近浸入邊界外側(cè)的一層單元為藍色,接著標記與藍色單元相鄰的位于內(nèi)側(cè)的一層單元為黃色,已標記的單元不再參與標記。以此類推,逐層向內(nèi)染色。

圖4 流體標記和輸運方向,曲線為浸入邊界Fig.4 Fluid markers and direction of transportation,the curve is the immersed boundary
第2個步驟是根據(jù)染色順序移動流體。考慮一個紅色單元向橙色單元輸運的過程。記單元的守恒量w=(ρ,ρu,ρHt)T。假設(shè)輸運前紅色單元的守恒量為wr,橙色單元的守恒量為wo。僅考慮均勻正交網(wǎng)格上的有限體積法格式,則輸運導致單元的守恒量變?yōu)椋?/p>

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

輸運后紅色單元的守恒量變?yōu)榱恪R源祟愅疲凑杖旧樞蛞来螌⒘黧w從內(nèi)側(cè)單元輸運到外側(cè),直到內(nèi)側(cè)單元流體完全轉(zhuǎn)移到藍色單元中,如圖5所示。圖中曲線為浸入邊界, λ1和 λ2為輸運比例系數(shù)。

圖5 依據(jù)染色順序逐層輸運流體Fig.5 Transport the fluid in sequence of colors
第3個步驟是校正藍色單元的速度。單元的守恒量w與物理量( ρ,u,e)的域之間存在雙射關(guān)系。假設(shè)當前藍色單元的守恒量為,則唯一導出當前密度、速度和內(nèi)能為:

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

至此,所有單元的校正結(jié)束。
上述的推導表明,該方法基于網(wǎng)格拓撲變化,無需構(gòu)造虛擬網(wǎng)格。式(28)~(31)描述了校正步驟中的流體輸運規(guī)則,本質(zhì)上是守恒律和無滑移條件的應用,因此是一種嚴格保證質(zhì)量守恒的方法。與同樣保證質(zhì)量守恒的截斷單元法相比,該方法的第一個步驟中的標記操作僅需要獲取網(wǎng)格單元與浸入邊界的相對距離,無需計算截斷單元的體積變化,簡化了算法實現(xiàn)。
為了考察方法的準確性和收斂性,對Monasse等[31]提出的一維活塞問題進行了數(shù)值模擬。假設(shè)沿x軸方向在[0,7]區(qū)間放置一個一維管道,并放置一個密度為2,長度為0.5的剛性活塞,其中心位于x=2處。在(2.25,5)內(nèi)填充理想氣體,密度為1,壓力為1 05,速度為零,多方指數(shù)為1.4。在(0,1.75)∪(5,7)內(nèi)填充相同的理想氣體,密度和壓力分別改為10和1 06。計算域兩側(cè)設(shè)定為周期邊界條件。網(wǎng)格數(shù)N=1 440。
初始時刻氣體和活塞速度均為零。兩側(cè)不平衡的氣體壓力推動活塞向右運動,產(chǎn)生向右的壓縮波。利用移動最小二乘虛擬網(wǎng)格法和預估-校正方法分別計算,并與Monasse使用截斷單元法的結(jié)果進行對比。圖6展示了t=0.003時刻的壓力分布曲線,表明3種方法計算結(jié)果較接近。虛擬網(wǎng)格法在左側(cè)靠近浸入邊界的壓力和右側(cè)波峰壓力均偏高。預估-校正方法在左側(cè)與截斷單元法非常吻合,在右側(cè)波峰則偏低一些。
圖7展示了t=0.003之前的流體相對質(zhì)量的變化。虛擬網(wǎng)格法計算的流場質(zhì)量經(jīng)歷了先減少后增加的過程,并且由于網(wǎng)格拓撲的突變而產(chǎn)生了振蕩。預估-校正方法計算的流場質(zhì)量始終不變。通過對比,證明了預估-校正方法具有精確維持質(zhì)量守恒的優(yōu)勢。

圖7 流體相對質(zhì)量的歷史曲線Fig.7 History of the relative mass of the fluid
為檢驗收斂性,使用如下公式計算壓力和速度分布的無量綱L2范數(shù)誤差:

圖6t=0.003時刻的壓力分布(網(wǎng)格數(shù)為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表示分布曲線的數(shù)據(jù)點編號,為參考分布, ?0為 特征量。此處采用網(wǎng)格數(shù)N=7 290的計算結(jié)果作為參考分布,選取壓力和速度的特征量分別為p0=106和u0=300 。圖8展示了預估-校正方法計算的流體壓力和速度的誤差收斂曲線,證明壓力和速度分布曲線關(guān)于流體網(wǎng)格數(shù)N呈一階收斂。

圖8 壓力和速度的無量綱L2范數(shù)誤差Fig.8 Dimensionless L2 norms of error of the pressure and velocity
圖9展示了預估-校正方法計算的流場密度關(guān)于時空坐標(x,t)的云圖。圖9中深藍色條形區(qū)域代表活塞的歷史軌跡,反映了活塞加速和減速的過程,與Monasse使用截斷單元浸入邊界法的計算結(jié)果也基本符合。

圖9 兩種方法獲得的流場密度 ρ (x,t)云圖Fig.9 Mass density ρ (x,t)contours of the fluid by two methods
總體而言,該算例證明了本方法在一維剛體的流固耦合問題中能取得較準確的結(jié)果,并證明了該方法是一階收斂的。
為了檢驗該方法在彈性結(jié)構(gòu)流固耦合問題中的準確性,開展了相應的實驗和數(shù)值模擬。考慮一個激波與彈性平板耦合問題。如圖10所示,假設(shè)存在靜止的二維均勻流場,計算域為(?120,60)×(0,65) ,單位為mm,下同。在( 0,4)×(0,50)區(qū)域內(nèi)放置一塊有機玻璃板,下方邊緣固定。流場介質(zhì)為多方指數(shù) γ =1.4 的理想氣體,初始流場的密度為1.29 kg/ m3,壓力為101.3 kPa。有機玻璃板的密度為1 180 kg/ m3,彈性模量為3.1 GPa,泊松比為0.4。在流場右邊界輸入一個馬赫數(shù)為1.2的激波。

圖10 激波管初始條件(藍色部分為有機玻璃板,灰色部分為靜止流場,右側(cè)紅色部分為輸入邊界)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所示。在激波管實驗段的下壁面安裝一個剛性基座,并在基座靠近前端位置插入平板。實驗段內(nèi)初始氣體為空氣,激波的驅(qū)動氣體為氮氣。實驗期間使用高速相機拍攝光學紋影。

圖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.)
在數(shù)值模擬中,對模型進行了簡化,略去了基座部分。流場使用均勻正交網(wǎng)格,單元密度為3 mm?1。平板的有限元網(wǎng)格包含253個Q4單元。局部網(wǎng)格如圖12所示。

圖12 局部網(wǎng)格(灰色部分為流場,藍色部分為平板)Fig.12 Local grids(The fluid is grey,and the panel is blue.)
圖13展示了實驗光學紋影與數(shù)值模擬紋影的對比。各時刻對比圖中上方為實驗紋影,下方為模擬紋影,即密度梯度云圖。模擬結(jié)果表明,激波向左運動一段時間后接觸平板,產(chǎn)生反射和繞射。繞流在平板上端右側(cè)頂點形成一個較大的渦。繞射后的激波改變方向,在平板后方的下壁面發(fā)生反射。在190~456 μs期間,繞流產(chǎn)生的渦逐漸擴大半徑并顯示出脫落的趨勢。在456 μs時刻下壁面反射波與逐漸擴大的渦接觸。通過對比實驗紋影,證明了模擬結(jié)果基本正確。圖14展示了在0~456 μs時間內(nèi)實驗測量的平板的最大撓度與數(shù)值模擬結(jié)果的對比,證明使用的數(shù)值模擬方法能夠較準確地預測平板在激波作用下的變形。

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

圖14 平板的最大撓度Fig.14 Maximum deflections of the panel
值得注意的是,圖13中實驗拍攝和模擬的激波速度存在小幅度的偏差,推測形成這一偏差的主要原因可能來自于實驗的誤差。理想條件下,在實驗過程中可以通過調(diào)節(jié)壓縮參數(shù)改變激波管產(chǎn)生的激波,達到預設(shè)的馬赫數(shù)。但實際操作時會有一些不穩(wěn)定因素,比如當天的氣溫較低、導致空氣介質(zhì)密度和壓力變化、管道輕微漏氣等,這些因素使得真實的馬赫數(shù)偏離預設(shè)的馬赫數(shù)。此外,由于實驗紋影照片分辨率的限制,根據(jù)實驗紋影照片來確定時間零點也會造成測量誤差,因此該實驗方案仍有改進的空間。
討論了銳利浸入邊界法中的質(zhì)量不守恒問題,指出守恒性與網(wǎng)格拓撲相關(guān),提出了一種基于網(wǎng)格拓撲的預估-校正浸入邊界法。首先建立了流固耦合系統(tǒng)的演化模型,推導了強耦合的預估-校正格式。接著闡述了求解器和邊界處理方法,分析了流固耦合系統(tǒng)各介質(zhì)的子域及其邊界的變化過程,給出了校正算法的具體步驟。對校正算法的推導過程表明,該算法不需要構(gòu)造虛擬網(wǎng)格,因此避免了虛擬網(wǎng)格法引起的質(zhì)量誤差。同時,該算法不計算截斷單元的體積變化,顯著簡化了算法的實現(xiàn)。
用該方法分別計算了一維剛體和二維線彈性平板的流固耦合問題,并與虛擬網(wǎng)格法和截斷單元法的結(jié)果進行對比。計算結(jié)果表明,該方法在一維剛體和二維線彈性平板的流固耦合問題中均能取得較準確的結(jié)果。在一維問題中,該方法嚴格地維持了全流場質(zhì)量守恒,與理論預測一致;在二維問題中,該方法準確預測了激波與平板的相互作用過程,與實驗測量結(jié)果基本吻合。總體而言,該預估-校正浸入邊界法相比于傳統(tǒng)的虛擬網(wǎng)格法和截斷單元法具有精確質(zhì)量守恒和易于實現(xiàn)的優(yōu)勢,可作為一種替代方法。其主要缺陷在于一階精度相對較低,有待后續(xù)研究完善。