盧鳳翎， 陳小前， 禹彩輝， 苗萌

1.國防科技大學(xué) 航天科學(xué)與工程學(xué)院， 長沙 410073 2.空間物理重點實驗室， 北京 100076

一種基于求解橢圓型方程的結(jié)構(gòu)動網(wǎng)格生成方法

盧鳳翎1,2,*， 陳小前1， 禹彩輝2， 苗萌2

1.國防科技大學(xué) 航天科學(xué)與工程學(xué)院， 長沙 410073 2.空間物理重點實驗室， 北京 100076

基于橢圓型網(wǎng)格生成法，實現(xiàn)了一種簡單高效的貼體結(jié)構(gòu)動網(wǎng)格生成方法，可用于具有移動邊界問題的非定常流動數(shù)值模擬。該方法提出，在網(wǎng)格變形過程中，Poisson方程需要的控制網(wǎng)格間距和正交性的源項可以通過提取已知的靜態(tài)網(wǎng)格源項直接得到，并在整個動網(wǎng)格生成過程中保持不變。因此，在橢圓型網(wǎng)格生成中需要通過外迭代確定源項的過程可以得到省略，而且該方法不需要人工指定參數(shù)。這使得方法具有高效和易于嵌入到已有程序中的特點。數(shù)值模擬結(jié)果證明，采用這種方法獲得的網(wǎng)格能夠較好地保持靜態(tài)網(wǎng)格原有的正交性和光滑性，在相同迭代步數(shù)約束下，網(wǎng)格求解效率低于傳統(tǒng)彈簧模擬法,但魯棒性優(yōu)于彈簧模擬法。

動網(wǎng)格生成； 橢圓型網(wǎng)格生成； 計算流體力學(xué)； 數(shù)值網(wǎng)格生成； 網(wǎng)格變形方法

工程應(yīng)用中經(jīng)常碰到具有移動邊界的非定常流體動力學(xué)問題，例如火箭級間分離、飛機(jī)機(jī)彈分離、飛機(jī)俯仰運動和渦輪機(jī)內(nèi)流動等。要處理這些問題，一般需要涉及到動網(wǎng)格生成技術(shù)。也就是說，網(wǎng)格需要跟隨物體運動或產(chǎn)生變形以適應(yīng)物體的運動。

對某些問題，例如進(jìn)行飛機(jī)俯仰運動的非定常模擬，可以考慮采用網(wǎng)格不變形但隨物體同時移動的所謂剛性網(wǎng)格方法去解決；但在處理具有相對運動部件的問題時，例如多段翼中的襟翼運動模擬以及火箭級間分離等問題的模擬，就不能采用剛性網(wǎng)格方法，而需要考慮其他的方法，例如重疊網(wǎng)格 (Overlapping Grid) ，也稱為Chimera 網(wǎng)格或者嵌套 (Overset) 網(wǎng)格的方法[1-4]。這種方法允許網(wǎng)格間重疊、嵌套。當(dāng)物體運動時，貼體的部件網(wǎng)格隨之運動，數(shù)值計算在各個分塊子網(wǎng)格上分別進(jìn)行，并在重疊或嵌套的區(qū)域間通過插值來實現(xiàn)流場信息的傳遞。重疊網(wǎng)格雖然降低了網(wǎng)格生成難度，但是挖洞、找點使得網(wǎng)格構(gòu)造變得極其耗時，這種方法的計算效率成為一個非常突出的問題。

此外,在求解上述問題時,還可以考慮采用動網(wǎng)格技術(shù),即在每一時間步對網(wǎng)格進(jìn)行部分或整體重新生成。網(wǎng)格重新生成時，如果網(wǎng)格規(guī)模和拓?fù)浣Y(jié)構(gòu)發(fā)生了變化，則稱之為網(wǎng)格重構(gòu)法；如果各個時間步上的網(wǎng)格與原始網(wǎng)格在規(guī)模和拓?fù)浣Y(jié)構(gòu)上均保持一致，則稱之為網(wǎng)格變形法。按此分類，網(wǎng)格變形法只是將物面的剛性運動或者局部變形傳遞到空間網(wǎng)格，對空間網(wǎng)格的坐標(biāo)進(jìn)行適當(dāng)更新，而不改變網(wǎng)格的規(guī)模和拓?fù)浣Y(jié)構(gòu)，效率更高，也更容易與流場求解相結(jié)合，因而更受歡迎，也是動網(wǎng)格技術(shù)研究的重要內(nèi)容。本文以結(jié)構(gòu)網(wǎng)格變形技術(shù)為研究對象，但只涉及物面作剛性運動情況下的網(wǎng)格變形方法。

