彈簧近似法在二維非結(jié)構(gòu)動網(wǎng)格生成技術(shù)中的應(yīng)用

霍世慧， 王富生， 岳珠峰

（西北工業(yè)大學(xué) 力學(xué)與土木建筑學(xué)院，西安 710129）

在工程實際問題中，存在著很多包括運動邊界的非定常復(fù)雜流場，如自由液面流動問題、流體與結(jié)構(gòu)耦合問題、氣動彈性問題、強迫振動及多體分離問題等，這些問題中運動邊界的處理是計算流體力學(xué)的一個難題。在這些問題的求解過程中，計算區(qū)域隨著時間的變化而變化，相應(yīng)區(qū)域網(wǎng)格也需發(fā)生改變。在網(wǎng)格的變形中，常用到的有網(wǎng)格重生成和網(wǎng)格變形技術(shù)。網(wǎng)格重生成技術(shù)是根據(jù)結(jié)構(gòu)變形修改幾何模型，并全部或部分重新生成模型網(wǎng)格；網(wǎng)格變形技術(shù)是基于原有網(wǎng)格進行網(wǎng)格位置、形狀的修改生成新的網(wǎng)格，而網(wǎng)格的拓?fù)浣Y(jié)構(gòu)沒有發(fā)生改變。網(wǎng)格重生產(chǎn)技術(shù)生成模型網(wǎng)格的過程復(fù)雜，耗費大量的時間，而且在模型局部特殊部位需要手動修改才能生成較好的網(wǎng)格，限制了動網(wǎng)格的應(yīng)用，而網(wǎng)格變形技術(shù)只是對網(wǎng)格進行位置、形狀的修改，很好地克服了網(wǎng)格重生成技術(shù)中的這些缺點。本文將開展網(wǎng)格變形技術(shù)在動網(wǎng)格中的應(yīng)用研究。

目前國內(nèi)外發(fā)展的可靠性較高的網(wǎng)格變形技術(shù)主要有彈性體近似法［1－3］（Elastic solid method，ESM）和彈簧近似法［4－8］（Spring analogy method，SAM）。彈性體近似法是把整個計算區(qū)域看成一個彈性體，通過求解控制網(wǎng)格點移動的基本方程來求出網(wǎng)格內(nèi)部點的位移；彈簧近似法是把整個計算區(qū)域看成一個由彈簧組成的系統(tǒng)，當(dāng)邊界發(fā)生運動時，內(nèi)部網(wǎng)格點按當(dāng)?shù)氐倪\動強度重新布置達到新的平衡。

國內(nèi)外已經(jīng)開展了很多基于彈簧近似法的研究，Batina［4］提出了彈簧近似法在三角形非結(jié)構(gòu)動網(wǎng)格中的應(yīng)用，并將其應(yīng)用到機翼顫振問題的分析研究中；Farhat［5］在彈簧近似法中引入了扭轉(zhuǎn)系數(shù)，在一定程度上避免了網(wǎng)格邊的交叉，擴大了彈簧近似法的應(yīng)用范圍；Blom［6］定義了彈簧近似法的兩種描述方法及其倔強系數(shù)的選取，并研究了其在二維可動邊界問題中的應(yīng)用；Mitsuhiro［7］實現(xiàn)了彈簧近似法在三維動網(wǎng)格中的應(yīng)用，并提出了根據(jù)單元形狀對倔強系數(shù)的修正；Zhang［8］將彈簧近似法和局部網(wǎng)格重生成聯(lián)合運用到動網(wǎng)格中，提高了網(wǎng)格的變形能力。上述文獻分別從扭轉(zhuǎn)系數(shù)、單元形狀和局部網(wǎng)格重生成等幾個方面對彈簧近似法進行了一定的改進，并將其應(yīng)用到實際問題的研究中，得到了很多有用的結(jié)論，但都沒有對變形后的網(wǎng)格質(zhì)量進行系統(tǒng)的分析。本文將開展變形后網(wǎng)格質(zhì)量的研究，并討論將網(wǎng)格質(zhì)量引入倔強系數(shù)后的網(wǎng)格變形能力。

1 彈簧近似法

