胡 偶，趙 寧，劉劍明，王東紅

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

0 引言

眾所周知，計算流體力學(xué)中一個持續(xù)的障礙是復(fù)雜幾何外形的網(wǎng)格生成。自適應(yīng)笛卡爾網(wǎng)格作為一種新興的網(wǎng)格生成方法具有結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格方法所無法比擬的優(yōu)勢。自適應(yīng)笛卡爾網(wǎng)格生成簡單，省時，同時網(wǎng)格結(jié)構(gòu)簡單，易于實現(xiàn)自動化、具有更強的自適應(yīng)能力，可提高計算精度。尤其對于復(fù)雜幾何外形問題，如果使用傳統(tǒng)的結(jié)構(gòu)網(wǎng)格，通常使用嵌套網(wǎng)格和多塊對接網(wǎng)格等，雖然這些方法已經(jīng)得到了成功的使用，但它們需要交換不同網(wǎng)格間的數(shù)據(jù)信息，會遇到一些復(fù)雜的插值和幾何的問題;非結(jié)構(gòu)網(wǎng)格雖然易于生成復(fù)雜外形的網(wǎng)格，但非結(jié)構(gòu)網(wǎng)格生成費時，計算時間存儲量一般都比結(jié)構(gòu)網(wǎng)格高;而使用自適應(yīng)笛卡爾網(wǎng)格則可以避免這些問題，使得網(wǎng)格生成變成一個簡單而省時的過程。

對于復(fù)雜幾何外形的物體，使用笛卡爾網(wǎng)格方法，在物面邊界處網(wǎng)格單元必然與物面相交，從而產(chǎn)生不同形狀的切割網(wǎng)格單元。切割網(wǎng)格單元存在的一個普遍問題就是在切割的過程中會產(chǎn)生非常小的切割單元，這就導(dǎo)致了方程系統(tǒng)的剛性以及在物面處會產(chǎn)生解的非物理振蕩。并且對于時間依賴問題的數(shù)值模擬，會限制時間步長且影響穩(wěn)定性。為了克服上述問題，很多學(xué)者給出了多種處理方法，主要有混合網(wǎng)格方法(hybrid grid)［1］，融合單元方法(merged cell approach)［2］，嵌入單元方法(embedded cell method)［3］，以及虛擬單元方法(ghost-cell method)［4］或浸入邊界方法(immersed boundary method)［5］等。

自適應(yīng)笛卡爾網(wǎng)格與虛擬單元方法結(jié)合處理復(fù)雜物面邊界是一個新的發(fā)展方向。最近，Dadone和Grossman 在他們的一系列文章［6，7，8］中，系統(tǒng)研究了結(jié)構(gòu)網(wǎng)格中的虛擬單元方法，并把他們命名為CCST(Curvature Corrected Symmetry Technique)方法，而且這些有效的邊界處理方法也被他們推廣用于笛卡爾網(wǎng)格［9，10，11］，并取名為GBCM方法。本文基于有限體積格式，將虛擬單元方法與自適應(yīng)笛卡爾網(wǎng)格結(jié)合，采用了按邊掃描的有限體積方法。此種方法不同于通常的叉樹結(jié)構(gòu)的自適應(yīng)笛卡爾網(wǎng)格方法，雖然增加了計算存儲量，但計算效率更高。同時，由于采用了有限體積方法而不是通常的有限差分方法，能更好地滿足守恒性要求。此外，本文利用基于速度的散度和旋度的判據(jù)，對流場進行解的自適應(yīng)，能夠更為準(zhǔn)確地捕捉激波等流動特征，數(shù)值結(jié)果顯示本文所使用的自適應(yīng)方法具有很好的分辨率，且易于實現(xiàn)。

1 控制方程與數(shù)值方法

采用自適應(yīng)笛卡爾網(wǎng)格最突出的問題就是在不同層次的網(wǎng)格過渡區(qū)會出現(xiàn)懸掛節(jié)點情況，這對代碼的編制以及數(shù)值模擬的精度都有很大的影響。本文采用一種基于自適應(yīng)笛卡兒網(wǎng)格的有限體積方法，將這些帶有懸掛節(jié)點的網(wǎng)格邊當(dāng)作兩條邊來處理，從而有效地解決了懸掛節(jié)點問題。

二維Euler方程的積分守恒型的表達式如下:

式中W為守恒變量，F(xiàn)，G分別為笛卡爾坐標(biāo)系下的x，y方向的通量?；谟邢摅w積格式的流動求解器主要包括以下幾個方面的內(nèi)容:解的重構(gòu)，Riemann求解器，時間方向離散等。

對于有限體積格式，解的重構(gòu)主要有兩種方法:Green-Gauss方法和最小二乘方法。由于本文所采用的是自適應(yīng)笛卡爾網(wǎng)格，網(wǎng)格的結(jié)構(gòu)相對簡單，但又存在不同層次網(wǎng)格過渡的問題，過渡區(qū)網(wǎng)格結(jié)構(gòu)不一致，同時本文采用的叉樹網(wǎng)格數(shù)據(jù)結(jié)構(gòu)能夠十分便捷地給出網(wǎng)格單元的鄰居的信息，因此本文采用了最小二乘法進行解的重構(gòu)，具體表達形式參見文獻［12，13］。

