基于TFI的局部網格變形方法研究

楊維龍，楊永

(西北工業大學 航空學院，西安　710072)

基于TFI的局部網格變形方法研究

楊維龍，楊永

(西北工業大學 航空學院，西安710072)

摘要：整體網格變形技術對其他網格技術的使用產生影響，難以相互配合來解決復雜的流動問題。發展一套基于TFI(超限插值法)技術的局部網格變形方法，該方法不同于已有的整體網格變形思路，將網格變形控制在有限的區域內，在保證網格質量的同時，不影響其他區域的網格，因此不會對其他動網格方法的使用造成影響。以帶后緣襟翼的NACA0012翼型為例，對翼型俯仰振蕩耦合襟翼偏轉運動的動態流場進行數值模擬，結果表明：計算結果與實驗值吻合良好，證明了本文發展的局部網格變形方法的可行性。

關鍵詞：襟翼偏轉；網格變形；俯仰振蕩；非定常流動；超限插值法

0引言

現代飛機設計對計算流體力學(CFD)的要求越來越高，在數值模擬飛機部件具有相對運動的復雜繞流時，需要合適的動網格技術使數值模擬過程更加高效準確，網格變形技術是經常使用的一種動網格方法。

網格變形技術是在數值模擬的過程中，隨著時間的推進，根據物體的運動規律，使用數值方法更新網格。常用的網格變形方法有超限插值法(TFI)[1]和彈簧類推法[2]。TFI方法是一種代數方法，因此網格變形效率很高。雖然該方法能夠在一定程度上保證初始網格的質量，但是只適用于結構網格，并且只適用于網格小幅度變形的情況。彈簧類推法能夠解決大幅度的變形問題，但是由于需要迭代求解靜力平衡方程，因此需要消耗較多的計算時間。

綜上所述，相對于彈簧類推法，TFI方法的變形效率高，適應能力強，易于編程，因此該方法得到了廣泛的應用和發展。L.Dubuc等[3]基于TFI方法開發了一套網格變形算法，并通過數值模擬Williams翼型襟翼的強迫振蕩驗證了該方法的可行性。袁先旭[4]在TFI的基礎上，提出了加權TFI動網格生成方法，生成了OREX飛船返回艙外形俯仰振蕩的結構動網格。張兵等[5]通過引入旋轉矢量的思想，提出了一種帶旋轉修正的彈簧-TFI混合動網格技術，解決了大幅度變形情況下物面附近網格的正交性問題。

對于飛機具有多個部件相對運動的數值模擬，僅使用網格變形方法有時是無法完成的，需要多種動網格方法配合使用。但是，對于多塊網格，由于物面運動帶來的擾動是在整個流場區域網格中傳播的，通常的網格變形程序是將整個求解區域的網格作為變形對象。然而，在網格變形過程中，物面網格的運動對離該物面較遠的網格帶來的影響較小，甚至可以忽略。若對這部分網格進行網格變形處理，會造成計算機資源的浪費。并且由于整個求解區域的網格都發生變形，也制約了該方法和其他動網格方法的配合使用。例如，當模擬某螺旋槳飛機的螺旋槳旋轉和副翼偏轉時，可以使用網格變形技術模擬副翼的偏轉[6]，使用動態面搭接技術模擬螺旋槳的旋轉[7]。但是，動態面搭接技術要求搭接面形狀保持不變，整體網格變形的策略就會造成搭接面的變形，造成兩種方法無法配合使用。

整體網格變形策略會給其他動網格技術的使用帶來影響，為了解決該問題，本文首先發展一套基于TFI技術的局部網格變形方法，將網格變形的區域限定在指定的范圍內，從而不會影響變形域外的網格；其次對NACA0012翼型俯仰振蕩耦合襟翼偏轉運動的動態流場進行數值模擬。

1網格變形方法

1.1局部網格變形介紹

