王婷婷，呂宏強(qiáng)，安 慰

(南京航空航天大學(xué) 航空學(xué)院，南京 210016)

0 引言

浸沒邊界法（immersed boundary method,IBM）最初由Peskin 于1972 年[1]提出，用于心臟瓣膜周圍流動的模擬。在該方法中[1-2]，流體在笛卡爾網(wǎng)格上求解，幾何體被一組拉格朗日點(diǎn)所取代，通過在控制方程中添加源項(xiàng)實(shí)現(xiàn)邊界對流體的作用，如Lee[3]和Nicolaou[4]所述，這類IBM 稱為連續(xù)強(qiáng)迫IBM。自從IBM 提出以來，該方法得到了巨大的發(fā)展并廣泛應(yīng)用。Mohd-Yusof[5]和Fadlun 等[6]提出了直接強(qiáng)迫IBM，與連續(xù)強(qiáng)迫IBM 不同的是，不使用拉格朗日點(diǎn)并且強(qiáng)迫項(xiàng)在笛卡爾網(wǎng)格上直接進(jìn)行計(jì)算，控制方程僅僅在流體域的計(jì)算網(wǎng)格上求解。何躍龍等[7]引入Choi 等[8]提出的指數(shù)插值法對目標(biāo)點(diǎn)速度進(jìn)行插值，建立了適用于可壓縮流動的浸沒邊界法。Yang等[9]采用大渦模擬與浸沒邊界相結(jié)合的方法，對后向臺階流動進(jìn)行了數(shù)值模擬。Giannenas 和Laizet[10]提出了一種基于三次樣條重建的浸沒邊界方法，用于基于笛卡爾網(wǎng)格的湍流中浸沒物體的高保真模擬。Kou等[11]提出一種基于體積懲罰的結(jié)合高階通量重建方法的浸沒邊界法。盡管IBM 已得到廣泛應(yīng)用，但采用笛卡爾網(wǎng)格結(jié)合高階方法的探索仍較少。

近年來，間斷伽遼金（discontinuous Galerkin,DG）方法已成為CFD 中最有前途的高階離散技術(shù)之一[12-13]。1997 年，Bassi 和Rebay 等[14-15]采用DG 方法求解了歐拉方程和N-S 方程，隨后Bassi 等[16]又將DG 方法拓展到可壓縮流和不可壓縮流的DNS（direct numerical simulation）和ILES（implicit large eddy simulation）。Landmann[17]在其博士論文中發(fā)展了并行高階DG 方法并對二維N-S 方程及RANS（Reynolds-averaged Navier-Stokes）方程進(jìn)行了求解。DG 方法在國內(nèi)起步較晚，蔚喜軍和周鐵[18]在2005 年采用龍格-庫塔間斷伽遼金（Runge-Kutta discontinuous Galerkin,RKDG）方法求解了二維可壓縮歐拉方程。于劍和閻超[19]通過引入全局提升算子和局部提升算子，發(fā)展了求解N-S 方程的DG 方法的一般框架。張來平等[20]發(fā)展了DG/FV 混合方法，計(jì)算效率相較于同階精度的DG 方法更高。黨亞斌等[21]發(fā)展了一種基于解析精確Jacobian 矩陣的GMRES隱式方法，用于求解可壓縮層流及湍流問題，提高了隱式高精度DG 方法的計(jì)算穩(wěn)定性和計(jì)算效率。國外已有研究者將DG 方法與浸沒邊界法相結(jié)合，但大多采用切割網(wǎng)格的方法。針對基于笛卡爾網(wǎng)格的浸沒式幾何體，Qin 等[22]提出了一種求解歐拉方程的DG 方法。Müller等[23]提出了一種高階離散格式，用于求解具有浸沒邊界的可壓縮歐拉和N-S 方程，在一個(gè)由水平集函數(shù)隱式定義的域中使用DG 離散化，采用了一種非侵入式單元聚集技術(shù)。雖然切割網(wǎng)格法對于靜態(tài)邊界能夠較好的實(shí)現(xiàn)，然而切割網(wǎng)格法僅限于低精度，而且小切割單元的處理以及動態(tài)邊界的處理比較復(fù)雜。