彈簧近似法是將網(wǎng)格單元的各條棱邊假設(shè)為彈簧，邊界發(fā)生運動后，內(nèi)部網(wǎng)格由彈簧控制達到新的平衡，從而生成變形后網(wǎng)格。彈簧近似方法有兩種描述方法［6］：一種是將彈簧的平衡長度假設(shè)為零，稱為頂點彈簧（Vertex springs）；另一種是將彈簧的平衡長度假設(shè)為其變形前長度，稱為棱邊彈簧（Segment springs）。頂點彈簧近似法能夠用于網(wǎng)格變形和網(wǎng)格優(yōu)化中，網(wǎng)格變形中認(rèn)為網(wǎng)格節(jié)點的受力始終等于初始狀態(tài)所受的合力；網(wǎng)格優(yōu)化中則認(rèn)為網(wǎng)格節(jié)點始終處于平衡狀態(tài)。本文主要研究頂點彈簧近似法在網(wǎng)格變形中的應(yīng)用。

在頂點彈簧近似法中，節(jié)點i，j間的彈簧張力可以表示為式（1）：

其中：Kij為連接節(jié)點i，j彈簧的倔強系數(shù)；ri、rj分別為節(jié)點i，j的位置矢量。假設(shè)計算區(qū)域中共有N個節(jié)點與節(jié)點i相連，則其合力可表示為：

由式（2）可得整個計算區(qū)域中m個節(jié)點在初始狀態(tài)下的矩陣表達式為：

當(dāng)模型邊界發(fā)生變形時，通過改變模型邊界的位置矢量和增加適當(dāng)?shù)耐膺吔缦拗茥l件，如計算區(qū)域外邊界位置矢量保持不變、模型邊界位置矢量改為新位置下的矢量等，運用Gauss－seidel迭代法對式（3）進行迭代求解，其相應(yīng)的分量迭代格式可表示式為：

當(dāng)模型發(fā)生變形時，利用式（4）對式（3）進行數(shù)次迭代，達到所需要的精度時，便可得到計算區(qū)域中各網(wǎng)格節(jié)點變形后坐標(biāo)。

2 彈簧倔強系數(shù)的選取

在上述彈簧近似法的描述中，各單元棱邊都被假設(shè)為一彈簧，都具有各自的倔強系數(shù)K，下面將進行彈簧倔強系數(shù)選取的研究。

2.1 傳統(tǒng)方法

在傳統(tǒng)彈簧倔強系數(shù)的選取中一般都只把網(wǎng)格邊的邊長作為參考值來定義［6］。彈簧倔強系數(shù)選取的一般格式為：

式中：lij為網(wǎng)格邊的邊長，當(dāng)lij→0時，Kij→∞，這樣較好地避免了網(wǎng)格節(jié)點i與j發(fā)生相撞，使網(wǎng)格取得較好的疏密特征，但在處理一些邊界發(fā)生較大運動或變形問題時，容易產(chǎn)生體積為負(fù)的無效單元。近年來，許多文獻中也提出了對傳統(tǒng)彈簧倔強系數(shù)的改進方法。史愛明［9］在傳統(tǒng)彈簧倔強系數(shù)中引入了對單元幾何尺寸的控制，每個時間步更新彈簧倔強系數(shù)，從而保證體積最小單元不會成為體積為負(fù)的無效單元，取得了較好的效果，其具體表達式為：

其中，Vmin為擁有網(wǎng)格邊ij單元的最小體積。Farhat［5］在每個網(wǎng)格的頂點上加了一個扭轉(zhuǎn)彈簧，每一個連接在頂點i的三角形單元扭轉(zhuǎn)彈簧的扭轉(zhuǎn)剛度Ci如式（7）所示：

其中，θk為以i為頂點三角形面單元的內(nèi)角。當(dāng)邊界發(fā)生運動或變形時，對每個頂點均可得到關(guān)于拉伸彈簧的力平衡方程和扭轉(zhuǎn)彈簧的力矩平衡方程，聯(lián)立求解得出計算區(qū)域內(nèi)的新網(wǎng)格點。

2.2 方法改進

網(wǎng)格的質(zhì)量隨著邊界運動或變形發(fā)生變化，下面將開展網(wǎng)格質(zhì)量的研究并提出相應(yīng)的改進方案。圖1給出了在二維頂點彈簧近似法中，三角網(wǎng)格模型容易產(chǎn)生的狹長、扁平兩種典型畸形網(wǎng)格單元。