局部網格變形，指的是只有部分網格發生變形，其余的網格并不發生變形；也可以是整體網格剛性運動后部分網格塊發生變形。基于TFI技術的局部網格變形方法的基本思路是：首先確定網格的變形域及邊界，并給定邊界的運動形式，該邊界包括運動的物面網格和變形域的邊界面網格；然后根據需求劃分子網格，建立主控點和受控點的對應關系[8]，再根據給定的運動形式確定邊界點的位移，使用基于弧長的TFI公式計算子網格內部網格點的位移；最后更新網格，得到變形后的網格。

子網格是將網格塊每隔若干個網格點取出一個網格點所得到的，子網格劃分的示意圖如圖1所示，將Block1分割成了9個子網格塊，A、B、C、D四點為子網格塊5的角點。

圖1　子網格劃分示意圖

1.2局部網格變形具體過程

令運動的物面和變形域邊界上的網格點為主控點，子網格的角點為受控點，每個受控點可以根據實際情況具有一個或者多個主控點，受控點的位移受到主控點的位移以及與其最近的主控點的距離影響。設xci為子網格塊上第i個受控點坐標，xsj為物面上第j個主控點坐標，則兩點的距離為

|rij|=‖xci-xsj‖2

(1)

受控點與主控點的距離越近，受控點的變形量受到的影響就越大。通過建立衰減函數確定兩者之間的關系，該衰減函數為

(2)

(3)

(4)

(5)

在確定子網格角點位移之后，先對內部網格進行坐標歸一化以便使用基于弧長的TFI方法計算內部網格點的位移，弧長歸一化的公式為

(6)

(7)

二維TFI的公式為

Δxi,j=U+V-UV

(8)

U=(1-αi,j)Δx1,j+αi,jΔximax,j

(9)

V=(1-βi,j)Δxi,1+βi,jΔxi,jmax

(10)

UV=αi,j(1-βi,j)Δximax,1+βi,j(1-αi,j)Δx1,jmax+

αi,jβi,jΔximax,jmax+(1-αi,j)(1-βi,j)Δx1,1

(11)

1.3流體力學控制方程

CFD在流體域求解守恒型積分的形式適用于動網格的N-S方程[9]

(12)

式中：t為時間；q為守恒變量矢量，表示單位體積的密度、動量和總能量；Ω(t)和?Ω(t)為控制體積及其邊界；n為控制邊界?Ω(t)的外法向；Fc和Fv分別為運輸和粘性通量；H為源項。

流動控制方程空間離散采用有限體積方法，無粘通量項采用二階精度的Roe通量差分分裂迎風格式離散，粘性項采用二階中心格式進行離散。采用全湍流假設，Sparlart-Allmas湍流模型進行湍流計算。全隱式的方程求解，非定常時間推進采用雙時間步法。

2算例驗證

2.1NACA0012翼型俯仰振蕩

本文使用的動網格方案包括剛性動網格技術和局部網格變形技術，為了驗證兩種方案的可行性，分別使用上述兩種網格技術模擬NACA0012翼型在跨音速流場中繞1/4弦點做俯仰振蕩的非定常流場[10]，來流馬赫數Ma=0.755。翼型迎角的變化規律為

(13)

本算例使用整體網格，網格為C型網格(如圖2所示)，網格量為16 016，物面附近的網格如圖3所示。翼型偏轉2.51°后的網格狀態如圖4所示。

圖2　整體網格

圖3　物面網格

圖4　變形后的網格

從圖3～圖4可以看出：此時的物面網格仍保持很好的正交性。在小角度的變形幅度下，基于TFI的網格變形技術能夠很好地滿足變形的需求，并保證了網格質量。由于TFI是一種代數的方法，數值求解時計算網格變形的速度非常快。

計算結果與實驗值的對比如圖5所示。

(a) CL

(b) CM

從圖5可以看出：使用兩種方法得到的升力系數和俯仰力矩系數均與實驗值吻合良好，表明本文使用的這兩種動網格方法均是可行的。

