改進(jìn)的動彈網(wǎng)格方法在航空氣彈計(jì)算中的應(yīng)用

胡 凡 ,范銳軍

(1.陜西省科技資源統(tǒng)籌中心，陜西 西安 710077)

(2.陜西西安現(xiàn)代控制技術(shù)研究所，陜西 西安 710065)

改進(jìn)的動彈網(wǎng)格方法在航空氣彈計(jì)算中的應(yīng)用

胡 凡1,范銳軍2

(1.陜西省科技資源統(tǒng)籌中心，陜西 西安 710077)

(2.陜西西安現(xiàn)代控制技術(shù)研究所，陜西 西安 710065)

氣動彈性是現(xiàn)代航空氣動力計(jì)算中一個突出的問題。主要研究基于Delaunay圖映射方法的動彈網(wǎng)格的歐拉方程CFD計(jì)算及其在航空標(biāo)模M6機(jī)翼上的靜氣動彈性應(yīng)用。以Delaunay圖映射方法為基礎(chǔ)，針對三維非結(jié)構(gòu)運(yùn)動網(wǎng)格技術(shù)進(jìn)行了研究、開發(fā)和改進(jìn)，同時利用計(jì)算流體力學(xué)的方法, 開發(fā)了一套適用性較好的非結(jié)構(gòu)網(wǎng)格歐拉方程流場求解器，進(jìn)一步通過流固耦合的力學(xué)方法，對航空標(biāo)模M6機(jī)翼的靜氣動彈性問題進(jìn)行了研究和分析，給出了CFD并行計(jì)算的設(shè)計(jì)方法及算例。

Delaunay圖映射動彈網(wǎng)格方法；歐拉方程；CFD/CSD耦合技術(shù)；靜氣動彈性

目前針對基于Euler方程或Navier-Stokes(N-S)方程等流動控制方程的CFD計(jì)算流固耦合問題，普遍采用有限體積的空間離散方法。在眾多的網(wǎng)格方法中，非結(jié)構(gòu)網(wǎng)格的應(yīng)用非常廣泛，它對于不同拓?fù)浣Y(jié)構(gòu)的復(fù)雜氣動外形具有公認(rèn)的普適性和靈活性。由于在航空氣彈計(jì)算中，氣動模型的外形隨著氣動彈性現(xiàn)象不斷變化，網(wǎng)格邊界上的網(wǎng)格點(diǎn)也需要根據(jù)氣動模型的外形變化加以變形和調(diào)整，甚至適時進(jìn)行重構(gòu)，所以網(wǎng)格的運(yùn)動技術(shù)成為了提高計(jì)算方法適應(yīng)性的關(guān)鍵。

本文以航空標(biāo)模M6機(jī)翼作為算例，通過自行開發(fā)的氣動彈性程序進(jìn)行靜氣動彈性的計(jì)算及分析。通過使用基于Delaunay圖映射的動彈網(wǎng)格技術(shù)，成功地對多塊混合非結(jié)構(gòu)網(wǎng)格系統(tǒng)進(jìn)行了運(yùn)動變形，實(shí)現(xiàn)了三維非結(jié)構(gòu)網(wǎng)格的合理運(yùn)動。在此基礎(chǔ)上，結(jié)合結(jié)構(gòu)有限元分析，耦合結(jié)構(gòu)靜平衡方程，進(jìn)行了數(shù)值模擬流固耦合計(jì)算。

1 動彈網(wǎng)格技術(shù)

隨著計(jì)算模型幾何外形復(fù)雜度的增加，網(wǎng)格的變形和重構(gòu)、氣彈耦合計(jì)算的效率以及計(jì)算程序的代數(shù)魯棒性，是動彈網(wǎng)格技術(shù)中同時作用又相互制約的幾個主要因素。針對大形變和多塊混合網(wǎng)格(例如復(fù)雜外形的多塊劃分網(wǎng)格)的運(yùn)動變形，現(xiàn)有的幾種網(wǎng)格運(yùn)動方法，例如TFI代數(shù)網(wǎng)格重構(gòu)生成方法、彈簧近似方法Spring Analogy Method等[1]，都具有計(jì)算效率低、魯棒性差等缺陷，而且在運(yùn)動過程中還易出現(xiàn)網(wǎng)格的交叉和產(chǎn)生負(fù)體積[2]，嚴(yán)重限制了航空氣動彈性求解方法的進(jìn)步。