實際上，有多種方法可以實施網(wǎng)格變形[5]，例如超限插值(Trans-Finite Interpolation，TFI)法、彈簧模擬法、彈性體法、Delaunay圖映射法和徑向基函數(shù)(Radial Basis Functions, RBF)法等。超限插值[6-7]法屬于代數(shù)插值類方法，一般只用于結(jié)構(gòu)網(wǎng)格，該方法計算量小，可以實現(xiàn)相對復(fù)雜的網(wǎng)格變形，但難以保證變形后的網(wǎng)格質(zhì)量，特別是物面網(wǎng)格的正交性。而其他幾種方法既可用于結(jié)構(gòu)網(wǎng)格，也可用于非結(jié)構(gòu)網(wǎng)格。其中，彈簧模擬法最初是由Batina[8]在1990年提出的，他把這種方法應(yīng)用于求解一個受到強(qiáng)迫振動的翼型。在其原始形式中，單元的每一邊被假設(shè)為一個具有剛度的彈簧，并且剛度與其邊長成反比。物體運動引起的網(wǎng)格變形反映在網(wǎng)格節(jié)點的位移上，此位移可以通過求解基于靜平衡狀態(tài)給出的一個大型非線性系統(tǒng)方程獲得。把節(jié)點位移與原節(jié)點位置相加即可得到變形后的網(wǎng)格。原始的彈簧模擬法魯棒性較差，為此，文獻(xiàn)[9-12]提出增加抗扭彈簧、張力彈簧等改進(jìn)措施，取得了一些好的結(jié)果。彈簧模擬法的不足之處在于，這種方法不能對網(wǎng)格大小、形狀進(jìn)行控制，在位移較大的地方容易出現(xiàn)網(wǎng)格交錯，產(chǎn)生負(fù)體積。對于結(jié)構(gòu)網(wǎng)格，這種方法不能提供可靠的、自動保證網(wǎng)格正交性和光滑性的手段。彈性體法[13]是將計算域比擬為彈性體，通過求解彈性力學(xué)平衡方程得到新的網(wǎng)格坐標(biāo)，這種方法具有較強(qiáng)的網(wǎng)格變形能力，但求解大型方程組帶來較大計算量，限制了它的使用。Delaunay圖映射法[14]是一種簡單、高效、無迭代過程的網(wǎng)格變形技術(shù)，它是通過在物面與外邊界之間生成Delaunay圖，也稱為背景網(wǎng)格的方式，將變形量映射到Delaunay背景網(wǎng)格上，然后再通過背景網(wǎng)格插值得到新網(wǎng)格坐標(biāo)。該方法計算量小，在凸邊界小變形問題中具有較大優(yōu)勢，但對于凹邊界網(wǎng)格變形處理比較復(fù)雜，而且當(dāng)背景網(wǎng)格出現(xiàn)交叉時，會導(dǎo)致此后的變形網(wǎng)格質(zhì)量越來越差。RBF法[15-16]是近期變形網(wǎng)格方法研究的一個熱點, 其原理是運用徑向基函數(shù)對物面網(wǎng)格點的位移進(jìn)行擬合，并以此來建立徑向基函數(shù)插值模型，然后利用該模型計算空間網(wǎng)格節(jié)點變形量并更新網(wǎng)格坐標(biāo)，此方法的優(yōu)點是能夠適應(yīng)復(fù)雜外形的大變形問題，網(wǎng)格質(zhì)量高。但在應(yīng)用于大規(guī)模網(wǎng)格變形時，如果直接采用所有物面節(jié)點來構(gòu)造徑向基插值函數(shù)的話，也會存在計算量過大的問題。為此，一些作者提出采用數(shù)據(jù)精簡的辦法，例如，利用貪婪算法減少徑向基基函數(shù)的方法[16-17]、峰值選點法[18]、RBF與Delaunay圖映射法相結(jié)合的方法[19]等，以提高網(wǎng)格變形效率，取得了好的效果。

可以看到, 上述各種網(wǎng)格變形法幾乎都是以代數(shù)或幾何的方法去解決動網(wǎng)格生成問題，而不像在靜態(tài)結(jié)構(gòu)網(wǎng)格生成時那樣，多以求解偏微分方程的方法來解決網(wǎng)格生成問題。這是否意味著偏微分方程法不適用于動網(wǎng)格生成呢，當(dāng)然不是。事實上，早有研究人員嘗試使用這種方法來解決網(wǎng)格變形問題，只是還沒有獲得實用的結(jié)果。例如，Johnson和Tezduyar[20]曾嘗試通過求解由網(wǎng)格速度導(dǎo)出的線彈性泊松方程，以產(chǎn)生隨流體邊界或界面運動的網(wǎng)格。Ushijima[21]則嘗試采用橢圓型網(wǎng)格生成法在每一時間步生成貼體網(wǎng)格，這種方法即便在較復(fù)雜的物體運動下也能生成高質(zhì)量的網(wǎng)格。但在每一時間步求解泊松方程的計算代價實在太高，遠(yuǎn)遠(yuǎn)高于彈簧模擬法。為了減少求解泊松方程的計算開銷，Grüber 和Cartens[22]采用了一種簡化策略，不是在每一時間步，而是只在選定的有限步求解泊松方程。比如只在物體運動的最大和最小時間步兩個位置給出由偏微分方程法生成的網(wǎng)格，其他中間位置的網(wǎng)格則用這兩個網(wǎng)格進(jìn)行插值。嚴(yán)格的講，這種方法只能算是偏微分方程法與代數(shù)插值法的組合。

本文提出了一種基于求解橢圓型微分方程的網(wǎng)格變形法，可用于結(jié)構(gòu)動網(wǎng)格的生成。數(shù)值算例表明，用該方法生成的動網(wǎng)格可以保持原始網(wǎng)格的光滑性、正交性，具有良好的魯棒性和較高的計算效率。

1 一般橢圓型網(wǎng)格生成法

橢圓型網(wǎng)格生成法是靜態(tài)結(jié)構(gòu)網(wǎng)格生成中最受歡迎、使用最普遍的一種方法。這種方法可以生成高質(zhì)量的網(wǎng)格，具有網(wǎng)格正交性可控、光滑的C2連續(xù)性、魯棒性較好的優(yōu)點。

用笛卡兒坐標(biāo)r=[xyz]T代表物理空間Ω內(nèi)的點，曲線坐標(biāo)ρ=[ξηζ]T代表計算空間D內(nèi)的點。計算空間內(nèi)的點ρ與物理空間內(nèi)的點r通過可逆映射T相聯(lián)系：

T∶r→ρ?T-1∶ρ→r

因此，網(wǎng)格生成其實就是一種由式(1)向量值函數(shù)確定的映射：

(1)

映射函數(shù)被要求滿足泊松方程：

(2)

式中：P、Q、R為強(qiáng)迫函數(shù)，也稱為源項。

通過交換自變量(x，y，z)與因變量(ξ，η，ζ)的位置，在變換后的計算域D內(nèi)，式(1)可寫為

α1rξ ξ+α2rη η+α3rζζ+2(β12rξη+β23rη ζ+

β31rζ ξ)+J2(rξP+rηQ+rζR)=0

(3)

式中：

(4)

另有，物理空間邊界?D上的邊界條件為

(5)

式(3)經(jīng)整理后可寫為

α1(rξ ξ+φPrξ)+α2(rηη+φQrη)+α3(rζζ+φRrζ)+