2.2翼型與襟翼共同振蕩的數值模擬

2.2.1模型與計算網格

模型的基本翼型是NACA0012，弦長(包含襟翼)c=180mm，其中襟翼的長度是40.69mm，為26.15 c，翼型的旋轉軸位于0.35 c處，襟翼的旋轉軸位于0.8 c處，實驗模型的展長為600mm。計算狀態為Ma=0.4，Re=1.63×106，實驗模型中，襟翼和翼型之間縫隙的最大距離是0.5mm。該模型的示意圖如圖6所示[11]。

圖6　帶有后緣襟翼的翼型

本文使用的網格整體拓撲結構是C形(如圖7所示)，為了保證襟翼在振蕩過程中網格變形后的質量，在翼型周圍生成一個O型網格，網格數量為32 192，物面網格點數為263個。

實驗中，翼型和襟翼同時在做不同頻率的俯仰振蕩，因此，本文使用整體網格剛性旋轉實現翼型的振蕩，而襟翼與翼型之間的偏轉使用網格變形來實現。數值模擬的過程中，襟翼的偏轉使其附近的流場空間網格發生了變形，然而，在距離襟翼比較遠的流場區域內，襟翼偏轉引起的網格變形量十分微小，可以忽略。因此，本文采用局部變形的策略(如圖8所示)，只需圖中顯示網格的區域發生變形，足以滿足襟翼在偏轉過程中的網格變形需求。在數值模擬開始前，在參數文件中給定相應的網格變形邊界及運動形式，以實現數值模擬過程中的網格局部變形。

圖8　局部網格變形策略

2.2.2數值模擬結果及分析

通過求解N-S方程，對帶襟翼偏轉運動的翼型俯仰振蕩進行非定常求解。其中，翼型的俯仰振蕩控制方程為

α(t)=α0+Δαsin(ωαt+φα)

(14)

式中：α(t)為翼型整體的實時迎角；α0為翼型整體振動的初始迎角；Δα為翼型整體振蕩的幅值迎角；ωα為振動圓頻率；φα為相位差。

襟翼偏轉的運動控制方程為

δ(t)=δ0+Δδsin(ωδt+φδ)

(15)

式中：δ(t)為襟翼的實時偏轉角；δ0為襟翼的初始偏轉角；Δδ為襟翼的偏轉角幅值；ωδ為振蕩頻率；φδ為相位差。

該系統的真實運動形式為方程(14)和方程(15)的疊加。

對帶襟翼翼型的兩種運動狀態進行計算，并與實驗值進行對比，其中兩種狀態的計算使用同一套網格。

(1) 狀態一

該狀態下，帶襟翼偏轉運動的翼型俯仰振蕩的運動方程為

(16)

該運動方程如圖9所示。通過非定常求解，升力系數遲滯環曲線和俯仰力矩系數遲滯環曲線如圖10所示，并與文獻[11]中的實驗值進行對比。從圖10可以看出：計算結果與實驗值吻合良好；初始迎角為0°，幅值迎角為6°，迎角范圍為(-6°，6°)，在此基礎上，襟翼偏角為(-5°，6°)，因此本文所使用的方法對于小迎角范圍內的襟翼動態偏轉問題可以進行較好地數值模擬，具有一定的計算精度；由于迎角振蕩頻率為襟翼偏轉頻率的一半，因此遲滯環曲線具有明顯的交叉。

圖9　翼型迎角和襟翼偏角運動規律(狀態一)

(a) CL

(b) CM

一個周期內帶襟翼翼型的動態變化流場壓力系數云圖及襟翼處的局面網格細節圖如圖11所示。襟翼偏轉后，仍然保有較好的局部網格質量，因此本文所使用的網格變形方法對襟翼偏轉問題具有較好的適用性。

(a1) 壓力系數云圖 (a2) 局部網格細節

(a)t=0.005T

(b1) 壓力系數云圖 (b2) 局部網格細節