單元的這兩種畸形情況嚴(yán)重影響著網(wǎng)格整體質(zhì)量，并最終導(dǎo)致網(wǎng)格邊的交叉。為了使動網(wǎng)格達到最大限度的變形情況、推遲畸形網(wǎng)格的出現(xiàn)，就必須在彈簧倔強系數(shù)中引入對網(wǎng)格質(zhì)量的控制。對網(wǎng)格質(zhì)量的評定，主要有傾斜度和縱橫比兩個指標(biāo)［10］，其計算公式分別為：

圖1 典型畸形網(wǎng)格單元Fig.1 Classical abnormal meshes

在網(wǎng)格的變形過程中，通過兩個網(wǎng)格質(zhì)量評定指標(biāo)的修正，逐漸增大網(wǎng)格質(zhì)量較差單元各棱邊的倔強系數(shù)，使網(wǎng)格質(zhì)量較好單元承擔(dān)一定的變形。

3 方法應(yīng)用

在進行彈簧近似法的應(yīng)用之前，首先引入兩種評定網(wǎng)格變形指標(biāo)［11］：單元面積變化和單元形狀變化。通過這些指標(biāo)來量化網(wǎng)格的變形，分析網(wǎng)格的變形程度。其具體評價方式為：

本文在算例一中采用NACA23012翼型，對其計算區(qū)域采用非結(jié)構(gòu)化離散，圖2（a）為翼型初始二維氣動網(wǎng)格，該網(wǎng)格是由非結(jié)構(gòu)三角形單元組成，共計有67 160個節(jié)點、133 966個單元，圖2（b）為翼型表面網(wǎng)格局部圖。

圖2 NACA23012初始氣動網(wǎng)格Fig.2 Original mesh of NACA23012

圖3 網(wǎng)格變形指標(biāo)變化情況Fig.3 The change of mesh deformation value

機翼的變形表現(xiàn)在二維翼型中為翼型的旋轉(zhuǎn)，下面將針對NACA23012翼型通過軟件Fortran編寫傳統(tǒng)和改進的彈簧近似法程序生成其不同旋轉(zhuǎn)角度時網(wǎng)格。為了清楚地比較兩程序生成網(wǎng)格中各單元的變形情況，圖3給出了在發(fā)生相同旋轉(zhuǎn)角度時，兩方法生成翼型周圍各網(wǎng)格的變形指標(biāo)值。圖3中區(qū)域顏色的深淺直接表示該處網(wǎng)格變形指標(biāo)值的大小，隨著顏色的加深值增大，區(qū)域網(wǎng)格變形加劇。

從圖3可以明顯地看到，傳統(tǒng)彈簧近似法在翼型前緣和后緣處網(wǎng)格變形指標(biāo)較大，而其它區(qū)域網(wǎng)格變形指標(biāo)較小，網(wǎng)格的變形主要集中在翼型前后緣處；改進彈簧近似法生成網(wǎng)格的變形指標(biāo)呈現(xiàn)以前緣、后緣點為中心向外擴展的形式，網(wǎng)格的變形較為分散，緩解了翼型前后緣的變形壓力。

圖4給出了傳統(tǒng)和改進彈簧近似法下，計算區(qū)域網(wǎng)格整體變形指標(biāo)fA、fAR隨翼型旋轉(zhuǎn)角度的變化情況。

圖4 變形指標(biāo)隨旋轉(zhuǎn)角度變化曲線Fig.4 Deformation value with the increase of rotation angle

從圖4可以發(fā)現(xiàn)，隨著翼型旋轉(zhuǎn)角度的增加，fA、fAR均逐漸增大，網(wǎng)格變形隨旋轉(zhuǎn)角度的增加而增大。比較圖4中相同旋轉(zhuǎn)角度兩種彈簧近似法生成網(wǎng)格的變形指標(biāo)可以發(fā)現(xiàn)，改進的彈簧近似法fA、fAR值始終低于傳統(tǒng)方法，改進的彈簧近似法將邊界的變形壓力分散到了整個計算區(qū)域，降低了整體網(wǎng)格的變形指標(biāo)值。