2(β12rξη+β23rη ζ+β31rζ ξ)=0

(6)

數(shù)值求解式(6)，方程的解即為笛卡兒坐標(biāo)下的網(wǎng)格點。內(nèi)點單元的大小和形狀通過調(diào)節(jié)泊松方程中的源項來實現(xiàn)。顯然，源項在生成適宜精確模擬流動問題的網(wǎng)格時起到了關(guān)鍵性的作用。

源項控制方法有多種形式，一般可分成2種類型：第1種類型是由式(6)導(dǎo)出源項表達(dá)式，然后在給定的坐標(biāo)邊界上以正交條件和指定的第一層網(wǎng)格與壁面的距離作為約束條件求出初始邊界源項，初始內(nèi)點源項則通過插值獲得，然后迭代求解方程；第2種類型是在迭代過程中根據(jù)源項的變化情況，采用“人工”控制源項值，來得到所期望的方程。Thomas和Middlecoff[23]以及Sorenson 和 Steger[24-25]方法屬于第一類，Hilgenstock[26]方法屬于第二類。這些方法都要求將物面和外邊界處的源項值內(nèi)插至內(nèi)部網(wǎng)格點處。

求解式(6)最實用的方法是高斯-賽德爾連續(xù)超松弛(Gauss-Seidel Successive Overrelaxation，GSSOR)迭代法。迭代循環(huán)一般由外迭代和內(nèi)迭代組成，在外迭代循環(huán)中改變源項，在內(nèi)迭代循環(huán)中源項保持不變。

另外，在結(jié)構(gòu)網(wǎng)格生成中，還需要用到以下幾個概念。譬如，一對一的映射也被稱為規(guī)則映射。而采用規(guī)則映射得到的網(wǎng)格往往被稱為規(guī)則網(wǎng)格。

從更實際的角度出發(fā)，規(guī)則網(wǎng)格也可以采用以下的定義[27]：如果一個網(wǎng)格的所有網(wǎng)格線均處于一個物理域的規(guī)定域內(nèi)，并且同族網(wǎng)格線不相交，不同族的兩條網(wǎng)格線如果相交，且只相交一次，則稱這樣的網(wǎng)格是規(guī)則的。

對于由式(6)生成的網(wǎng)格，可以做出以下論斷[27]：假設(shè)已由式(6)生成了一個規(guī)則網(wǎng)格，那么以此網(wǎng)格為初始網(wǎng)格, 改變源項值，但保持新源項的符號分布與初始網(wǎng)格相同，則所生成的新網(wǎng)格也是規(guī)則的。

2 本文動網(wǎng)格生成算法

假設(shè)已有一個網(wǎng)格，準(zhǔn)備用于帶有移動邊界的某非定常流動問題的模擬。該網(wǎng)格可以是由任何網(wǎng)格生成方法獲得的網(wǎng)格，例如，該網(wǎng)格可以是由代數(shù)方法或者雙曲方法生成，也可以直接由橢圓型網(wǎng)格生成法生成。以此網(wǎng)格為出發(fā)點構(gòu)造后續(xù)每一時間步的動網(wǎng)格。

2.1 源項的獲取

顯然，無論采用何種方法生成的網(wǎng)格，它首先應(yīng)該是規(guī)則的。否則，它將不適宜被用來模擬一個實際的流動問題。剩下的問題是，是否可以在泊松類方程的解空間中找到與此相同的一個網(wǎng)格。對這個問題先不做正面回答。觀察式(6)的網(wǎng)格生成系統(tǒng)，可以看到其中所有導(dǎo)數(shù)項以及度量系數(shù)都可以由已知網(wǎng)格通過差分方法計算得到，系統(tǒng)中未知的是源項φP、φQ和φR。式(6)可以寫為

α1rξφP+α2rηφQ+α3rζφR=-α1rξ ξ-

α2rηη-α3rζζ-2(β12rξη+β23rη ζ+β31rζξ)

(7)

如果令

b=-α1rξξ-α2rηη-α3rζζ-2(β12rξη+β23rη ζ+β31rζξ)

式(7)又可寫成標(biāo)量形式：

(8)

如果用i、j、k分別表示相應(yīng)于坐標(biāo)ξ、η、ζ方向的整型指標(biāo)，式(7)和式(8)中的導(dǎo)數(shù)項可以表示成中心差分的形式：

(9)

其他2個坐標(biāo)方向η和ζ的導(dǎo)數(shù)也有類似的表達(dá)。

式(8)的邊界條件可以十分方便地采用有限體積法中的啞元(Dummy Cell)技術(shù)，或者稱為虛擬單元技術(shù)來確定。簡單說，就是通過在每個網(wǎng)格塊的邊界處附加另外的虛擬網(wǎng)格來實現(xiàn),以i=Imax邊界為例,第一層虛擬單元為

rImax+1,j,k=rImax,j,k+rImax,j,k-rImax-1,j,k

第二層虛擬單元可在第一層虛擬單元的基礎(chǔ)上生成。當(dāng)然，如果邊界導(dǎo)數(shù)項采用單邊差分算法求解的話，甚至可以不用虛擬單元，但是會造成編程上麻煩。

通過求解式(8)可以在初始網(wǎng)格的每一點上都獲得一個源項φP、φQ和φR。同時，也可以看到，對于任意一個規(guī)則網(wǎng)格，有且只有唯一一個與式(8)的解相對應(yīng)的源項。因此，可以作出以下論斷：每一個規(guī)則網(wǎng)格都可視為是由具有某種源項分布的泊松類方程生成的結(jié)果。

顯然，源項分布是橢圓型網(wǎng)格生成法中的一個關(guān)鍵因素。它包含了描述一個網(wǎng)格的重要信息，但又不是網(wǎng)格本身。觀察物理域內(nèi)的泊松方程式(2)，可以看到源項滿足的是曲線坐標(biāo)對笛卡兒坐標(biāo)二階導(dǎo)數(shù)項組成的方程，因此，源項的變化是通過影響這些導(dǎo)數(shù)項來改變網(wǎng)格位置的，不是直接對網(wǎng)格位置進(jìn)行控制。網(wǎng)格本身也就是式(2)的離散解只有在每個坐標(biāo)變換方向的離散點數(shù)、正交性要求、網(wǎng)格間距要求(可轉(zhuǎn)化為對源項分布的要求)確定后，才能被確定。