(b)t=0.125T

(c1) 壓力系數云圖 (c2) 局部網格細節

(c)t=0.250T

(d1) 壓力系數云圖 (d2) 局部網格細節

(d)t=0.375T

(e1) 壓力系數云圖 (e2) 局部網格細節

(e)t=0.500T

(f1) 壓力系數云圖 (f2) 局部網格細節

(f)t=0.625T

(g1) 壓力系數云圖 (g2) 局部網格細節

(g)t=0.750T

(h1) 壓力系數云圖 (h2) 局部網格細節

(h)t=0.875T

圖11各時間節點的壓力系數云圖及局部網格細節

Fig.11Cpcontours and local mesh details of each time node

從圖11可以看出：在t=0.005T時，翼型為正迎角，具有抬頭角速度，襟翼偏轉角為負值(根據飛行力學定義，襟翼下偏為正)，且具有向下偏轉的角速度，此時翼型上表面低壓區較強，但由于襟翼上偏，因此翼型整體具有負升力；在t=0.125T時，翼型為正迎角，具有抬頭角速度，襟翼偏轉角為正值，此時翼型上表面低壓區較強，翼型整體具有正升力；在t=0.250T時，翼型為正迎角，且達到峰值，襟翼偏轉角為正值，且具有向上偏轉的角速度，翼型整體具有正升力；在t=0.375T時，翼型為正迎角，具有低頭角速度，襟翼偏轉角為負值；在t=0.500T時，翼型為負迎角，具有低頭角速度，襟翼偏轉角為負值，翼型整體升力為負值；在t=0.625T時，翼型為負迎角，翼型上表面高壓區較強，盡管襟翼偏轉角為正值，但翼型整體升力為負值；在t=0.750T時，翼型為負迎角，達到迎角峰值，襟翼偏轉角為正值，達到襟翼偏角峰值，翼型整體升力為正值；在t=0.875T時，翼型為負迎角，襟翼偏轉角為負值，翼型整體升力為負值。

(2) 狀態二

該狀態下，帶襟翼偏轉運動的翼型俯仰振動的運動方程為

(17)

該運動方程如圖12所示。

圖12　翼型迎角和襟翼偏角運動規律圖(狀態二)

通過非定常求解，升力系數遲滯環曲線和俯仰力矩系數遲滯環曲線如圖13所示，并與文獻[11]中的實驗值進行對比。

(a) CL

(b) CM

從圖13可以看出：計算結果與實驗值吻合良好；初始迎角為4.25°，幅值迎角為5.75°，迎角范圍為(-1.5°，10°)，在此基礎上，襟翼偏角為(-6°，5°),因此本文所使用的方法對于中等迎角范圍內的襟翼動態偏轉問題可以進行較好地數值模擬，具有一定的計算精度。

3結論

(1) 本文發展了一套基于TFI技術的局部網格變形方法，該方法能夠在保證計算精度的情況下提高計算效率，并且可以較好地配合其他動網格技術以模擬有多部件相對運動的復雜流場。

(2) 本文結合局部網格變形方法和剛性動網格方法數值模擬了NACA0012翼型與襟翼以不同頻率進行俯仰振蕩的復雜流場，變形后的網格依然具有很好的網格質量，將計算結果與實驗值進行對比，驗證了本文發展方法的可行性。

(3) 本文僅介紹了局部網格變形方法在二維問題上的應用，對于帶副翼動態偏轉的三維問題，也可以使用該方法結合動態面搭接技術進行數值模擬，從而有效地解決由于副翼偏轉引起的“剪刀差”問題。

參考文獻

[1] 楊國偉, 錢衛. 飛行器跨聲速氣動彈性數值分析[J]. 力學學報, 2005, 37(6)： 769-776.

Yang Guowei, Qian Wei. Numerical analyses of transonic flutter on an aircraft[J]. Chinese Journal of Theoretical and Applied Mechanics, 2005, 37(6): 769-776. (in Chinese)