本文引入Kou 等[11]提出的基于體積懲罰的浸沒邊界實(shí)現(xiàn)思路，提出基于懲罰項(xiàng)模式的高階DG 浸沒邊界法，采用典型算例對所提出方法進(jìn)行測試驗(yàn)證。

1 控制方程

可壓縮黏性流體的控制方程[17]寫成：

式中：U為守恒變量；Fi為無黏通量；Fv為黏性通量。

2 間斷伽遼金空間離散

因?yàn)槭剑?）中存在2 階偏導(dǎo)數(shù)項(xiàng)，在求解前必須對方程進(jìn)行降階處理。采用Landmann[17]使用的混合方法，在計(jì)算中把1 階偏導(dǎo)數(shù)項(xiàng)當(dāng)作額外變量，引入輔助變量 Θ。式（1）可以改寫成如下形式：

采用間斷伽遼金有限元法對式（2）和式（3）進(jìn)行空間離散。在式（2）和式（3）的兩邊同乘一個(gè)任意的測試函數(shù)W，并將測試函數(shù)的高階表達(dá)式代入，在計(jì)算域單元邊界上的積分中引入逆風(fēng)格式，得到如下方程組：

3 隱式時(shí)間離散

對式（4）和式（5）進(jìn)行數(shù)值積分之后，得到的方程組如下所示：

式中：u=[u1,u2,···,uk,···,uNele]T為全局自由度矢量，Nele為全局單元的總數(shù)，uk為單元k的自由度；R(u)=[R1,R2,···,Rk,···,RNele]T為全局殘值矢量，Rk為單元k的殘值矢量；M為全局塊質(zhì)量矩陣。

采用隱式時(shí)間離散方法，式（6）右端殘值項(xiàng)取新的時(shí)間步tn+1，計(jì)算得到：

采用向后差分方法處理之后，式（7）可以寫成：

由于求解的是非線性的控制方程，所以對式（8）的右端項(xiàng)進(jìn)行線性化：

式中，(?R/?u)k為需要求解的大型線性方程組的雅克比矩陣。在每個(gè)時(shí)間步tn+1內(nèi)采用牛頓迭代法進(jìn)行求解，對于定常算例，只需進(jìn)行一個(gè)牛頓步就可以達(dá)到需要的精度：

其中：k=(0,1,···,kmax)；c∈[0,1]為阻尼因子。在一個(gè)時(shí)間步內(nèi)即一個(gè)牛頓迭代循環(huán)內(nèi)，以第tn步的數(shù)值解作為初始值進(jìn)行牛頓迭代，直到式（11）中的殘值下降到設(shè)定閾值，將第kmax個(gè)牛頓步的結(jié)果作為新的時(shí)間步tn+1的數(shù)值解并進(jìn)入下一個(gè)時(shí)間步進(jìn)行計(jì)算。

4 浸沒邊界法

IBM 的基本思想是通過適當(dāng)?shù)臄?shù)值處理向非貼體網(wǎng)格施加邊界條件。體積懲罰是一種典型的方法，它通過懲罰固體中積分點(diǎn)的速度來施加邊界條件。

4.1 體積懲罰法

體積懲罰法通過在控制方程中引入懲罰源項(xiàng)來施加邊界條件。在這種方法中，首先定義一個(gè)用于區(qū)分流體區(qū)域和固體區(qū)域的界面函數(shù)：

該界面函數(shù)用于確定是否對當(dāng)前積分點(diǎn)施加IBM力。對于移動的邊界，χ(x,t)是隨時(shí)間而變化的。對于高階笛卡爾網(wǎng)格的體積懲罰方法（圖1），所有被固體區(qū)域覆蓋的積分點(diǎn)都需要進(jìn)行懲罰以施加邊界條件。

用IBM 處理的N-S 方程可以寫成：

式中，“S”為IBM 源項(xiàng)。

考慮固體內(nèi)部速度us=(us,vs)T狄利克雷邊界條件，源項(xiàng)為：