通過上面的討論，知道一個規(guī)則網(wǎng)格必有一個屬于這個網(wǎng)格本身的源項分布?；蛘哒f，對任意一個規(guī)則網(wǎng)格總存在滿足式(8)的一個源項分布。本文的重點在于如何在動網(wǎng)格生成中利用初始網(wǎng)格中已知的源項分布。

對某些移動邊界問題，如飛機(jī)外掛物分離，雖然存在相對運動，但物面邊界的外形是不變的。因此，如果采用一個物理空間內(nèi)的解析函數(shù)f(x,y,z)來定義這個物面邊界，并且假定它對應(yīng)計算空間內(nèi)的坐標(biāo)面ξ=ξmin，如圖1所示，盡管物面邊界經(jīng)過了位移和/或轉(zhuǎn)動，但向量值函數(shù)f(x,y,z)仍然保持不變。這意味著此物面邊界上的點在經(jīng)歷了從t0到t1的運動后，對曲線坐標(biāo)的一階和二階以及交叉導(dǎo)數(shù)也沒有變化。因此，有

(10)

(11)

式中：i，j=1，2，3；ξ1=ξ;ξ2=η;ξ3=ζ。

如前所述，在時間t=t0時刻的初始網(wǎng)格源項(或稱為靜態(tài)網(wǎng)格)可由式(8)確定。再考慮到式(10)和式(11)時，可知在ξ=ξmin邊界面上的源項是時不變的。因此可以得到

(12)

式中：L=P，Q，R。事實上，式(12)對于其他外形不變的靜止或運動的邊界也同樣成立。

除了物面邊界和遠(yuǎn)場邊界外，網(wǎng)格中一般還有內(nèi)部邊界，即塊與塊之間的對接面，或者是O形以及C形拓?fù)渲芯哂邢嗤奈锢砜臻g坐標(biāo)，但分屬不同的計算坐標(biāo)面而具有邊界含義的網(wǎng)格面(有的資料稱為坐標(biāo)切口, Coordinate Cut)。這種邊界存在于網(wǎng)格內(nèi)部，在新的網(wǎng)格生成過程中，不能證明式(12)的存在。但從源項控制的角度看，源項直接控制的是邊界附近網(wǎng)格的角度和網(wǎng)格間距，如果只要求內(nèi)部邊界附近的網(wǎng)格角度和網(wǎng)格間距得到保持，而不關(guān)心內(nèi)部邊界的實際位置的話，完全可以在內(nèi)部邊界上人為指定式(12)成立。

這意味著，所有用于動網(wǎng)格生成的邊界源項都可以通過繼承靜態(tài)網(wǎng)格的邊界源項而被決定下來。至于內(nèi)部源項，一般可以通過在同一時刻的兩兩邊界之間插值獲得。但在本文中，考慮到由式(8)獲得的解不但包括了源項本身，而且也隱含地包括了它們的插值關(guān)系。因此，不需要重新插值，而是直接使用由式(8)獲得的解作為動網(wǎng)格生成的源項。

2.2 動網(wǎng)格生成求解過程

為了更容易理解本文提出的動網(wǎng)格生成方法，有必要先回顧一下通常用于靜態(tài)橢圓型網(wǎng)格生成的數(shù)值求解過程。這個迭代求解過程可以歸納為

1) 假定邊界源項初值。

2) 改變邊界源項值，即實施“源項控制”，使其更接近于對邊界附近網(wǎng)格角度和間距控制所需要的源項值。

3) 將邊界上新的源項值插值到內(nèi)部網(wǎng)格點處。

4) 迭代求解出離散網(wǎng)格生成系統(tǒng)。

5) 重復(fù)步驟3)和4)，直到獲得收斂解。

稱步驟4)中的迭代過程為內(nèi)迭代，而步驟2)～4)的迭代為外迭代。源項只在外迭代中改變，在內(nèi)迭代中不變。

顯然，在靜態(tài)網(wǎng)格生成中，外迭代的主要功能是決定出源項值。在本文中，由于用于動網(wǎng)格生成的源項假定為從靜態(tài)網(wǎng)格繼承而來，因而是已知的。所以，新的動網(wǎng)格生成過程為

1) 讀入已知靜網(wǎng)格。

2) 通過求解式(8)，從已知靜網(wǎng)格提取源項。

3) 在移動邊界和其他邊界上施加與靜網(wǎng)格相同的網(wǎng)格點分布。

4) 采用上一時間步的網(wǎng)格坐標(biāo)值作為當(dāng)前時間步的初始解。

5) 迭代求解離散方程式(6)直至收斂。

在步驟5)的迭代收斂過程中，源項始終保持不變。因此不像靜態(tài)網(wǎng)格生成中那樣需要外迭代過程，這使得動態(tài)網(wǎng)格生成更加高效。

以上動網(wǎng)格生成過程，新的源項不僅繼承了靜態(tài)網(wǎng)格源項的符號和數(shù)值，也繼承了對網(wǎng)格特性的控制。因此邊界附近的網(wǎng)格正交性以及網(wǎng)格間距也得到了保持。

3 算法驗證

3.1 帶平移和轉(zhuǎn)動的NACA0012翼型動網(wǎng)格生成

本算例中，采用NACA0012翼型模擬邊界運動。所用靜態(tài)網(wǎng)格是一個準(zhǔn)二維單塊結(jié)構(gòu)網(wǎng)格，如圖2所示,圖中橫縱坐標(biāo)均為無量綱長度。該網(wǎng)格由大小為106×69的二維O型網(wǎng)格沿展向拉伸一個弦長得到，二維O型網(wǎng)格采用雙曲型網(wǎng)格生成法生成，第一層網(wǎng)格與翼表面的距離為弦長的10-3。