本文中作者自主開發(fā)并使用Delaunay圖映射方法來進(jìn)行動彈網(wǎng)格的構(gòu)造和變形[3]，在一定程度上，保證了變形前后網(wǎng)格的拓?fù)浣Y(jié)構(gòu)和密度分布的一致性，同時也避免了網(wǎng)格交叉重疊現(xiàn)象的出現(xiàn)，運(yùn)動變形之后的網(wǎng)格具有良好的網(wǎng)格品質(zhì)，具備一定的高效性和良好的魯棒性。

1.1動彈網(wǎng)格系數(shù)的確定

圖1解釋了動彈網(wǎng)格中一個重要的系數(shù)——面積比系數(shù)的確定方法。為了簡化說明，本文采用二維情況，三維情況可以類推得到。

如圖1中所示，網(wǎng)格節(jié)點(diǎn)G包含在三角形ABC中，面積比系數(shù)定義如下：

ei=Si/Si=1,2,3

式中：Si和S為三角形面積。

1.2動彈網(wǎng)格的計(jì)算方法

動彈網(wǎng)格的計(jì)算流程，主要包括以下幾步[4]：

步驟1，將位于網(wǎng)格邊界上的網(wǎng)格節(jié)點(diǎn)，定義為動彈網(wǎng)格特征點(diǎn)，并基于這些特征點(diǎn)來定義計(jì)算區(qū)域。對一組給定的網(wǎng)格特征點(diǎn)，總是存在著一組與之唯一對應(yīng)的Delaunay圖三角形，邊界會在映射回來時自動保持，原始網(wǎng)格的拓?fù)浣Y(jié)構(gòu)會在移動和映射過程中保持不變。

步驟2，步驟1中最終生成的Delaunay圖三角形集合與氣彈計(jì)算的計(jì)算域能夠完全重合，因此三角形區(qū)域能夠包含網(wǎng)格中的任意一個節(jié)點(diǎn)。對于一個給定三角形單元內(nèi)的網(wǎng)格節(jié)點(diǎn)，只有當(dāng)它所有的面積比系數(shù)大于等于零時，該節(jié)點(diǎn)才包含于這個三角形單元。根據(jù)1.1節(jié)中定義的面積比系數(shù)，可建立Delaunay圖三角形與任意網(wǎng)格節(jié)點(diǎn)的映射關(guān)系。此時，Delaunay圖中的特征點(diǎn)能夠真實(shí)準(zhǔn)確地描述模型物面的邊界運(yùn)動。

步驟3，根據(jù)氣動彈性計(jì)算中氣動模型的形變規(guī)律，首先對與網(wǎng)格邊界節(jié)點(diǎn)一一對應(yīng)的Delaunay圖三角形進(jìn)行變形，然后再根據(jù)步驟2中建立的Delaunay圖三角形單元頂點(diǎn)坐標(biāo)與邊界網(wǎng)格節(jié)點(diǎn)坐標(biāo)之間的對應(yīng)關(guān)系，更新變形后網(wǎng)格節(jié)點(diǎn)的坐標(biāo)。注意：在此過程中，動彈網(wǎng)格的面積比系數(shù)不變。

如圖2～圖4所示，利用Delaunay圖映射方法對大位移的混合網(wǎng)格進(jìn)行變形(注：Delaunay圖映射網(wǎng)格變形方法是非數(shù)值迭代方法)，變形后的效果較好，即在網(wǎng)格質(zhì)量得到了較好的保證的同時，保證了計(jì)算效率和計(jì)算的穩(wěn)定性。

2 非結(jié)構(gòu)運(yùn)動網(wǎng)格在靜氣彈求解中的應(yīng)用

2.1基于非結(jié)構(gòu)網(wǎng)格的流場求解

本文CFD解算控制方程是基于ALE描述的守恒型三維可壓縮歐拉方程，其積分形式如下:

式中：Ω為控制體，?Ω表示控制體單元的邊界；dV是體積微元；n是控制體邊界外法向單位向量；dS是面積微元；守恒變量項(xiàng)Q為

其中：ρ表示流體的密度;u,v,w分別為x,y和z軸方向的流體速度;e為單位體積流體的總內(nèi)能。

對流通量矢量項(xiàng)F(Q)在x,y和z軸的分量為：

逆變速度矢量Ψ=Ui+Vj+Wk,沿x,y和z軸的分量U,V和W定義為：

