一種基于求解橢圓型方程的結構動網格生成方法

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

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

一種基于求解橢圓型方程的結構動網格生成方法

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

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

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

動網格生成； 橢圓型網格生成； 計算流體力學； 數值網格生成； 網格變形方法

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

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

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

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

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

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

1 一般橢圓型網格生成法

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

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

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

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

(1)

映射函數被要求滿足泊松方程：

(2)

式中：P、Q、R為強迫函數，也稱為源項。

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

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

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

(3)

式中：

(4)

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

(5)

式(3)經整理后可寫為

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

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

(6)

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

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

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

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

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

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

2 本文動網格生成算法

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

2.1 源項的獲取

顯然，無論采用何種方法生成的網格，它首先應該是規(guī)則的。否則，它將不適宜被用來模擬一個實際的流動問題。剩下的問題是，是否可以在泊松類方程的解空間中找到與此相同的一個網格。對這個問題先不做正面回答。觀察式(6)的網格生成系統(tǒ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)又可寫成標量形式：

(8)

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

(9)

其他2個坐標方向η和ζ的導數也有類似的表達。

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

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

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

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

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

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

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

(10)

(11)

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

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

(12)

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

除了物面邊界和遠場邊界外，網格中一般還有內部邊界，即塊與塊之間的對接面，或者是O形以及C形拓撲中具有相同的物理空間坐標，但分屬不同的計算坐標面而具有邊界含義的網格面(有的資料稱為坐標切口, Coordinate Cut)。這種邊界存在于網格內部，在新的網格生成過程中，不能證明式(12)的存在。但從源項控制的角度看，源項直接控制的是邊界附近網格的角度和網格間距，如果只要求內部邊界附近的網格角度和網格間距得到保持，而不關心內部邊界的實際位置的話，完全可以在內部邊界上人為指定式(12)成立。

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

2.2 動網格生成求解過程

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

1) 假定邊界源項初值。

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

3) 將邊界上新的源項值插值到內部網格點處。

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

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

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

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

1) 讀入已知靜網格。

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

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

4) 采用上一時間步的網格坐標值作為當前時間步的初始解。

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

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

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

3 算法驗證

3.1 帶平移和轉動的NACA0012翼型動網格生成

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

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

(13)

(14)

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

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

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

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

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)網格與3.1節(jié)所用靜態(tài)網格完全一致。翼型繞1/4弦線的俯仰振蕩由式(15)決定：

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

(15)

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

k=ωc/(2U∞)

(16)

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

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

工況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給出了瞬時升力系數CL和俯仰力矩系數Cm與試驗結果[28-29]的比較。除了工況B中的Cm曲線外，2種動網格方法的預測結果與試驗結果都具有很好的一致性。根據文獻[28]的分析，Cm曲線的差異是由于實驗模型的實際俯仰力矩中心與理論名義中心之間存在微小差異而引起的。對于實際力矩中心，文獻[28]認為應該是27.3%而不是25%。圖9也給出了采用實際力矩中心對數值結果進行修正后的情況，修正后的結果與實驗結果具有更好的一致性。

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

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

3.3 再入返回艙的俯仰振蕩

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

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

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

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是返回艙旋轉7.0° 時在xOy平面內的網格。表3給出的是靜導數Cmθ實驗結果與計算結果的比較?？紤]到實驗結果沒有明確給出平衡迎角αm，而且采用的是自由振動法，因此也沒有對應的強迫振動迎角振幅α0以及減縮頻率k，這幾個參數是本文指定的, 可能與實驗參數之間并沒有準確的對應關系，表中的結果可用于定性比較。

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

圖14是返回艙旋轉45.0° 時在xOy平面內的網格。說明這種方法在三維大位移情況也能得到較好的結果。

表3 靜導數Cmθ實驗結果與計算結果的比較

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種動網格生成方法執(zhí)行時間比較(迭代步數=200,網格點數=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 黏性網格算例

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

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

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

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

4 結 論

1) 基于求解橢圓型偏微分方程方法，實現了一種簡單高效的結構動網格重構方法，數值模擬結果證明了方法的正確性。

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

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

[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] 王文, 閻超. 新型重疊網格洞面優(yōu)化方法及其應用[J]. 航空學報, 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] 張偉偉, 高傳強, 葉正寅. 氣動彈性計算中網格變形方法研究進展[J]. 航空學報, 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]. 長沙：國防科學技術大學， 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] 郭正, 劉君, 翟章華. 用非結構動網格方法模擬有相對運動的多體繞流[J]. 空氣動力學學報, 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] 魏其, 李春娜, 谷良賢, 等. 一種基于徑向函數和峰值選擇法的高效網格變形技術[J]. 航空學報, 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] 孫巖, 鄧小剛, 王光學, 等. 基于徑向基函數改進的Delaunay圖映射動網格方法[J]. 航空學報, 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.

(責任編輯：李明敏)

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； 網絡出版時間：2016-11-24 10:34

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

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

盧鳳翎， 陳小前， 禹彩輝， 等． 一種基于求解橢圓型方程的結構動網格生成方法[J]． 航空學報, 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