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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

(2)

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

通過交換自變量(x，y，z)與因變量(ξ，η，ζ)的位置，在變換后的計(jì)算域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)格點(diǎn)。內(nèi)點(diǎn)單元的大小和形狀通過調(diào)節(jié)泊松方程中的源項(xiàng)來實(shí)現(xiàn)。顯然，源項(xiàng)在生成適宜精確模擬流動(dòng)問題的網(wǎng)格時(shí)起到了關(guān)鍵性的作用。

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

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

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

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

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

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

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

2.1 源項(xiàng)的獲取

顯然，無論采用何種方法生成的網(wǎng)格，它首先應(yīng)該是規(guī)則的。否則，它將不適宜被用來模擬一個(gè)實(shí)際的流動(dòng)問題。剩下的問題是，是否可以在泊松類方程的解空間中找到與此相同的一個(gè)網(wǎng)格。對這個(gè)問題先不做正面回答。觀察式(6)的網(wǎng)格生成系統(tǒng)，可以看到其中所有導(dǎo)數(shù)項(xiàng)以及度量系數(shù)都可以由已知網(wǎng)格通過差分方法計(jì)算得到，系統(tǒng)中未知的是源項(xià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ù)項(xiàng)可以表示成中心差分的形式：

(9)

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

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

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

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

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

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

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

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

(10)

(11)

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

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

(12)

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

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

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

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

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

1) 假定邊界源項(xiàng)初值。

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

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

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

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

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

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

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

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

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

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

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

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

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

3 算法驗(yàn)證

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

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

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

(13)

(14)

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

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

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

圖3給出翼型運(yùn)動(dòng)過程中4個(gè)不同時(shí)刻翼面附近網(wǎng)格的局部視圖，可以看到，翼面附近網(wǎng)格的質(zhì)量得到了很好的保留。從定量角度，給出了最大長寬比和平均長寬比的迭代收斂歷程，如圖4和圖5所示；并給出了網(wǎng)格單元角對正交角偏離的統(tǒng)計(jì)結(jié)果，如表1所示。由圖4和圖5可以看出，網(wǎng)格最大長寬比和平均長寬比的迭代過程是穩(wěn)定的，在迭代600步時(shí)已收斂到與靜態(tài)網(wǎng)格最大長寬比和平均長寬比相當(dāng)?shù)某潭?。從網(wǎng)格角度偏離的統(tǒng)計(jì)結(jié)果表1來看，內(nèi)點(diǎn)以及外邊界上的角度偏離比靜態(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的非定常無黏跨聲速流動(dòng)

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

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

(15)

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

k=ωc/(2U∞)

(16)

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

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

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

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

3.3 再入返回艙的俯仰振蕩

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

OREX的幾何尺寸來自文獻(xiàn)[30]。艙體附近的三維靜態(tài)網(wǎng)格如圖12所示，整個(gè)網(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)指示周向。計(jì)算區(qū)域?yàn)橹睆降?0倍，第1層網(wǎng)格距離壁面為直徑的10-3。

算例采用無黏流模型，對返回艙施加式(15)形式的強(qiáng)迫運(yùn)動(dòng)，旋轉(zhuǎn)軸為通過質(zhì)心的z軸。根據(jù)Yoshinaga等[31]的實(shí)驗(yàn)結(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° 時(shí)在xOy平面內(nèi)的網(wǎng)格。表3給出的是靜導(dǎo)數(shù)Cmθ實(shí)驗(yàn)結(jié)果與計(jì)算結(jié)果的比較?？紤]到實(shí)驗(yàn)結(jié)果沒有明確給出平衡迎角αm，而且采用的是自由振動(dòng)法，因此也沒有對應(yīng)的強(qiáng)迫振動(dòng)迎角振幅α0以及減縮頻率k，這幾個(gè)參數(shù)是本文指定的, 可能與實(shí)驗(yàn)參數(shù)之間并沒有準(zhǔn)確的對應(yīng)關(guān)系，表中的結(jié)果可用于定性比較。

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

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