式中，η為體積懲罰參數(shù)?？紤]無滑移邊界條件時(shí)，設(shè)條件us=(0,0)T。本文懲罰項(xiàng)遵循文獻(xiàn)[25]中的公式，以施加絕熱壁的邊界條件。上述方程也可用于運(yùn)動物體，其中固體速度根據(jù)運(yùn)動方程進(jìn)行更新。

通過結(jié)合算子分裂和隱式時(shí)間離散方法，控制方程等式（14）在時(shí)間上進(jìn)行。使用2 階Strang 分裂[26]方法添加源項(xiàng)，正如Piquet 等[27]所討論的，通過Strang分裂可以精確計(jì)算動量和能量方程的懲罰項(xiàng)。在時(shí)間步驟n，執(zhí)行以下操作：

在數(shù)值計(jì)算中，將控制方程的變量進(jìn)行無量綱化處理[17]。下文中如無特別說明，懲罰參數(shù)均為η=1 × 10-3。在時(shí)間步驟n，執(zhí)行第2 步，即對式（17）進(jìn)行牛頓法迭代時(shí)，式（11）中的 無量綱時(shí)間步長 Δt=0.1。

4.2 邊界表示和界面函數(shù)

固體邊界的表示是IBM 方法的關(guān)鍵之一，其主要用來進(jìn)行界面函數(shù)的定義和氣動系數(shù)的計(jì)算。

浸沒邊界的離散化應(yīng)足夠靈活，以處理復(fù)雜的幾何形狀。在本文中，選擇使用一組拉格朗日標(biāo)記點(diǎn)來表示實(shí)體邊界，定義為浸沒邊界點(diǎn)。標(biāo)記點(diǎn)由線性元素（即二維線段）連接。使用該表示法可以有效地計(jì)算幾何量，包括曲面法線、用于數(shù)據(jù)重建的插值模板、曲面距離[28-29]。浸沒邊界由離散的點(diǎn)定義，點(diǎn)在網(wǎng)格中的位置用水平集方法可描述為滿足 φ=0，φ是到浸沒邊界的有符號距離[30]。

對于圓柱體、球體或機(jī)翼等簡單幾何體，存在解析形狀函數(shù)，可用于定義界面函數(shù)。例如，對于直徑為1 的圓，有：

4.3 曲面數(shù)據(jù)重建

為了從模擬結(jié)果中獲得相關(guān)物理量（如壓力、摩擦力）的分布，固體表面物理量的重建是一個(gè)重要工作。

在高階框架中，可以直接考慮使用每個(gè)單元中定義的高階多項(xiàng)式來執(zhí)行數(shù)據(jù)重建。然而，這種方法有一些潛在的缺點(diǎn)，可能導(dǎo)致結(jié)果不準(zhǔn)確。這主要是因?yàn)樵诖蠖鄶?shù)情況下，此類插值方法將涉及浸入固體中的積分點(diǎn)，而這些積分點(diǎn)通常具有非物理值。

除多項(xiàng)式插值外，另一種方法是從幾個(gè)最近的積分點(diǎn)插值數(shù)據(jù)。與基于高階多項(xiàng)式的插值相比，其具有更大的靈活性來選擇插值點(diǎn)（interpolation point,IP），而不涉及具有非物理值的固體點(diǎn)。這里，采用文獻(xiàn)[31]中描述的適用于高階框架的插值方法。通常，文獻(xiàn)[31]中的插值方法基于IB 點(diǎn)與插值點(diǎn)之間的逆距離。首先從IB 點(diǎn)周圍最近的流體點(diǎn)中選擇候選插值點(diǎn)，然后插值保守變量的值和梯度。本文采用插值點(diǎn)處的逆距離權(quán)重（inverse distance weight at interpolation point,IDW-IP）方法進(jìn)行插值，如圖2 所示。

圖2 IDW-IP 方法的插值模板示意圖Fig.2 Interpolation stencil for IDW-IP method

