柳兆偉，侯中喜，陳立立

適用于臨近空間飛行器大變形的動網(wǎng)格策略*

柳兆偉，侯中喜，陳立立

(國防科技大學(xué) 航天科學(xué)與工程學(xué)院， 湖南 長沙 410073)

對于超大展弦比構(gòu)型的低速臨近空間飛行器而言，由于其在飛行過程中結(jié)構(gòu)變形非常顯著，因此基于計算流體力學(xué)的分析方法對于動網(wǎng)格提出了非常高的要求。為此，提出了一種適用于邊界大變形的動網(wǎng)格策略，該種動網(wǎng)格基于映射的思想，將邊界網(wǎng)格的位置變化以某種權(quán)重反映到流場網(wǎng)格，并更新網(wǎng)格節(jié)點位置。選取距離倒數(shù)的n次方作為權(quán)重，研究不同的權(quán)重指數(shù)n對網(wǎng)格變形的影響規(guī)律，然后開展了二維與三維動網(wǎng)格實例分析。結(jié)果表明，這種動網(wǎng)格方法能夠很好地適用于大變形的情形，并能很好地保證變形后的網(wǎng)格質(zhì)量。

動網(wǎng)格；大變形；變形策略

(College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China)

在軍用和民用領(lǐng)域巨大需求的牽引下，高空長航時(High Altitude Long Endurance,HALE)飛行器得到快速的發(fā)展，特別的以“太陽神”[1]“微風(fēng)”“陽光動力”等為代表的一系列太陽能飛機(jī)的發(fā)展，大大促進(jìn)了該技術(shù)的提升。如圖1所示，展示了幾種高空長航時太陽能飛行器。為了實現(xiàn)高空長航時這一目標(biāo)，該類型飛行器結(jié)構(gòu)面密度通常較低，由此導(dǎo)致剛度明顯不足，在飛行過程中受到氣動載荷時，結(jié)構(gòu)變形非常顯著，同時結(jié)構(gòu)變形又使得氣動性能發(fā)生改變，因而結(jié)構(gòu)與氣動出現(xiàn)較強(qiáng)耦合。

圖1 高空長航時大展弦比飛行器Fig1 HALE high-aspect-ratio aircrafts

近年來，基于計算流體力學(xué)的流固耦合分析方法得到快速的發(fā)展與廣泛的應(yīng)用[2-6]。這種方法要求結(jié)構(gòu)和氣動兩個學(xué)科獨立建模，并采用交錯求解的方式進(jìn)行計算。流體計算中一般采用基于空間位置的Euler網(wǎng)格，因而在耦合過程中，流體網(wǎng)格需要根據(jù)結(jié)構(gòu)邊界的移動而變化。特別的，對于超大展弦比構(gòu)型的低速臨近空間飛行器而言，在飛行過程中其結(jié)構(gòu)變形非常顯著，因而需要發(fā)展適用于大變形的動網(wǎng)格策略。

目前，動網(wǎng)格主要的實現(xiàn)方法可分為兩大類[7]：網(wǎng)格重構(gòu)和網(wǎng)格變形。網(wǎng)格重構(gòu)技術(shù)基于超限插值(TransFinite Interpolation,TFI)技術(shù)，在每次邊界變化時，重新劃分網(wǎng)格，這種方法計算量較大，效率較低。另一種常用的動網(wǎng)格方式是網(wǎng)格變形技術(shù)，該方法將網(wǎng)格作為彈簧或彈性體來處理，并根據(jù)靜力平衡計算得到新的網(wǎng)格點位置。但是彈簧算法在處理大變形或者較密的網(wǎng)格時，常會出現(xiàn)交叉而產(chǎn)生負(fù)體積網(wǎng)格，導(dǎo)致網(wǎng)格更新失敗。近年來， Liu[7]等采用Delaunay背景網(wǎng)格的變形方法實現(xiàn)網(wǎng)格更新，這種動網(wǎng)格方法能夠適應(yīng)于大變形的情形，特別是對于大的位移情況，然而當(dāng)有大的扭轉(zhuǎn)變形時，將會出現(xiàn)Delaunay圖的交叉，而導(dǎo)致網(wǎng)格更新失敗。