本文采用的Riemann求解器是由Liou和Steffen提出來AUSM格式［14］，AUSM格式具有結(jié)構(gòu)簡單，對激波具有很好的捕捉等特性。時間方向離散采用四步Runge-Kutta顯示時間推進方法。

同時，針對本文所采用的自適應(yīng)笛卡爾網(wǎng)格，由于網(wǎng)格的非貼體性，需要采用特殊的方法來處理物體壁面邊界，給出合適的壁面邊界條件，本文采用一種稱作虛擬單元方法的浸入邊界處理方法。

2 虛擬單元方法

采用笛卡爾網(wǎng)格，物體與網(wǎng)格相交，如圖1所示。圖中圓圈和三角形表示物體內(nèi)部的網(wǎng)格點，即虛擬單元(Ghost Cells);實心正方形表示虛擬單元點關(guān)于物面邊界的對稱點。本文使用虛擬單元方法需要給出浸入邊界(Immersed Boundary)內(nèi)部網(wǎng)格點處的物理量，例如點A處的值。

在文獻［6］中，Dadone等人針對貼體結(jié)構(gòu)網(wǎng)格提出了一種新的邊界條件處理方法，并命名為CCST方法(Curvature-Corrected Symmetry Technique)。虛擬單元點上的流場值是通過在物面附近假設(shè)的一個法向具有局部對稱分布的熵(S)與總焓(H)的渦流場得到的。這種流動模型滿足如下的法向動量方程:

圖1 笛卡爾網(wǎng)格Fig.1 Cartesian grid

其中R是帶符號的局部曲率半徑，如果曲率中心在物體內(nèi)部為正，反之為負;是切向速度，滿足無穿透邊界條件·=0。由于沿著物體的表面法向強加反對稱的熵與總焓的法向?qū)?shù) ?S/?n，?H/?n。這種熵與總焓分布使得在無旋流動時，熵與總焓的法向?qū)?shù)為零，甚至在物面出現(xiàn)渦時也能滿足Crocco關(guān)系。對于二維問題，為了求虛擬單元點A的流場值，給出了如下的表達式:

上式中B點為虛擬單元點A關(guān)于物面的法向?qū)ΨQ點。Dadone在文章中比較了GBCM方法與貼體網(wǎng)格下的CCST方法，得出GBCM方法與貼體網(wǎng)格下的CCST方法具有相同的優(yōu)越性。

在數(shù)值計算中虛擬單元中心A所對應(yīng)的法向?qū)ΨQ點B變量的值可以通過一般插值得到，例如圖1中的B點，我們可以通過其周圍的四個流場網(wǎng)格點C、D、E、F采用雙線性插值［5，6］方法(三維使用三線性插值［7］)得到，而對于那種不能找到四個流場網(wǎng)格點的情況，我們可以采用逆距離插值方法，這樣在我們的代碼中就能保證始終使用流體內(nèi)部網(wǎng)格點來插值。

3 解自適應(yīng)

為了能夠自動地捕捉流場特征，有必要針對流場解

4 數(shù)值驗證與應(yīng)用

為驗證本文基于有限體積格式的自適應(yīng)笛卡爾網(wǎng)格虛擬單元方法的有效性，本文分別以NACA0012翼型和RAE2822翼型繞流問題為例，進行了數(shù)值模擬，最后采用本文的方法對雙錯位NACA0012翼型進行了數(shù)值模擬。

4.1 NACA0012翼型的守恒性驗證

本文以NACA0012翼型繞流為例來比較基于有限體積格式的自適應(yīng)笛卡爾網(wǎng)格方法相對于傳統(tǒng)的基于有限差分方法的自適應(yīng)笛卡爾網(wǎng)格的差別。來流馬赫數(shù) M∞=0.8，攻角 α =0.0°，計算網(wǎng)格如圖 2(a)所示，網(wǎng)格進行了7次局部加密，在計算的過程中網(wǎng)格不進行解自適應(yīng)，以避免網(wǎng)格的不同影響計算結(jié)果的比較。翼型的表面壓強系數(shù)分布如圖2(b)所示，從圖中我們可以看出，使用傳統(tǒng)的有限差分方法激波位置明顯偏離正確的位置，所以為了克服守恒性問題，Berger［16］等人提出需要在粗細網(wǎng)格界面處對粗網(wǎng)格單元作守恒修正，使得流入到細網(wǎng)格的通量等于從粗網(wǎng)格流出的通量，但這個過程增加了計算復(fù)雜度，而使用本文的有限體積方法，守恒性能很好地滿足，激波位置與參考文獻［17］的結(jié)果吻合的很好。

4.2 RAE2822翼型解自適應(yīng)驗證