插值點(diǎn)在這種情況下被定義為一個(gè)點(diǎn)，位于IB 點(diǎn)的法線上，在這個(gè)點(diǎn)上，屬性從周圍的單元內(nèi)插，然后轉(zhuǎn)移到IB 點(diǎn)。關(guān)于插值點(diǎn)，首先要確定它的位置，其次是它的屬性。為了確定插值點(diǎn)的位置，使用如下關(guān)系式：

式中：d1為 A 點(diǎn)到IB 點(diǎn)法線的垂直距離；d2為A 點(diǎn)到IB 點(diǎn)的法線的距離投影；i為模板中場單元和帶單元的固體外積分點(diǎn)。為了確定插值點(diǎn)的壓力值，使用以下關(guān)系式：

在這種方法中，在插值點(diǎn)確定的壓力被復(fù)制到IB 點(diǎn)。

5 算例測試

在本節(jié)中，應(yīng)用不同測試算例驗(yàn)證所提出的基于高階間斷伽遼金數(shù)值方法的浸沒邊界法。

5.1 定常圓柱繞流

第一個(gè)測試案例為二維定常圓柱繞流，雷諾數(shù)和馬赫數(shù)分別設(shè)置為40 和0.2。壁面采用無滑移絕熱壁面邊界條件，因此在固體中施加us=(0,0)T。矩形計(jì)算域的大小為x∈[-5D,10D]、y∈[-5D,5D]，其 中D是圓柱的直徑，設(shè)置為1。在正方形區(qū)域x∈[-D,D]、y∈[-D,D]，使用均勻網(wǎng)格。均勻網(wǎng)格尺寸h分別選擇為0.03D（網(wǎng)格1）、0.05D（網(wǎng)格2）和0.1D（網(wǎng)格3）。網(wǎng)格1 如圖3 所示，局部均勻網(wǎng)格大小h=0.03D。

圖3 計(jì)算網(wǎng)格（網(wǎng)格1）Fig.3 Computational grid (grid 1)

在網(wǎng)格1 情況下，由0 階至3 階DG 計(jì)算得到的流線如圖4、物面附近馬赫數(shù)云圖如圖5所示。從圖中可以看出，隨著階數(shù)的提高，同一網(wǎng)格下計(jì)算結(jié)果的精度顯著提高，尤其是物面附近的數(shù)值解更加光滑。在網(wǎng)格1 情況下，不同階數(shù)計(jì)算得到的壓力系數(shù)Cp與文獻(xiàn)[11]中貼體網(wǎng)格模擬結(jié)果的對比如圖6 所示，圖中實(shí)線表示文獻(xiàn)[11]貼體網(wǎng)格的計(jì)算結(jié)果，符號表示本文方法數(shù)值模擬得到的結(jié)果（不同符號代表不同的階數(shù)）。結(jié)果表明，壓力系數(shù)分布隨著計(jì)算階數(shù)的提高與貼體網(wǎng)格結(jié)果具有良好的一致性。

圖4 網(wǎng)格1 下不同階數(shù)流線圖和馬赫數(shù)云圖Fig.4 Streamlines and Mach number contours calculated by DG with different orders (grid 1)

圖5 網(wǎng)格1 情況下不同階數(shù)物面附近馬赫云圖Fig.5 Mach number contours calculated by DG with different orders (grid 1)

不同網(wǎng)格尺度下壓力系數(shù)Cp的 對比如圖7 所示，圖中分別給出了0 階至3 階DG、物面附近網(wǎng)格尺度為0.03D、0.05D、0.1D時(shí)的結(jié)果對比。由圖7（a-d）所示，在同1 階數(shù)下，隨著網(wǎng)格尺度的減小Cp曲線更光滑，且總體曲線更接近文獻(xiàn)[11]貼體網(wǎng)格的測試結(jié)果；無論網(wǎng)格2 還是網(wǎng)格3，再次驗(yàn)證了同一網(wǎng)格下，計(jì)算階數(shù)越高時(shí)精度越高，模擬計(jì)算得到的結(jié)果與文獻(xiàn)結(jié)果吻合越好。