本算例中，移動邊界的位置隨時間的變化關(guān)系為

(13)

(14)

式中：α(t)為繞1/4弦線的轉(zhuǎn)動；α0為轉(zhuǎn)動運動幅值；Tm為轉(zhuǎn)動運動周期；xc(t)、yc(t)為翼型形心平動位移；hm為水平方向的最大位移。顯然，這種平動與轉(zhuǎn)動位移的組合運動可用于模擬飛機(jī)外掛物分離問題。

本算例的目的在于驗證本文方法的正確性和收斂性，并對生成網(wǎng)格的質(zhì)量進(jìn)行評估。因此，只生成網(wǎng)格而不求解流場。

設(shè)α0=10.0°,Tm=0.15 s,hm=2.0c，c為翼型的無量綱弦長。在一個周期Tm內(nèi)取300個點，因此時間步長Δt=0.000 5 s。

圖3給出翼型運動過程中4個不同時刻翼面附近網(wǎng)格的局部視圖，可以看到，翼面附近網(wǎng)格的質(zhì)量得到了很好的保留。從定量角度，給出了最大長寬比和平均長寬比的迭代收斂歷程，如圖4和圖5所示；并給出了網(wǎng)格單元角對正交角偏離的統(tǒng)計結(jié)果，如表1所示。由圖4和圖5可以看出，網(wǎng)格最大長寬比和平均長寬比的迭代過程是穩(wěn)定的，在迭代600步時已收斂到與靜態(tài)網(wǎng)格最大長寬比和平均長寬比相當(dāng)?shù)某潭?。從網(wǎng)格角度偏離的統(tǒng)計結(jié)果表1來看，內(nèi)點以及外邊界上的角度偏離比靜態(tài)網(wǎng)格的角度偏離要大，但仍然處于可以接受的程度。情況最好的是翼型表面，幾乎達(dá)到了和靜態(tài)網(wǎng)格一樣的程度。

GridInteriorOutboundarySurfaceofairfoilMaximumAverageMaximumAverageMaximumAverageStaticgrid(t=0s)12.82000.62180.78950.00870.00710.0002Dynamicgrid(t=0.1125s)15.16002.414011.45000.29490.00730.0002Dynamicgrid(t=0.15s)12.79000.89696.35800.11730.00710.0002

3.2 繞俯仰振蕩翼型NACA0012的非定常無黏跨聲速流動

本算例所用靜態(tài)網(wǎng)格與3.1節(jié)所用靜態(tài)網(wǎng)格完全一致。翼型繞1/4弦線的俯仰振蕩由式(15)決定：

α(t)=αm+α0sin(ωt)

(15)

式中：α(t)為瞬時迎角；αm、α0和ω分別為平均迎角、迎角振蕩幅度和角運動圓頻率。角頻率ω與減縮頻率k的關(guān)系定義為

k=ωc/(2U∞)

(16)

式中：U∞為自由來流速度。

同時采用本文提出的動網(wǎng)格生成方法和線段彈簧(Segment Spring)模擬法對以下工況進(jìn)行模擬。

工況AMa∞=0.76,αm=0.016°,α0=2.51°,k=0.081 4。

工況BMa∞=0.60,αm=3.16°,α0=4.59°,k=0.081 4。

圖6～圖9給出了瞬時升力系數(shù)CL和俯仰力矩系數(shù)Cm與試驗結(jié)果[28-29]的比較。除了工況B中的Cm曲線外，2種動網(wǎng)格方法的預(yù)測結(jié)果與試驗結(jié)果都具有很好的一致性。根據(jù)文獻(xiàn)[28]的分析，Cm曲線的差異是由于實驗?zāi)Ｐ偷膶嶋H俯仰力矩中心與理論名義中心之間存在微小差異而引起的。對于實際力矩中心，文獻(xiàn)[28]認(rèn)為應(yīng)該是27.3%而不是25%。圖9也給出了采用實際力矩中心對數(shù)值結(jié)果進(jìn)行修正后的情況，修正后的結(jié)果與實驗結(jié)果具有更好的一致性。

為比較本文方法和線段彈簧模擬法在網(wǎng)格生成上的效率，采用工況B的運動參數(shù)讓翼型做同樣的俯仰振蕩運動，但屏蔽流場求解代碼，在每一時間步長上只生成網(wǎng)格，不求解流場。并指定2種方法的外迭代次數(shù)都為500次,內(nèi)迭代次數(shù)都為200次。網(wǎng)格求解都采用Gauss-Seidel點松弛法。對2種算法花費的總CPU時間進(jìn)行了統(tǒng)計,同時給出了每個網(wǎng)格點每次外迭代所用時間，結(jié)果如表2 所示。實測數(shù)據(jù)說明在同樣的迭代次數(shù)下本文方法在計算效率上落后于彈簧模擬法。但需要注意,這是在指定相同迭代步數(shù)的前提下作出的比較，如果換為以某種網(wǎng)格質(zhì)量指標(biāo)作為約束條件的話，結(jié)果可能會有所不同。也就是說，如果是在達(dá)到同樣網(wǎng)格質(zhì)量的情況下來比較，橢圓型方法就可能因為需要的迭代步數(shù)更少，而表現(xiàn)出更高的計算效率。

表2 2種動網(wǎng)格生成方法執(zhí)行時間比較(迭代步數(shù)=500,網(wǎng)格點數(shù)=34 845)

Table 2 Comparison of execution time for two dynamic grid generation methods (iterations=500, number of grid points=34 845)

MethodTotalcomputationtime/sGridpointiterationtime/sThepresent454.92.61×10-5Springanalogy317.51.82×10-5