式中：xt,yt和zt分別為網(wǎng)格沿x,y和z軸方向的網(wǎng)格速度。

此外對于理想氣體壓強(qiáng)p為：

式中：γ為理想氣體比熱比，這里取γ=1.4。

本文采用格心格式的有限體積法對CFD解算控制方程進(jìn)行空間上的離散，采用經(jīng)典四步龍格庫塔方法進(jìn)行虛擬時間上的離散，使用壁面函數(shù)對無粘求解進(jìn)行了粘性修正，并且使用焓阻尼方法、當(dāng)?shù)貢r間步長、隱式殘值光順等加速收斂技術(shù)[5]，實(shí)現(xiàn)對三維氣動彈性問題的流場求解。

2.2靜氣動彈性求解流程

本文使用CFD控制方程和結(jié)構(gòu)靜平衡方程耦合求解的方法對靜氣動彈性問題進(jìn)行求解。利用氣動力方程和結(jié)構(gòu)靜平衡方程耦合迭代求解，得到收斂的結(jié)果。

結(jié)構(gòu)影響系數(shù)法是普遍適用于靜氣動彈性問題求解的一種方法[6]。結(jié)構(gòu)影響系數(shù)法的主要原理描述如下。

結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)(x,y)處的變形為

式中：W(x,y)為節(jié)點(diǎn)(x,y)處的法向變形；Czz(x,y,ξ,η)為機(jī)翼(ξ,η)處作用單位的力在節(jié)點(diǎn)(x,y)處發(fā)生的位移；F(ξ,η)為作用在(ξ,η)處的載荷，分為氣動力和慣性力兩部分，即

F(ξ,η)=L(ξ,η)-Nm(ξ,η)g

其中:L(ξ,η)為作用在(ξ,η)處的氣動力，由CFD控制方程求解出的氣動網(wǎng)格流場值進(jìn)行三維面插值得到；N為過載系數(shù)；m(ξ,η)為(ξ,η)的質(zhì)量；g為重力加速度。

靜氣動彈性計(jì)算流程如圖5所示。

本文所采用的氣彈計(jì)算算例，是跨音速航空標(biāo)模M6機(jī)翼的靜氣動彈性變形計(jì)算。計(jì)算條件：馬赫數(shù)Ma=0.84，初始攻角α0=3.06°，飛行高度20 000m。

靜氣動彈性的計(jì)算收斂結(jié)果如圖6～7所示。從圖6中可以看出，M6機(jī)翼剛軸沿展向均有由于靜氣動彈性所引起的垂向形變，這主要是由于機(jī)翼的升、阻力聯(lián)合作用使得M6機(jī)翼產(chǎn)生了向上彎曲(根部固連于對稱面)；從圖7可以看出，由于氣動彈性的作用，機(jī)翼在展向各個位置均產(chǎn)生了逐漸增大的負(fù)扭轉(zhuǎn)角，這一負(fù)扭轉(zhuǎn)角使展向各位置剖面的當(dāng)?shù)赜菧p小，如圖8和圖9所示，機(jī)翼的升力分布在機(jī)翼產(chǎn)生靜氣動彈性變形后變化很大，同時機(jī)翼表面特別是機(jī)翼前、后緣的壓力分布也發(fā)生了較大變化，彈性機(jī)翼相對于剛性機(jī)翼來說，展向各剖面的升力系數(shù)均有顯著下降(主要由于機(jī)翼在氣動力作用下發(fā)生彎曲和扭轉(zhuǎn)變形，從而產(chǎn)生負(fù)的附加攻角)。靜氣動彈性變形給M6彈性機(jī)翼的跨音速升力特性帶來了負(fù)面的影響。

圖中，X/C為歸一化當(dāng)?shù)匾硇拖议L，Y/b為歸一化機(jī)翼展長。

3 并行化計(jì)算

通常航空氣動彈性問題的計(jì)算，效率都比較低，需要耗費(fèi)大量的計(jì)算時間和占用大量CPU[7]。本文使用MPI并行編程環(huán)境對氣動彈性解算程序進(jìn)行并行化程序改造，并使用4臺單核8線程IBM工作站和8臺單核8線程IBM工作站進(jìn)行集群化并行計(jì)算，與8核CPU的IBM刀片服務(wù)器的并行計(jì)算精度以及并行效率對比，對比結(jié)果見表1。