設(shè)置無量綱時(shí)間步長 Δt=1 × 10-3。在網(wǎng)格1 情況下，1 階DG 至3 階DG 計(jì)算得到的壓力系數(shù)Cp與文獻(xiàn)[11]中貼體網(wǎng)格模擬結(jié)果的對比見圖8。結(jié)果表明，本文方法計(jì)算所得壓力系數(shù)結(jié)果與文獻(xiàn)結(jié)果具有良好的一致性，且階數(shù)越高與文獻(xiàn)結(jié)果的吻合度越高。由圖6 和圖8 可以看出，同為網(wǎng)格1 且懲罰參數(shù)相同時(shí)，時(shí)間步長的取值會對計(jì)算結(jié)果產(chǎn)生影響。與無量綱時(shí)間步長為0.1 時(shí)相比，當(dāng)無量綱時(shí)間步長設(shè)置為與懲罰參數(shù)相同時(shí)（即1 × 10-3）），1 階DG 至3 階DG 得到的計(jì)算結(jié)果均更光滑且與文獻(xiàn)結(jié)果吻合更好。

本文方法懲罰項(xiàng)的添加基于Brinkman 懲罰法[25]，該方法中將固體視為滲透率極低的多孔介質(zhì)。Engels等[32]研究發(fā)現(xiàn)，當(dāng)使用顯式時(shí)間離散時(shí)，懲罰參數(shù)的選擇受CFL 的限制，顯式時(shí)間步長 Δt必須小于懲罰參數(shù)η；當(dāng)采用隱式時(shí)間離散時(shí)，允許懲罰參數(shù) η小于時(shí)間步長 Δt。Labert[33]對N-S 方程進(jìn)行了研究，揭示了懲罰參數(shù)與時(shí)間步長之間的密切耦合—誤差在懲罰參數(shù) η約等于時(shí)間步長 Δt時(shí)達(dá)到飽和。因此在體積懲罰法中通常建議 η=Δt，一種解釋為：懲罰項(xiàng)作為速度上的強(qiáng)阻尼，其引入了一個(gè)特征時(shí)間尺度，該特征時(shí)間尺度階數(shù)為η，必須用時(shí)間推進(jìn)來解決。

5.2 非定常圓柱繞流

為了驗(yàn)證對非定常解的有效性，測試二維圓柱周圍的低雷諾數(shù)渦脫落。當(dāng)Ma=0.2、Re=200時(shí)，計(jì)算域大小為x∈[-5D,10D]、y∈[-5D,5D]，其中D是圓柱的直徑，設(shè)置為1。在正方形區(qū)域x∈[-D,D]、y∈[-D,D]，使用均勻網(wǎng)格，均勻網(wǎng)格尺寸h選為0.05D（圖9）。

圖9 0.05D 非定常計(jì)算網(wǎng)格Fig.9 Computational grid for unsteady flow simulation(h=0.05D)

對不同階數(shù)下的圓柱繞流進(jìn)行數(shù)值模擬，1 階DG 和3 階DG 的馬赫數(shù)云圖分別如圖10 和圖11 所示。不論1 階DG 還是3 階DG 結(jié)果，從全流場的圖上均可以看到圓柱尾跡中有周期性的旋渦脫落。從物面附近流場圖可以看出，1 階DG 的物面附近存在階梯狀，3 階DG 的物面附近較1 階DG 的更為光滑。計(jì)算描述渦脫落的特征參數(shù)St，并與文獻(xiàn)[34-35]對比。1 階DG 的St=0.193，3 階DG 的St=0.20，與文獻(xiàn)[34]中St=0.194 以及文獻(xiàn)[35]中St=0.194較為吻合。本文所采用的體積懲罰法對邊界條件的實(shí)現(xiàn)是間接的，采用的是一種弱約束方式，該方法對邊界條件的實(shí)現(xiàn)精度與貼體網(wǎng)格相比較低，邊界附近的計(jì)算精度依賴于網(wǎng)格密度、數(shù)值方法的精度以及懲罰參數(shù)η的設(shè)置[32-33]。

圖10 Re = 200 時(shí)圓柱馬赫數(shù)云圖（p1）Fig.10 Mach number contours at Re=200 (p1)

