彈簧近似法在二維非結構動網格生成技術中的應用

霍世慧， 王富生， 岳珠峰

（西北工業大學 力學與土木建筑學院，西安 710129）

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

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

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

1 彈簧近似法

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

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

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

由式（2）可得整個計算區域中m個節點在初始狀態下的矩陣表達式為：

當模型邊界發生變形時，通過改變模型邊界的位置矢量和增加適當的外邊界限制條件，如計算區域外邊界位置矢量保持不變、模型邊界位置矢量改為新位置下的矢量等，運用Gauss－seidel迭代法對式（3）進行迭代求解，其相應的分量迭代格式可表示式為：

當模型發生變形時，利用式（4）對式（3）進行數次迭代，達到所需要的精度時，便可得到計算區域中各網格節點變形后坐標。

2 彈簧倔強系數的選取

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

2.1 傳統方法

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

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

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

其中，θk為以i為頂點三角形面單元的內角。當邊界發生運動或變形時，對每個頂點均可得到關于拉伸彈簧的力平衡方程和扭轉彈簧的力矩平衡方程，聯立求解得出計算區域內的新網格點。

2.2 方法改進

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

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

圖1 典型畸形網格單元Fig.1 Classical abnormal meshes

在網格的變形過程中，通過兩個網格質量評定指標的修正，逐漸增大網格質量較差單元各棱邊的倔強系數，使網格質量較好單元承擔一定的變形。

3 方法應用

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

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

圖2 NACA23012初始氣動網格Fig.2 Original mesh of NACA23012

圖3 網格變形指標變化情況Fig.3 The change of mesh deformation value

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

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

圖4給出了傳統和改進彈簧近似法下，計算區域網格整體變形指標fA、fAR隨翼型旋轉角度的變化情況。

圖4 變形指標隨旋轉角度變化曲線Fig.4 Deformation value with the increase of rotation angle

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

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

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

圖5 網格質量隨旋轉角度變化曲線Fig.5 Mesh quality with the increase of rotation angle

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

4 算例

4.1 算例一

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

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

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

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

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

4.2 算例二

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

由圖13可以發現，改進彈簧近似法生成的動網格在二維諧和振動翼型非定常氣動力的求解中能夠得到較傳統方法更接近于實驗數據的結果，該方法能夠較好地運用于實際的數值模擬中。

5 結論

通過本文分析可以得到以下結論：

（1）隨著模型邊界旋轉角度的增大，網格質量逐漸下降；邊界旋轉角度增大到一定程度時，在翼尖周圍網格出現嚴重變形甚至發生網格邊交叉現象。

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

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

（4）在模型發生較大旋轉角度時，傳統彈簧近似法無法生成可用的網格，而通過對倔強系數的修正使最大旋轉角度有了較大的提高。

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

［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］史愛明.非結構動網格下的非定常氣動力計算和嗡鳴研究［D］.西安：西北工業大學，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，Florida，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.