本文提出一種基于映射的網(wǎng)格變形策略，其基本思想是：不改變網(wǎng)格的拓?fù)浣Y(jié)構(gòu)，將邊界的變形量按照一定的權(quán)重映射到流域中的網(wǎng)格點，從而確定流場中網(wǎng)格的位移。選取待移動網(wǎng)格節(jié)點到邊界上點的距離倒數(shù)的n次方作為權(quán)重，其出發(fā)點是，要調(diào)整的網(wǎng)格點到邊界網(wǎng)格點距離越近，受到邊界網(wǎng)格點的影響也越大，通過調(diào)整指數(shù)n可以調(diào)整影響的擴(kuò)散范圍。

1 網(wǎng)格變形

1.1 網(wǎng)格變形方法

網(wǎng)格更新方法可分為四個步驟：

第一步，計算網(wǎng)格節(jié)點至邊界點的距離。設(shè)計算域為D，節(jié)點個數(shù)為p，邊界區(qū)域為B，節(jié)點個數(shù)為q，其中D包括B。則計算域D內(nèi)任意一點i到邊界B上一點j的距離為dij。且：

(1)

第二步，計算邊界網(wǎng)格點位移量。按照一定的要求變化邊界區(qū)域B，得到網(wǎng)格邊界B上任意一點j的位移量：

(2)

第三步，選取距離的倒數(shù)的n次方作為網(wǎng)格更新的權(quán)重，計算網(wǎng)格位移量。對于區(qū)域D內(nèi)的任意一點i而言，根據(jù)選取的n，按照映射關(guān)系，可確定其位移量為：

(3)

有一點需要特別注意的是，當(dāng)點為邊界上的點時，其位移變化則不受到其他點的影響。即若dij=0，則該點的位移為：

Δri=Δrj

(4)

第四步，根據(jù)上述位移量，重新確定網(wǎng)格點坐標(biāo)。對于計算域D內(nèi)的任意一點i，其新的坐標(biāo)位置為：

(5)

1.2 網(wǎng)格更新的實現(xiàn)

整個網(wǎng)格變形的思想是：將邊界的運動按照某種權(quán)重映射到每個網(wǎng)格點上。因而在實現(xiàn)過程中，可分為內(nèi)外層兩個循環(huán)，內(nèi)層是計算邊界位移對于某一網(wǎng)格節(jié)點的權(quán)重，得到網(wǎng)格點位移；外層則是循環(huán)所有網(wǎng)格節(jié)點，更新坐標(biāo)位置。

具體過程可寫為如下的偽代碼：

for i=1：number of domain points

{

Δrsum=0;

dback=0;

for j=1：number of boundary points

{

更新邊界點j坐標(biāo)，得到Δrj;

計算該點到邊界點j的距離dij;

if dij= 0

{

Δrsum=Δrj；

dback=1；

break；

}

}

end

}

end

2 權(quán)重指數(shù)的影響

選取(1/d)n作為權(quán)重，出發(fā)點是：要調(diào)整的網(wǎng)格點到邊界網(wǎng)格點距離越近，受到邊界網(wǎng)格點的影響也越大，通過調(diào)整指數(shù)n可以調(diào)整影響的擴(kuò)散范圍。

下面結(jié)合二維大變形的動網(wǎng)格算例，通過觀察不同的n值得到網(wǎng)格的差異，研究n的取值對于網(wǎng)格變形的影響規(guī)律。設(shè)原始網(wǎng)格如圖2所示，將大變形分為內(nèi)部物體的扭轉(zhuǎn)變形和平移變形兩種情況。則不同的n值對應(yīng)的變形網(wǎng)格如圖3和圖4所示。

圖2 原始網(wǎng)格Fig.2 The initial mesh

(a)n=2時網(wǎng)格(a) The mesh while n=2(b) n=4時的網(wǎng)格(b) The mesh while n=4

(c) n=6的網(wǎng)格(c) The mesh while n=6(d) n=8的網(wǎng)格(d) The mesh while n=8圖3 物體轉(zhuǎn)動時不同n值對應(yīng)的物體周圍網(wǎng)格Fig.3 The mesh versus different n with torsion deformation