圖12 為中軸線y=0上時(shí)均流向速度分布。由圖可看出，1 階DG 和3 階DG 的模擬結(jié)果中，圓柱內(nèi)部的速度為0，同時(shí)沿來流方向在圓柱后方近壁面流速為負(fù)值，存在先減后增的趨勢，故圓柱后方存在回流區(qū)，在遠(yuǎn)處流速趨于穩(wěn)定。整個(gè)回流區(qū)及近尾區(qū)的流速分布與文獻(xiàn)[34]結(jié)果吻合良好，但在恢復(fù)遠(yuǎn)區(qū)，可能由于遠(yuǎn)處網(wǎng)格尺度較大導(dǎo)致本文計(jì)算的結(jié)果偏小。總體上，3 階DG 的結(jié)果相較于1 階DG 更接近文獻(xiàn)[34]的結(jié)果。

圖12 沿中軸線y=0的時(shí)均流向速度分布Fig.12 Mean streamwise velocity in the y=0 direction

圖13 至圖15 分別給出了圓柱后的三個(gè)不同剖面x/D=1.0、2.0、4.5 上的時(shí)均流向速度分布。圖中實(shí)線為文獻(xiàn)[34]計(jì)算結(jié)果，實(shí)心三角和空心方形分別表示1 階DG 和3 階DG 的計(jì)算結(jié)果。可以看出，1 階DG 和3 階DG 的計(jì)算結(jié)果與文獻(xiàn)[34]結(jié)果在不同剖面處流向速度隨坐標(biāo)的變化規(guī)律一致，均關(guān)于圓柱中心線（y/D=0）對稱，且在兩側(cè)趨于穩(wěn)定數(shù)值。通過三個(gè)剖面處的結(jié)果對比，可以看出，同一網(wǎng)格下3 階DG 的計(jì)算結(jié)果相較于1 階DG 的計(jì)算結(jié)果與文獻(xiàn)更為吻合?？傮w來說，本文方法成功模擬了低雷諾數(shù)下非定常問題，且階數(shù)越高模擬的結(jié)果精度越高。

圖13 x/D=1.0截面的時(shí)均流速分布Fig.13 Mean streamwise velocity atx/D=1.0

圖14 x/D=2.0截面的時(shí)均流速分布Fig.14 Mean streamwise velocity atx/D=2.0

圖15 x/D=4.5截面的時(shí)均流速分布Fig.15 Mean streamwise velocity atx/D=4.5

6 結(jié)論

本文提出了適用于可壓縮流動數(shù)值模擬的結(jié)合高階DG 的浸沒邊界法，用體積懲罰方法來實(shí)現(xiàn)物面邊界條件，避免了切割網(wǎng)格。相較于顯式迭代，對懲罰項(xiàng)采用隱式時(shí)間離散允許較小的懲罰參數(shù)。采用IDW-IP 方法對物面附近的數(shù)據(jù)進(jìn)行重建，以獲得物面表面壓力。研究得出以下結(jié)論：

1）通過不同的流動模擬，包括靜態(tài)圓柱定常和非定常的流動，驗(yàn)證了本文方法計(jì)算可壓縮流動的有效性和準(zhǔn)確性。

2）高階DG 的高精度優(yōu)勢在本文的浸沒邊界方法中得到了顯著體現(xiàn)。隨著計(jì)算階數(shù)的提高，邊界條件的實(shí)現(xiàn)精度以及整體的數(shù)值精度都得到了顯著提高。

3）與貼體網(wǎng)格方法相比，IBM 最大的優(yōu)勢在于網(wǎng)格生成十分簡單，缺點(diǎn)在于邊界條件的實(shí)現(xiàn)是采用弱約束方式，因此邊界處理的精度不如貼體網(wǎng)格方法，計(jì)算結(jié)果依賴于網(wǎng)格密度、高階DG 的階數(shù)以及懲罰參數(shù) η的設(shè)置。

4）本文方法為未來復(fù)雜外形及移動邊界等要求高質(zhì)量貼體網(wǎng)格的計(jì)算提供了更為簡便的思路。但對移動邊界的處理還需進(jìn)一步研究驗(yàn)證。