[2] Eriksson L E. Generation of boundary-conforming grids around wing-body configurations using transfinite interpolation[J]. AIAA Journal, 1982(20): 1313-1320.

[3] Dubuc L, Cantariti F, Woodgate M, et al. A grid deformation technique for unsteady flow computations[J]. International Journal for Numerical Methods in Fluids, 2000, 32(3): 285-311.

[4] 袁先旭. 非定常流動數值模擬及飛行器動態特性分析研究[D]. 綿陽:中國空氣動力研究與發展中心, 2004.

Yuan Xianxu. Numerical simulation of unsteady flow and dynamic characteristics research of aircraft[D]. Mianyang: Computational Aerodynamics Institute of China Aerodynamics Research and Development Center, 2004.(in Chinese)

[5] 張兵, 韓景龍. 帶旋轉修正的彈簧TFI-混合動網格方法[J]. 航空學報, 2011, 32(10): 1815-1823.

Zhang Bing, Han Jinglong. Spring-TFI hybrid dynamic mesh method with rotation correction[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(10): 1815-1823.(in Chinese)

[6] 李喜樂， 楊永. 帶副翼偏轉的三角翼自由滾轉運動數值模擬[J]. 航空學報, 2012， 33(3)： 453-462.

Li Xile, Yang Yong. Numerical simulation of the free rolling motion of a delta wing configuration with aileron deflection[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33(3): 453-462.(in Chinese)

[7] 夏貞鋒, 楊永. 螺旋槳滑流與機翼氣動干擾的非定常數值模擬[J]. 航空學報, 2011, 32(7): 1195-1201.

Xia Zhenfeng, Yang Yong. Unsteady numerical simulation of interaction effects of propeller and wing[J]. Acta Aeronauticaet Astronautica Sinica, 2011, 32(7): 1195-1201.(in Chinese)

[8] Peter M Hartwich, Shreekant Agrawal. Method for perturbing multiblock patched grids in aeroelastic and design optimization applications[R]. AIAA-97-2038, 1997.

[9] Anderson John D, Jr. Fundamentals of aerodynamics[M]. New York: McGraw-Hill Book Company, 2007.

[10] Landon R H. NACA0012 oscillation and transient pitching[R]. Compendium of Unsteady Aerodynamic Measurements, Data Set 3, AGARD-R-702, 1982.

[11] Andrzej Krzysiak, Janusz Narkiewicz. Aerodynamic loads on airfoil with trailing-edge flap pitching with different frequencies[J]. Journal of Aircraft, 2006, 43(2): 407-418.

Investigation on Local Grid Deformation Method Based on TFI

Yang Weilong, Yang Yong

(School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China)

Abstract:Global mesh deformation technique has influence on other dynamic grid methods, so it is difficult to be combined with other dynamic grid methods to solve complex flow problems. A local grid deformation method based on TFI(Transfinite Interpolation) is developed. This method, which is different from other global grid deformation methods, is that it can limit the deformation in specific zones without any influence on grids of other zones, so it will not hamper the use of other dynamic grid methods. For an airfoil NACA0012 with a flap, numerical simulations of pitching oscillations are conducted, and computational results agrees well with experimental results. Then it is concluded that the method developed in this paper is viable for flow solutions of complex geometries.

Key words:flap deflection; grid deformation; pitching oscillation; unsteady flow; transfinite interpolation

收稿日期：2016-04-01；修回日期：2016-04-25

通信作者：楊維龍,ywlvs@163.com

文章編號:1674-8190(2016)02-216-09

中圖分類號：V211.3

文獻標識碼：A

DOI：10.16615/j.cnki.1674-8190.2016.02.012

作者簡介：

楊維龍(1988-)，男，碩士研究生。主要研究方向：動網格技術。

楊永(1962-)，男，博士，教授，博導。主要研究方向：空氣動力學、計算流體力學和設計空氣動力學。

(編輯：趙毓梅)