為了研究方法的魯棒性，選擇工況B作為研究實例。讓角運動幅度α0按照每次4.59° 的間隔遞增，然后采用2種方法分別求解動網(wǎng)格，直到網(wǎng)格線出現(xiàn)交叉或負(fù)體積網(wǎng)格。當(dāng)α0增加到7×4.59°=32.13° 時，彈簧模擬法首先出現(xiàn)了網(wǎng)格線交叉的情況，如圖10所示。該圖也同時展示了翼型前緣網(wǎng)格的變形情況。在同樣條件下采用本文方法得到的網(wǎng)格如圖11所示，可見，無論在前緣還是后緣附近，網(wǎng)格質(zhì)量明顯高于彈簧模擬法得到的結(jié)果。2種方法的迭代次數(shù)都是500次。這說明本文方法確實具有較高的魯棒性，能夠適應(yīng)較大的邊界運動。

3.3 再入返回艙的俯仰振蕩

通過對再入返回艙(OREX)俯仰振蕩的模擬，演示本文動網(wǎng)格生成方法處理三維多塊結(jié)構(gòu)網(wǎng)格的能力。

OREX的幾何尺寸來自文獻(xiàn)[30]。艙體附近的三維靜態(tài)網(wǎng)格如圖12所示，整個網(wǎng)格被劃分為4塊，各網(wǎng)格塊的規(guī)模分別為51×21×21、51×21×21、67×41×51、67×41×51，其中第1指標(biāo)指示法向，第2指標(biāo)指示流向，第3指標(biāo)指示周向。計算區(qū)域為直徑的30倍，第1層網(wǎng)格距離壁面為直徑的10-3。

算例采用無黏流模型，對返回艙施加式(15)形式的強(qiáng)迫運動，旋轉(zhuǎn)軸為通過質(zhì)心的z軸。根據(jù)Yoshinaga等[31]的實驗結(jié)果，選擇以下兩組參數(shù)進(jìn)行非定常流模擬：

1)Ma∞=0.9,αm=7.0°,α0=1.0°,7.0°，k=0.047 7。

2)Ma∞=1.0,αm=7.0°,α0=1.0°,7.0°，k=0.047 7。

圖13是返回艙旋轉(zhuǎn)7.0° 時在xOy平面內(nèi)的網(wǎng)格。表3給出的是靜導(dǎo)數(shù)Cmθ實驗結(jié)果與計算結(jié)果的比較。考慮到實驗結(jié)果沒有明確給出平衡迎角αm，而且采用的是自由振動法，因此也沒有對應(yīng)的強(qiáng)迫振動迎角振幅α0以及減縮頻率k，這幾個參數(shù)是本文指定的, 可能與實驗參數(shù)之間并沒有準(zhǔn)確的對應(yīng)關(guān)系，表中的結(jié)果可用于定性比較。

此處，利用三維算例再次對網(wǎng)格更新時間進(jìn)行了測算。同樣，屏蔽流場求解，在每一時間步長上只生成網(wǎng)格，不求解流場。并指定本文算法和彈簧模擬算法的外迭代次數(shù)都為200次,內(nèi)迭代次數(shù)仍然為200次。結(jié)果如表4所示，其中，每個網(wǎng)格點單次迭代的時間開銷與表2的結(jié)果十分接近，驗證了表2的結(jié)果。

圖14是返回艙旋轉(zhuǎn)45.0° 時在xOy平面內(nèi)的網(wǎng)格。說明這種方法在三維大位移情況也能得到較好的結(jié)果。

表3 靜導(dǎo)數(shù)Cmθ實驗結(jié)果與計算結(jié)果的比較

Table 3 Comparison of computational and experimental results of static derivativeCmθ

NominalMachnumberα0/(°)StaticderivativeCmθ/rad-1ComputationExperiment[31]0.91.0-0.1437.0-0.127-0.1541.01.0-0.1137.0-0.098-0.148

表4 2種動網(wǎng)格生成方法執(zhí)行時間比較(迭代步數(shù)=200,網(wǎng)格點數(shù)=325 176)

Table 4 Comparison of excution time for two dynamic grid generation methods (iterations=200, number of grid points=325 716)

MethodTotalcomputationtime/sGridpointiterationstime/sThepresent1783.62.73×10-5Springanalogy1345.52.06×10-6

3.4 黏性網(wǎng)格算例

原理上，使用式(8)提取靜態(tài)網(wǎng)格源項的方法并不僅限于無黏網(wǎng)格，黏性網(wǎng)格同樣適用。不同之處在于,黏性動網(wǎng)格的生成需要特別關(guān)注初始迭代向量,即每一時間步長下網(wǎng)格迭代初場的選取以及物面銳角邊界的處理問題[32]。這是容易導(dǎo)致黏性動網(wǎng)格收斂失敗的2個主要因素。第1個因素主要是由于黏性網(wǎng)格中高梯度區(qū)域的存在使求解網(wǎng)格的偏微分方程非線性程度加劇，收斂域縮小[27]而引起的。解決的辦法一般是通過減小時間步長來控制物面位移量，采用第n步已收斂的網(wǎng)格作為第n+1步的迭代初場，盡量使迭代初場位于第n+1步網(wǎng)格的收斂域內(nèi)。因此, 黏性動網(wǎng)格的生成往往需要考慮時間步長的約束,尤其是在物面存在凸銳角邊界時這個問題會更突出。

此處，以一個帶凸銳角邊界的RAE2822翼型為例,驗證本文方法對黏性網(wǎng)格的適應(yīng)性。翼型繞1/4弦線的俯仰振蕩關(guān)系仍然由式(15)決定。取振蕩幅度α0=5°,振蕩頻率f=10 Hz(ω=2πf)。C形靜態(tài)網(wǎng)格如圖15和圖16所示，網(wǎng)格數(shù)據(jù)來自文獻(xiàn)[33]。

圖17～圖19顯示了翼型繞1/4弦線旋轉(zhuǎn)5°后前緣及后緣部分的局部網(wǎng)格與原始靜態(tài)網(wǎng)格(圖20和圖21)相比，除后緣部分的網(wǎng)格質(zhì)量有所下降外，其他物面網(wǎng)格質(zhì)量仍得到了較好的保持。后緣部分網(wǎng)格質(zhì)量下降的原因，主要是由于凸銳角邊界(斜率間斷)引起。而對間斷較為敏感其實是微分方程法本身固有的弱點,這也是導(dǎo)致在間斷附近數(shù)值解收斂十分困難的原因。