對于純扭轉(zhuǎn)位移而言，假設(shè)扭轉(zhuǎn)角度為60°。當(dāng)n=2時變形網(wǎng)格如圖3(a)所示，可以看出，外層網(wǎng)格變形較小，而內(nèi)層網(wǎng)格扭曲較為嚴(yán)重，出現(xiàn)負(fù)體積網(wǎng)格。而隨著n值的增大，內(nèi)層網(wǎng)格對于變形壁面的跟隨性也越好，外層網(wǎng)格扭轉(zhuǎn)幅度增大，這也意味著壁面網(wǎng)格質(zhì)量更好，如圖3(c)、(d)。

平移變形情況下，不同的n對應(yīng)的網(wǎng)格變形結(jié)果如圖4所示?？梢钥闯觯琻取值越大，變形物體周圍內(nèi)層網(wǎng)格的變形越小，邊界的平移變形被傳播到更遠(yuǎn)的外層網(wǎng)格區(qū)域。然而當(dāng)外層網(wǎng)格需要承受非常大的變形時，則可能會出現(xiàn)交叉而使得網(wǎng)格更新失敗。

(a) n=2時網(wǎng)格(a) The mesh while n=2

(b) n=3時網(wǎng)格(b) The mesh while n=3

(c) n=4時網(wǎng)格(c) The mesh while n=4圖4 物體平移時不同n值對應(yīng)的物體周圍網(wǎng)格 Fig.4 The mesh versus different n with large displacement

通過設(shè)定不同的n值，觀察網(wǎng)格變形的規(guī)律，可以發(fā)現(xiàn)：n值越大，運動邊界周圍網(wǎng)格的剛性越強(qiáng)，對于扭轉(zhuǎn)變形的適應(yīng)能力也越強(qiáng)；而對于平移變形而言，當(dāng)n較大時，運動壁面的位移傳播到外層網(wǎng)格中，而使得遠(yuǎn)離壁面的區(qū)域的網(wǎng)格變形較大。綜上，n值的選取可采取如下的原則：一般n取2～6，且當(dāng)扭轉(zhuǎn)變形較大時，可取相對較大的n值，而當(dāng)平移較大時，可適當(dāng)取較小的n值。同時，對于初始網(wǎng)格本身而言，更大的流域以及較大的外部網(wǎng)格尺寸也會增大網(wǎng)格變形的空間，增強(qiáng)變形能力。

3 網(wǎng)格變形實例

下面結(jié)合不同的大變形動網(wǎng)格實例，分析變形后網(wǎng)格的質(zhì)量。算例包括二維翼型的旋轉(zhuǎn)與平移、二維不規(guī)則變形、三維球體的變形與移動以及三維機(jī)翼大幅縱向變形。

3.1 二維翼型網(wǎng)格變形

在飛行器的性能分析中，對典型截面的非線性氣動彈性研究具有重要的意義?；谟嬎懔黧w力學(xué)的分析手段能夠非常好地預(yù)測流體的分離等非線性現(xiàn)象等。下面考慮如下的翼型C形結(jié)構(gòu)化網(wǎng)格，如圖5所示，當(dāng)翼型具有大的平移和扭轉(zhuǎn)時，根據(jù)上述方法得到的網(wǎng)格分別如圖6、圖7、圖8所示。

圖5 翼型原始網(wǎng)格Fig.5 The initial mesh of airfoil

圖6 翼型扭轉(zhuǎn)60°網(wǎng)格變形圖(n=4)Fig.6 The mesh with 60° torsion deformation (n=4)

圖7 大范圍平移后網(wǎng)格變形圖(n=2)Fig.7 The mesh with large displacement (n=2)

圖8 翼型同時扭轉(zhuǎn)和平移后的變形網(wǎng)格(n=2)Fig.8 The mesh with large displacement and torsion (n=2)

由上述變形后的網(wǎng)格可以看出，在翼型出現(xiàn)非常大的扭轉(zhuǎn)或者平移時，可以通過調(diào)節(jié)權(quán)重指數(shù)n的值，將翼型邊界的位移很好的傳播到大尺寸的網(wǎng)格中，從而保證良好的壁面邊界層網(wǎng)格質(zhì)量。