改進的彈簧近似法緩解了網(wǎng)格變形在計算區(qū)域局部集中出現(xiàn)的情況，將變形壓力分散到整個計算區(qū)域，下面將開展變形壓力分散前后，兩種方法生成網(wǎng)格的質(zhì)量變化情況的研究。圖5給出了翼型邊界發(fā)生旋轉(zhuǎn)時兩種彈簧近似法生成網(wǎng)格的整體網(wǎng)格質(zhì)量變化曲線。

從圖5中可以看出：傳統(tǒng)和改進的彈簧近似法生成網(wǎng)格的qs、qa均隨旋轉(zhuǎn)角度的增加而增大，網(wǎng)格質(zhì)量下降；發(fā)生相同旋轉(zhuǎn)角度時，改進的彈簧近似法生成的網(wǎng)格質(zhì)量均優(yōu)于傳統(tǒng)彈簧近似法；傳統(tǒng)方法在邊界發(fā)生13°旋轉(zhuǎn)角度時，qs、qa值達到1，網(wǎng)格質(zhì)量最差，無法繼續(xù)使用，改進方法在邊界發(fā)生27°旋轉(zhuǎn)角度時qs、qa值才達到1，通過方法的改進，提高了網(wǎng)格的變形能力。圖6、圖7分別給出了翼型發(fā)生14°旋轉(zhuǎn)時運用兩種彈簧近似法生成的網(wǎng)格局部圖，右邊小圖為翼型后緣局部網(wǎng)格放大圖。

圖5 網(wǎng)格質(zhì)量隨旋轉(zhuǎn)角度變化曲線Fig.5 Mesh quality with the increase of rotation angle

對比圖6、圖7中兩種彈簧近似法生成的翼型發(fā)生14°旋轉(zhuǎn)時網(wǎng)格，可以很明顯地看出，傳統(tǒng)彈簧近似法在后緣局部發(fā)生了網(wǎng)格邊的交叉，網(wǎng)格不能繼續(xù)使用，而改進彈簧近似法網(wǎng)格卻能夠保持較好的質(zhì)量。

4 算例

4.1 算例一

運用改進的彈簧近似法生成非結(jié)構(gòu)動網(wǎng)格，分析二維NACA23012翼型不同攻角下的升力系數(shù)變化情況NACA23012翼型初始網(wǎng)格如圖2所示，圖8給出了運用改進彈簧近似法生成的翼型在 －5°、0°、5°和10°攻角下網(wǎng)格局部圖。

運用上面生成的翼型網(wǎng)格，選定雷諾數(shù)為6.0×106，分析了機翼在不同攻角下的升力系數(shù)。圖9給出了機翼升力系數(shù)隨攻角的變化情況，并與參考文獻［12］所得到的實驗數(shù)據(jù)進行比較。

觀察圖9中曲線可以發(fā)現(xiàn)：雷諾數(shù)為6.0×106時，NACA23012翼型升力系數(shù)隨攻角的增大而提高，在攻角達到12°左右時，翼型發(fā)生失速，升力系數(shù)開始下降。通過將模擬結(jié)果與實驗數(shù)據(jù)比較發(fā)現(xiàn)，該模擬方法能夠較好地得出翼型的升力系數(shù)變化情況。圖10給出了翼型在 －5°、0°、5°和 10°攻角下翼型周圍壓力分布情況。

圖8 翼型周圍網(wǎng)格局部圖Fig.8 Close-up views of mesh around airfoil

圖9 升力系數(shù)隨攻角變化情況Fig.9 Effect of Attack Angle on Lift Coefficient

4.2 算例二

比較圖11、圖12可以發(fā)現(xiàn)，翼型升力系數(shù)的變化趨勢與翼型攻角的變化完全吻合，其周期均為T=2π/ω=51.1 s，這與實際情況較為吻合。圖13給出了升力系數(shù)與翼型攻角的變化曲線，并與傳統(tǒng)彈簧近似法及實驗數(shù)據(jù)進行比較。

由圖13可以發(fā)現(xiàn)，改進彈簧近似法生成的動網(wǎng)格在二維諧和振動翼型非定常氣動力的求解中能夠得到較傳統(tǒng)方法更接近于實驗數(shù)據(jù)的結(jié)果，該方法能夠較好地運用于實際的數(shù)值模擬中。

5 結(jié)論

通過本文分析可以得到以下結(jié)論：