對這個問題，文獻(xiàn)[32]建議采用人工固定凸銳角邊界附近網(wǎng)格點的辦法來解決。但由于動網(wǎng)格生成過程需要與流場求解器相融合，難以提供有效的人機(jī)交互接口來指定需要“固定”的網(wǎng)格點，因而難以實施。本文建議在原始靜態(tài)網(wǎng)格生成時，對帶有凸銳角的物面邊界，首先應(yīng)考慮生成具有O型拓?fù)涞木W(wǎng)格，并對凸銳角采取“打鈍”的辦法，預(yù)先用1～2 mm左右的圓弧對銳角進(jìn)行圓角化，避免凸銳角的出現(xiàn)。對只能采用C型拓?fù)涞木W(wǎng)格，在凸銳角邊界附近無法采用鈍化處理?？煽紤]采用其他策略，例如，考慮在迭代初期采用最簡單的線段彈簧法生成動網(wǎng)格“初場”，然后由橢圓型方程方法繼續(xù)進(jìn)行迭代，獲得更加光順的網(wǎng)格。彈簧法與橢圓型方程法的迭代步數(shù)可按3∶7 進(jìn)行分配。圖22是采用以上策略后得到的后緣部分的網(wǎng)格，與圖19相比，網(wǎng)格質(zhì)量得到了明顯改善，甚至部分恢復(fù)到了與初始靜態(tài)網(wǎng)格(圖21)相當(dāng)?shù)某潭取?/p>

4 結(jié) 論

1) 基于求解橢圓型偏微分方程方法，實現(xiàn)了一種簡單高效的結(jié)構(gòu)動網(wǎng)格重構(gòu)方法，數(shù)值模擬結(jié)果證明了方法的正確性。

2) 在同樣迭代步數(shù)約束下，本文方法花費的網(wǎng)格求解時間略長于傳統(tǒng)彈簧模擬法，但網(wǎng)格質(zhì)量較好，方法的魯棒性優(yōu)于彈簧模擬法。

3) 文中采用繼承靜態(tài)網(wǎng)格源項，并在動網(wǎng)格生成過程中保持不變的策略，既解決了動網(wǎng)格生成過程中如何確定網(wǎng)格源項的問題，又節(jié)約了搜索源項需要的外迭代過程，是方法成功的關(guān)鍵。

[1] STEGER J L, DOUGHERTY F C, BENEK J A. A chimera grid scheme[C]//Mini Symposium on Advances in Grid Generation， 1982.

[2] CHAN W M, PANDYAY S A, ROGERS S E. Effcient creation of overset grid hole boundaries and effects of their locations on aerodynamic loads[C]//21st AIAA Computational Fluid Dynamics Conference. Reston: AIAA, 2013.

[3] KIM N, CHAN W M. Automation of hole-cutting for overset grids using the x-rays approach: AIAA-2011-3052[R]. Reston: AIAA, 2011.