3.2 二維壁面不規(guī)則變形

在飛行器外形設(shè)計及翼型優(yōu)化等方面，可能會涉及到曲線或曲面的不規(guī)則變形等情形。針對如圖9所示的網(wǎng)格區(qū)域，假設(shè)內(nèi)部邊界的上表面出現(xiàn)正弦波動變形，而下表面出現(xiàn)鋸齒狀變形。

圖9 初始網(wǎng)格以及邊界變化Fig.9 The initial mesh and the deformation of the wall

觀察上述邊界變化，可知這種變形主要是小區(qū)域內(nèi)的大扭轉(zhuǎn)變形，此時我們需要保證靠近壁面的網(wǎng)格質(zhì)量，并避免出現(xiàn)過度扭曲而導(dǎo)致網(wǎng)格更新失敗，根據(jù)上面n值的選擇原則，應(yīng)取較大的n值，以保證壁面附近網(wǎng)格的剛度，此處取n=6，如圖10所示為更新后的網(wǎng)格。

圖10 更新后的網(wǎng)格(n=6)Fig.10 The mesh after deformation (n=6)

對于上表面而言，邊界扭曲較為嚴(yán)重，網(wǎng)格受到一定的拉伸或壓縮；而對于下表面而言，網(wǎng)格變形后依然較為均勻；而側(cè)面由于未變形，網(wǎng)格變化很小?？梢钥闯?，當(dāng)壁面出現(xiàn)不規(guī)則變形時，采用這種動網(wǎng)格方法，通過選取合適的n值，能夠較好地保證變形后的網(wǎng)格質(zhì)量。

3.3 三維球體網(wǎng)格變形

對于臨近空間球形或者艇形的浮空器而言，當(dāng)內(nèi)外壓差變化或者受到外界氣流干擾時，其結(jié)構(gòu)會出現(xiàn)大幅變形。下面以三維球體的外流場計算網(wǎng)格為例，驗證上述動網(wǎng)格方法。

球體的原始網(wǎng)格如圖11所示，在球體的周圍網(wǎng)格較密，而外層網(wǎng)格較為稀疏。圖12(a)、(b)、(c)分別為球體發(fā)生平移、變形、扭轉(zhuǎn)后的變形網(wǎng)格結(jié)果?？梢钥闯觯捎蒙鲜鼍W(wǎng)格變形方法，能夠獲得高質(zhì)量的網(wǎng)格。

圖11 球體原始網(wǎng)格Fig.11 The initial mesh of sphere

(c) 內(nèi)部球體旋轉(zhuǎn)60°網(wǎng)格圖(n=6)(c) The mesh with 60° rotation of inside sphere (n=6)圖12 不同變形情形下的網(wǎng)格結(jié)果Fig.12 The mesh with different deformation case

3.4 三維機(jī)翼網(wǎng)格變形

對于臨近空間大展弦比長航時太陽能飛機(jī)而言，為了降低能源消耗，其重量通常相對較小，由此導(dǎo)致結(jié)構(gòu)剛度不足，在飛行過程中常常會出現(xiàn)非常大的縱向變形。下面假定一定幅度的機(jī)翼縱向變形，查看上述動網(wǎng)格方法對于變形的適應(yīng)性及變形后網(wǎng)格的效果。

在機(jī)翼的氣動力計算中，通常機(jī)翼周圍的網(wǎng)格較密(以提高計算精度)，而遠(yuǎn)離機(jī)翼壁面的區(qū)域網(wǎng)格非常稀疏，如圖13所示。

圖13 機(jī)翼原始網(wǎng)格Fig.13 The initial mesh of wing

希望在網(wǎng)格變形后，依然能夠保持良好的壁面網(wǎng)格質(zhì)量。假設(shè)機(jī)翼的翼尖縱向變形量達(dá)到展長的50%，如圖14所示。采用本文的網(wǎng)格變形策略得到變形后的網(wǎng)格如圖15、圖16所示。可以看出，大的縱向變形出現(xiàn)后，機(jī)翼壁面周圍網(wǎng)格質(zhì)量依然較好，而變形被傳播到遠(yuǎn)離機(jī)翼壁面的外層大尺度網(wǎng)格中，這對氣動計算而言非常重要。