（1）隨著模型邊界旋轉(zhuǎn)角度的增大，網(wǎng)格質(zhì)量逐漸下降；邊界旋轉(zhuǎn)角度增大到一定程度時，在翼尖周圍網(wǎng)格出現(xiàn)嚴(yán)重變形甚至發(fā)生網(wǎng)格邊交叉現(xiàn)象。

（2）通過引入網(wǎng)格質(zhì)量改進彈簧近似法，改進的彈簧近似法在網(wǎng)格整體質(zhì)量、局部畸形情況上都較傳統(tǒng)方法有所改善。

（3）在邊界發(fā)生旋轉(zhuǎn)時，傳統(tǒng)彈簧近似法生成網(wǎng)格變形主要集中在翼型前后緣處，容易發(fā)生局部網(wǎng)格變形集中；改進的彈簧近似法將變形分散到了整個計算區(qū)域，緩解了局部網(wǎng)格的變形壓力。

（4）在模型發(fā)生較大旋轉(zhuǎn)角度時，傳統(tǒng)彈簧近似法無法生成可用的網(wǎng)格，而通過對倔強系數(shù)的修正使最大旋轉(zhuǎn)角度有了較大的提高。

（5）運用改進的彈簧近似法對二維翼型進行數(shù)值模擬，并與實驗數(shù)據(jù)進行比較，得到了較滿意的結(jié)果，該方法能夠較好地運用于實際數(shù)值模擬中。

［1］ Rausch R D，Batina J.Spatial adaptation procedures on tetrahedral meshes for unsteady aerodynamic flow calculations［J］.AIAA Paper 93－0670，AIAA 31stAerospace Sciences Meeting.1993.

［2］ Baker T，Vassberg J.Tetrahedral mesh generation and optimization［J］.Proc.6thInternational Conference on Numerical Grid Generation，ISGG，1998：337－349.

［3］Baker T J，Cavallo P A.Dynamic adaptation for deforming tetrahedral meshes［J］.AIAA Paper，A99 － 3253，1999：19－29.

［4］Batina J T.Unsteady euler airfoil solutions using unstructured dynamic meshes［J］.AIAA Paper，AIAA 27thAerospace SciencesMeeting， Reno， Nevada， 1990，28（8）：1381－1388.

［5］Farhat C，Degand C，Koobus B，et al.An improved method of spring analogy for dynamic unstructured fluid meshes［J］.AIAA Paper，98 －2070.

［6］ Blom T J.Considerations on the spring analogy［J］.Journal of Aircraft，2000，32：647 －668.

［7］ Mitsuhiro Murayama，et al.Unstructured dynamic mesh for large movement and deformation［J］.AIAA 2002 －0122.

［8］ Zhang S J，Liu J，Chen Y S.A dynamic mesh method for unstructured flow solvers［J］.AIAA Computational Fluid Dynamics Conference，2003－3843.

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

［10］ Software M S C.MSC.patran user’s manual［M］.Los Angeles：The Macneal－ Schwendler Corporation，2005.

［11］ Singh K P，Newman J C，Baysal O.Dynamic unstructured method for flows past multiple objects in relative motion［J］.AIAA Journal，1995，33（4）：641 －649.

［12］ Abbott I H，Von Doenhoff A F.Theory of wing sections，including a summary of airfoil data［M］.New York：Dover Publications，1959.

［13］ Bohn D E，Moritz N.Algebraic method for efficient adaption of structured grids to fluctuating geometries［J］.Proceedings of the Institution of Mechanical Engineers，Part A：Journal of Power and Energy，2005，219（4）：303 －314.

［14］ Ivan M.Dynamic finite element meshes for 3D lagrangian CFD［J］.AIAA Paper，AIAA 16thComputational Fluid Dynamics Conference，Orlando，F(xiàn)lorida，June 2003.

［15］ Charbel Farhat，Ajaykumar Rajasekharan and Bruno Koobus.A dynamic variational multiscale method for large eddy simulations on unstructured meshes［J］.Comput.Methods Appl.Mech.Engrg，2006，195：1667 －1691.

［16］ Fernando de Goes，Siome Goldenstein，Luiz Velho.A simple and flexible framework to adapt dynamic meshes［J］.Computers ＆ Graphics，2008，32：141 －148.