表3 靜導(dǎo)數(shù)Cmθ實(shí)驗(yàn)結(jié)果與計(jì)算結(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種動(dòng)網(wǎng)格生成方法執(zhí)行時(shí)間比較(迭代步數(shù)=200,網(wǎng)格點(diǎn)數(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)格源項(xiàng)的方法并不僅限于無黏網(wǎng)格，黏性網(wǎng)格同樣適用。不同之處在于,黏性動(dòng)網(wǎng)格的生成需要特別關(guān)注初始迭代向量,即每一時(shí)間步長下網(wǎng)格迭代初場的選取以及物面銳角邊界的處理問題[32]。這是容易導(dǎo)致黏性動(dòng)網(wǎng)格收斂失敗的2個(gè)主要因素。第1個(gè)因素主要是由于黏性網(wǎng)格中高梯度區(qū)域的存在使求解網(wǎng)格的偏微分方程非線性程度加劇，收斂域縮小[27]而引起的。解決的辦法一般是通過減小時(shí)間步長來控制物面位移量，采用第n步已收斂的網(wǎng)格作為第n+1步的迭代初場，盡量使迭代初場位于第n+1步網(wǎng)格的收斂域內(nèi)。因此, 黏性動(dòng)網(wǎng)格的生成往往需要考慮時(shí)間步長的約束,尤其是在物面存在凸銳角邊界時(shí)這個(gè)問題會(huì)更突出。

此處，以一個(gè)帶凸銳角邊界的RAE2822翼型為例,驗(yàn)證本文方法對黏性網(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ì)量下降的原因，主要是由于凸銳角邊界(斜率間斷)引起。而對間斷較為敏感其實(shí)是微分方程法本身固有的弱點(diǎn),這也是導(dǎo)致在間斷附近數(shù)值解收斂十分困難的原因。

對這個(gè)問題，文獻(xiàn)[32]建議采用人工固定凸銳角邊界附近網(wǎng)格點(diǎn)的辦法來解決。但由于動(dòng)網(wǎng)格生成過程需要與流場求解器相融合，難以提供有效的人機(jī)交互接口來指定需要“固定”的網(wǎng)格點(diǎn)，因而難以實(shí)施。本文建議在原始靜態(tài)網(wǎng)格生成時(shí)，對帶有凸銳角的物面邊界，首先應(yīng)考慮生成具有O型拓?fù)涞木W(wǎng)格，并對凸銳角采取“打鈍”的辦法，預(yù)先用1～2 mm左右的圓弧對銳角進(jìn)行圓角化，避免凸銳角的出現(xiàn)。對只能采用C型拓?fù)涞木W(wǎng)格，在凸銳角邊界附近無法采用鈍化處理?？煽紤]采用其他策略，例如，考慮在迭代初期采用最簡單的線段彈簧法生成動(dò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) 基于求解橢圓型偏微分方程方法，實(shí)現(xiàn)了一種簡單高效的結(jié)構(gòu)動(dòng)網(wǎng)格重構(gòu)方法，數(shù)值模擬結(jié)果證明了方法的正確性。

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

3) 文中采用繼承靜態(tài)網(wǎng)格源項(xiàng)，并在動(dòng)網(wǎng)格生成過程中保持不變的策略，既解決了動(dòng)網(wǎng)格生成過程中如何確定網(wǎng)格源項(xiàng)的問題，又節(jié)約了搜索源項(xiàng)需要的外迭代過程，是方法成功的關(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é)報(bào), 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), 葉正寅. 氣動(dòng)彈性計(jì)算中網(wǎng)格變形方法研究進(jìn)展[J]. 航空學(xué)報(bào), 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] 劉楓. 動(dòng)網(wǎng)格技術(shù)研究及其在高超聲速流動(dòng)中的應(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)動(dòng)網(wǎng)格方法模擬有相對運(yùn)動(dòng)的多體繞流[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 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é)報(bào), 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圖映射動(dòng)網(wǎng)格方法[J]. 航空學(xué)報(bào), 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ò)出版時(shí)間：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)動(dòng)網(wǎng)格生成方法[J]． 航空學(xué)報(bào), 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