圖14 機(jī)翼縱向變形Fig.14 The longitudinal deformation of the wing

圖15 不同截面上的網(wǎng)格(n=4)Fig.15 The mesh at different sections (n=4)

圖16 變形后的機(jī)翼網(wǎng)格(n=4)Fig16 The wing mesh after deformation (n=4)

4 結(jié)論

基于映射的思想，提出一種網(wǎng)格變形策略，其基本思想是：將邊界的變形量以某種權(quán)重反映到流場區(qū)域網(wǎng)格中。文中以(1/d)n為權(quán)重，并結(jié)合二維算例研究了指數(shù)n對于網(wǎng)格變形的影響，結(jié)果表明，指數(shù)n值越大，邊界周圍網(wǎng)格的剛性越強(qiáng)，而變形更容易被傳播到遠(yuǎn)離邊界的區(qū)域。n值的選取可采取如下的原則：一般取2～6，當(dāng)扭轉(zhuǎn)變形較大時，可取相對較大的n值，而當(dāng)平移較大時，可適當(dāng)取較小的n值。

References)

[1]NollTE,BrownJM,Perez-DavisME,etal.Investigationoftheheliosprototypeaircraftmishap,volumeI:mishapreport[R].USA:NationalOceanicandAtmosphericAdministration, 2004.

[2]YangGW,ChenDW,CuiK.Responsesurfacetechniqueforstaticaeroelasticoptimizationonahigh-aspect-ratiowing[J].JournalofAircraft, 2009, 46(4): 1444-1450.

[3]PalaciosR,CesnikCES.Staticnonlinearaeroelasticityofflexibleslenderwingsincompressibleflow[C]//Proceedingsof46thAIAA/ASME/ASCE/AHS/ASCStructures,StructuralDynamicsandMaterialsConference, 2005.

[4]HallissyBP,CesnikCES.High-fidelityaeroelasticanalysisofveryflexibleaircraft[C]//Proceedingsof52ndAIAA/ASME/ASCE/AHS/ASCStructures,StructuralDynamicsandMaterialsConference, 2011.

[5]GarciaJA.Numericalinvestigationofnonlinearaeroelasticeffectsonflexiblehigh-aspect-ratiowings[J].JournalofAircraft, 2005, 42(4): 1025-1036.

[6]CarnieG,QinN.Fluid-structureinteractionofHALEwingconfigurationwithanefficientmovinggridmethod[C]//Proceedingsof46thAIAAAerospaceSciencesMeetingandExhibit, 2008.

[7]LiuXQ,QinN,XiaH.FastdynamicgriddeformationbasedonDelaunaygraphmapping[J].JournalofComputationalPhysics, 2006, 211(2): 405-423.

Moving mesh strategy for large deformation of near-space aircrafts

LIU Zhaowei, HOU Zhongxi，CHEN Lili

The high-aspect-ratio low-speed near-space aircrafts may undergo very large deformation during flight, so a high demand of moving mesh is required for the analysis method based on computational fluid dynamics. To this end, a moving mesh strategy for large deformation of the boundary was presented. The strategy which is based on the mapping interpolation method reflects the displacement of boundary mesh to flow field mesh using a certain kind of weight and then updates the position of mesh nodes. Inverse distance’snth-power was chosen as the weighting factor and the influence of different weight indexnon the mesh deformation was studied, then the analysis of some two-dimensional and three-dimensional moving mesh cases was carried out. The results suggest that this method is capable of handling the large deformation and ensuring the quality of deformed mesh.

moving mesh; large deformation; deformation strategy

2015-04-08

國家高分重大專項資助項目(GFZX04060103)

柳兆偉(1988—)，男，安徽臨泉人，博士研究生，E-mail：liuzhaowei@nudt.edu.cn；侯中喜(通信作者)，男，教授，博士，博士生導(dǎo)師，E-mail：hzx@nudt.edu.cn

10.11887/j.cn.201504004

http://journal.nudt.edu.cn

V211.3

A

1001-2486(2015)04-019-06