基于有限體積格式的自適應(yīng)笛卡爾網(wǎng)格虛擬單元方法及其應(yīng)用

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

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

0 引言

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

對(duì)于復(fù)雜幾何外形的物體，使用笛卡爾網(wǎng)格方法，在物面邊界處網(wǎng)格單元必然與物面相交，從而產(chǎn)生不同形狀的切割網(wǎng)格單元。切割網(wǎng)格單元存在的一個(gè)普遍問題就是在切割的過程中會(huì)產(chǎn)生非常小的切割單元，這就導(dǎo)致了方程系統(tǒng)的剛性以及在物面處會(huì)產(chǎn)生解的非物理振蕩。并且對(duì)于時(shí)間依賴問題的數(shù)值模擬，會(huì)限制時(shí)間步長(zhǎng)且影響穩(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ù)雜物面邊界是一個(gè)新的發(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ì)算存儲(chǔ)量，但計(jì)算效率更高。同時(shí)，由于采用了有限體積方法而不是通常的有限差分方法，能更好地滿足守恒性要求。此外，本文利用基于速度的散度和旋度的判據(jù)，對(duì)流場(chǎng)進(jìn)行解的自適應(yīng)，能夠更為準(zhǔn)確地捕捉激波等流動(dòng)特征，數(shù)值結(jié)果顯示本文所使用的自適應(yīng)方法具有很好的分辨率，且易于實(shí)現(xiàn)。

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

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

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

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

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

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

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

2 虛擬單元方法

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

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

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

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

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

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

3 解自適應(yīng)

為了能夠自動(dòng)地捕捉流場(chǎng)特征，有必要針對(duì)流場(chǎng)解

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

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

4.1 NACA0012翼型的守恒性驗(yàn)證

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

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

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

4.3 雙錯(cuò)位NACA0012翼型

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

5 結(jié)論

本文發(fā)展了基于有限體積格式的自適應(yīng)笛卡爾網(wǎng)格虛擬單元方法，并成功地將這種方法應(yīng)用于不同流動(dòng)問題，并使用本文的方法成功地進(jìn)行多物體流動(dòng)問題的求解，得到了理想的結(jié)果。數(shù)值實(shí)驗(yàn)顯示本文方法是十分有效的，為將這些方法推廣并應(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.