考慮來流馬赫數(shù) M∞=0.75及攻角 α=3.19°的RAE2822翼型繞流。在這種狀態(tài)下，翼型的上表面會形成一道強激波。在計算的過程中我們進行了三次解自適應(yīng)加密，三次解自適應(yīng)粗化，最終得到的計算網(wǎng)格如圖3(a)所示，計算網(wǎng)格的總數(shù)為10320個。圖3(b)和圖3(c)分別給出了壓強云圖和翼型表面的壓強系數(shù)分布。從圖中可以看出，壓力值在網(wǎng)格變化區(qū)域過渡光滑，激波位置捕捉精確，與參考文獻［18］中的結(jié)果吻合的比較好。

4.3 雙錯位NACA0012翼型

以雙錯位NACA0012翼型為例，驗證本文方法對多物體計算的有效性。其中來流馬赫數(shù)M∞=0.70，攻角α=0.0°。在該算例中，兩個翼型之間會形成一道強激波，類似于管道流動。在計算的過程中本文進行了三次解自適應(yīng)，最終的計算網(wǎng)格圖如圖4(a)所示，網(wǎng)格總數(shù)為16337個。圖4(b)和4(c)分別給出了壓強云圖和上下兩個翼型表面的壓強系數(shù)分布，從圖中可以看出，激波結(jié)構(gòu)清晰，位置準(zhǔn)確，與參考文獻［18］中的結(jié)果吻合的比較好。

5 結(jié)論

本文發(fā)展了基于有限體積格式的自適應(yīng)笛卡爾網(wǎng)格虛擬單元方法，并成功地將這種方法應(yīng)用于不同流動問題，并使用本文的方法成功地進行多物體流動問題的求解，得到了理想的結(jié)果。數(shù)值實驗顯示本文方法是十分有效的，為將這些方法推廣并應(yīng)用到三維問題奠定了良好的基礎(chǔ)。

［1］CHARLTON E F.An otree solution to conservation-laws over arbitrary regions(OSCAE)with applications to aircraft aerodynamics［D］.［Ph.D.Thesis］.University of Michigan，1997.

［2］UDAYKUMAR H S.Multiphase dynamics in arbitrary geometries on fixed cartesian grids［J］.Journal of Computational Physics，1997，137:366 －405.

［3］MARSHALL D D.Extending the functionalities of cartesian grid solvers:viscous effects modeling and MPI parallelization［D］.Ph.D.Thesis，Georgia Institute of Technology，2002.

［4］DADONE A and GROSSMAN B.An immersed body methodology for inviscid flows on cartesian grids［R］.AIAA Paper，2002 －1059，2002.

［5］FORRER H，JELTSCH R.A higher order boundary treatment for cartesian-grid method［J］.Journal of Computational Physics，1998，140(2):259 －277.

［6］DADONE A，GROSSMAN B.Surface boundary conditions for the numerical solution of the Euler equations［R］.AIAA 93－3334.1993.

［7］DADONE A.Symmetry technique for the numerical solution of the 2d Euler equations at impermeable boundaries［J］.Int.J.Numer.Meth.Fluids.，1998，28(7):1093 －1108.

［8］DADONE A，GROSSMAN B.Surface boundary conditions for the numerical solution of the Euler equations in three dimensions［J］.Lect.Notes.Phys.，1995，453:188 －94.

［9］DADONE A，GROSSMAN B.Ghost-cell method for inviscid two-dimensional flows on cartesian grids［J］.AIAA JOURNAL，2004，42(12):2499 －2507.

［10］DADONE A，GROSSMAN B.Ghost-cell method with far field coarsening and mesh adaptation for Cartesian grids［J］.Computers＆ Fluids 2006，35(7):676－687.

［11］DADONE A，GROSSMAN B.Ghost-cell method for inviscid three-dimensional flows on cartesian grids［J］.Computers ＆Fluids，2007，36(10):1513 －1528.

［12］WANG Z J.A quadtree-based adaptive Cartesian/quad grid flow solver for Navier-Stokes equations［J］.Computers ＆Fluids，1998，27(4):529 －549.

［13］VENKATAKRISHNAN V.On the accuracy of limiters and convergence to steady state solutions［R］.AIAA Paper 93 －0880，1993.

［14］LIOU M S，STEFFEN C J.A new flux splitting scheme［J］.J.Computional Physics，1993，107:23 －39.

［15］DARREN L.DE ZEEUW.A quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations［D］.［PhD thesis］.The University of Michigan，1993.

［16］BERGER M J.Local adaptive mesh refinement for shock hydrodynamics［J］.J.Comput.Phys.，1989，82:64 －84.

［17］STOLCIS L，JOHNSTON L J.Solution of the Euler equations on unstructured grids for two-dimensional compressible flow［J］.Aeronautical Journal，1990，94:181 －195.

［18］馬志華.自適應(yīng)無網(wǎng)格及網(wǎng)格和無網(wǎng)格混合算法研究［D］.［博士學(xué)位論文］.南京:南京航空航天大學(xué)，2008.