[4] 王文, 閻超. 新型重疊網(wǎng)格洞面優(yōu)化方法及其應(yīng)用[J]. 航空學(xué)報, 2016, 37(3): 826-835. WANG W, YAN C. Novel overlapping optimization algorithm of overlapping grid and its applications[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(3):826-835 (in Chinese).

[5] 張偉偉, 高傳強(qiáng), 葉正寅. 氣動彈性計算中網(wǎng)格變形方法研究進(jìn)展[J]. 航空學(xué)報, 2014, 35(2): 303-319. ZHANG W W, GAO C Q, YE Z Y. Reseach progress on mesh deformation in computational aeroelasticity[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(3): 303-319 (in Chinese).

[6] FOSTER N F. Accuracy of high-order cfd and overset interpolation in finite volumee/difference codes[C]//22nd AIAA Computational Fluid Dynamics Conference. Reston: AIAA, 2015.

[7] UCHIDA Y, KANDA H. Spacing control of 3-D transfinite interpolation grid generation: AIAA-2005-5243[R]. Reston: AIAA, 2005.

[8] BATINA J T. Unsteady Euler airfoil solutions using unstructured dynamic meshes[J]. AIAA Journal, 1990, 28(8): 1381-1388.

[9] BLOM F J. Considerations on the spring analogy[J]. International Journal of Numerical Methods in Fluid, 2000, 32(6): 647-668.

[10] CANTARITI D L, GRIBBEN W M, BADCOCK K J, et al. A grid deformation technique for unsteady flow computations[J]. International Journal of Numerical Methods in Fluids, 2000, 32(3): 285-311.

[11] 劉楓. 動網(wǎng)格技術(shù)研究及其在高超聲速流動中的應(yīng)用[D]. 長沙：國防科學(xué)技術(shù)大學(xué)， 2009: 54-55. LIU F. Investigations of dynamic grid generation and its applications for supersonic flow[D]. Changsha: National University of Defense Technology, 2009:54-55 (in Chinese).

[12] 郭正, 劉君, 翟章華. 用非結(jié)構(gòu)動網(wǎng)格方法模擬有相對運動的多體繞流[J]. 空氣動力學(xué)學(xué)報, 2001, 19(3): 310-316. GUO Z, LIU J, ZHAI Z H. Simulation of flows past multi-body in relative motion with dynamic unstructured method[J]. Acta Aerodynamic Sinica, 2001, 19(3): 310-316 (in Chinese).

[13] TEZDUYAR T E. Stability finite element formulations for incompressible flow computations[J]. Advances in Applied Mechanics, 1992, 28(1): 1-44.

[14] LIU X, QIN N, XIA H. Fast dynamic grid deformation based on delaunay graph mapping[J]. Journal of Computational Physics, 2006, 211(2): 405-423.

[15] BORE A, SCHOOT M S, FACULTY S H. Mesh deformation based on radial basis function interpolation[J]. Computers and Structures, 2007, 85(11): 784-795.

[16] RENDALL T C S, ALLEN C B. Efficient mesh motion using radial basis functions with data reduction algorithms[J]. Journal of Computational Physics, 2010, 229(8): 2810-2820.

[17] WANG G, HARIS H M, YE Z Y. An improved point selection method for hybrid unstructured mesh deformation using radial basis functions: AIAA-2013-3076[R]. Reston: AIAA, 2013.

[18] 魏其, 李春娜, 谷良賢, 等. 一種基于徑向函數(shù)和峰值選擇法的高效網(wǎng)格變形技術(shù)[J]. 航空學(xué)報, 2016, 37(7): 1401-1414. WEI Q, LI C N, GU L X, et al. An efficient mesh deformation method based on radial basis functions and peak-selection method[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(7): 1401-1414 (in Chinese).

[19] 孫巖, 鄧小剛, 王光學(xué), 等. 基于徑向基函數(shù)改進(jìn)的Delaunay圖映射動網(wǎng)格方法[J]. 航空學(xué)報, 2014, 35(3): 727-735. SUN Y, DENG X G, WANG G X, et al. An improvement on Delaunay graph mapping dynamic grid method based on radial basis functions[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(3): 727-735 (in Chinese).

[20] JOHNSON A A, TEZDUYAR T E. Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces[J]. Computer Methods in Applied Mechanics and Engineering, 1994, 119(1-2): 73-94.

[21] USHIJIMA S. Three-dimensional arbitrary Lagrangian-Eulerian numerical prediction method for non-linear free surface oscillation[J]. International Journal for Numerical Methods in Fluids, 1998, 26(5): 605-623.

[22] GRüBER B, CARTENS V. Computation of the unsteady transonic flow in harmonically oscillating turbine cascades taking into account viscous effects[J]. Journal of Turbomachinery, 1998, 120(1): 104-111.

[23] THOMAS P D, MIDDLECOFF J F. Direct control of the grid point distribution in meshes generated by elliptic equations[J]. AIAA Journal, 1979, 18(6): 652-656.

[24] SORENSON R L, STEGER J L. Numerical generation of 2D grids by the use of poisson equations with grids control: N81-14722[R]. 1981.

[25] SORENSON R L. A computer program to generate 2D grids about airfoil and other shapes by the use of poisson’s equation: N80-26266[R]. 1980.

[26] HILGENSTOCK A. A fast method for the elliptic generation of three-dimensional grids with full boundary control[C]//Numerical Grid Generation in Computational Fluid Mechanics, 1988: 137-146.

[27] SONAR T. Grid generation using elliptic partial differential equations: DFVLR-FB89-15[R]. 1989.

[28] GAITONDE A L, FIDDES S P. A moving mesh system for the calculation of unsteady flows: AIAA-1993-0641[R]. Reston: AIAA, 1993.

[29] LI J, HUANG S. Unsteady viscous flow simulations by a fully implicit method with deforming grid: AIAA-2005-1221[R]. Reston: AIAA, 2005.

[30] UNMEEL B M. Synthesis of contributed simulations for OREX test cases: NASA/TM-1998-112238[R]. Washington, D.C.: NASA, 1998.

[31] YOSHINAGA T, TATE A, WATANABE M. Dynamic test of the orbital reentry vehicle (OREX) in a transonic wind tunnel with comparison to flight data: AIAA-1995-1901-CP[R]. Reston: AIAA, 1995.

[32] THOMPSON J F. Numerical grid generation-foundations and applications[M]. North Holland: Amsterdam, 1985.

[33] SLATER J W. RAE2822 transonic airfoil: Study #1[EB/OL]. (2015-01-22)[2016-07-17]. http://www.grc.nasa.gov/WWW/wind/valid/raetaf/raetaf.html.

(責(zé)任編輯：李明敏)

URL：www.cnki.net/kcms/detail/11.1929.V.20161124.1034.002.html

＊Corresponding author. E-mail: lufengling126@126.com

A dynamic structured grid generation method based onsolving elliptic equations

LU Fengling1,2,*, CHEN Xiaoqian1, YU Caihui2, MIAO Meng2

1.CollegeofAerospaceScienceandEngineering,NationalUniversityofDefenceTechnology,Changsha410073,China2.ScienceandTechnologyonSpacePhysicsLaboratory,Beijing100076,China

A simple robust structured dynamic grid generation method based on the solution of elliptic partial differential equations is developed for computing the unsteady flows with moving boundaries. In the method, the source terms of the Poisson equation which can control the spacing and the orthogonality of the grid are inherited from the known static grid, and held fixed throughout the process of dynamic grid generation. With the process, the outer iterations for determining the source terms usually needed in the elliptic grid generation can thus be saved, and no adjustable parameters are required to be prescribed. This makes the method more efficient and easy to implement in an existing CFD code. The numerical results demonstrate that the proposed method can provide an efficient way of deforming the grid based on solving elliptic partial differential equations. The orthogonality and smoothness of original static grid can be maintained well by the proposed method. When the same number of iterations are given as the constraint condition, the grid generation efficiency of the method is lower than that of the spring analogy, but the robustness of the method is superior to the spring analogy.

dynamic grid generation; elliptic grid generation; computational fluid dynamics; numerical grid generation; grid deformation approach

2016-07-21； Revised：2016-08-11； Accepted：2016-11-20； Published online：2016-11-24 10:34

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0303

2016-07-21； 退修日期：2016-08-11； 錄用日期：2016-11-20； 網(wǎng)絡(luò)出版時間：2016-11-24 10:34

www.cnki.net/kcms/detail/11.1929.V.20161124.1034.002.html

＊通訊作者.E-mail: lufengling126@126.com

盧鳳翎， 陳小前， 禹彩輝， 等． 一種基于求解橢圓型方程的結(jié)構(gòu)動網(wǎng)格生成方法[J]． 航空學(xué)報, 2017, 38(3): 120632. LU F L, CHEN X Q, YU C H, et al. A dynamic structured grid generation method based on solving elliptic equations[J]. Acta Aeronau-tica et Astronautica Sinica, 2017, 38(3): 120632.

V211.3

A

1000-6893(2017)03-120632-14