其中， 定義網(wǎng)格并行加速比如下：假設(shè)tseq是使用一個求解單元計(jì)算所需的時間，tp是用p個求解單元計(jì)算所需時間；并行加速比定義為Sp=tseq/tp，并行效率定義為Ep=Sp/p。表1中的時間為計(jì)算程序進(jìn)行一步時間迭代所需要的時間，單位為s。

由表1中可以看出，想要使用并行編程來獲得較高的并行加速比和并行效率，需要從以下幾個方面考慮：

a.采用良好的并行算法，在編程中最大限度地增加程序并行度，比如說采用顯式方法來替代隱式方法。

b.減少通信規(guī)模，提高通信帶寬，根據(jù)表1中的數(shù)據(jù)可以看出，集群計(jì)算的網(wǎng)絡(luò)通訊速率要小于多核服務(wù)器內(nèi)部系統(tǒng)總線的傳輸速率，從而使得網(wǎng)絡(luò)通信成為集群并行計(jì)算實(shí)現(xiàn)高效計(jì)算的瓶頸。

4 結(jié)束語

本文所采用的Delaunay圖映射的動彈性網(wǎng)格方法，能夠大幅度提高動彈網(wǎng)格的變形能力，并在變形過程中較好地保證了網(wǎng)格質(zhì)量和計(jì)算效率。本文通過對航空標(biāo)模M6機(jī)翼的靜氣動彈性計(jì)算可以看出，該方法能夠較好解決航空氣動彈性計(jì)算中關(guān)鍵的網(wǎng)格變形問題，同時耦合氣動力方程和靜平衡方程的流固耦合結(jié)算方法也具有較高的結(jié)算精度和較好的適應(yīng)性。本文所開發(fā)的并行化靜氣彈求解程序，能夠較好地解決航空復(fù)雜構(gòu)型的靜氣動彈性問題，在靜氣動彈性耦合數(shù)值模擬和氣彈剪裁等方面，也擁有較為廣闊的應(yīng)用前景。

[1] Blom F J. Considerations on the spring analogy[J]．Journal of Aircraft, 2000,32: 647-668.

[2] Liu F, Cai J, Zhu Y．Calculation of wing flutter by a coupled CFD-CSD method [R]．AIAA-2000-0907,2000.

[3] 范銳軍，馮朝輝．大展弦比無人機(jī)的靜氣彈問題計(jì)算及分析[J]．力學(xué)季刊， 2009(4)：182-188．

[4] 沈克揚(yáng)，張錫華．彈性機(jī)翼的跨音速壓強(qiáng)分布計(jì)算[J]．空氣動力學(xué)學(xué)報，1984，1(3)：221～230．

[5] 王軍利．基于非結(jié)構(gòu)動網(wǎng)格的非定常氣動力計(jì)算[D]．西安：西北工業(yè)大學(xué)，2004．

[6] 史愛明．非結(jié)構(gòu)動網(wǎng)格下非定常氣動力計(jì)算和嗡鳴研究[D]．西安：西北工業(yè)大學(xué)，2003．

[7] 范銳軍，白俊強(qiáng)．氣動彈性及并行化求解研究[J]．彈箭與制導(dǎo)學(xué)報，2006，6(1)：95-102．

Application of improved dynamic unstructured grids in aeroelastic model

HU Fan1, FAN Ruijun2

(1.Shaanxi Science and Technology Resource Center, Shaanxi Xi'an, 710077, China)

(2.Xi'an Modern Control Technology Research Institute, Shaanxi Xi'an, 710065, China)

Aeroelastic model is significant for large amount of airplanes in modern aerodynamics computing. This paper presents a strategy for generating 3D unstructured grids and the dynamic grids based on Delaunay graph mapping method. On the base of above, it performs a set of flow field solver based on Euler equations into establishmeng of the static aeroelastic cases of the M6 standard model coupled with structure dynamic equation, shows the design of the parallel computing method and numerical example.

dynamics grids deformation technique based on Delaunay graph mapping dynamics unstructured grids; Euler equations; CFD/CSD coupled technology; static aeroelastics

10.3969/j.issn.2095-509X.2015.03.013

2015-02-10

國家自然科學(xué)基金/青年基金項(xiàng)目資助(11302175)

胡凡(1975—)，女，陜西咸陽人，陜西省科技資源統(tǒng)籌中心工程師，主要研究方向?yàn)楣I(yè)設(shè)計(jì)。

V211.3

A

2095-509X(2015)03-0060-05