李楠+倪原+李聚峰+牛佳慧+田華

摘 要： 課題采用湍流模型中的單方程模型，并借助于商業(yè)CFD軟件Fluent，對某型號飛行器在給定攻角和馬赫數下進行氣動力仿真。首先，采用湍流模型中的單方程模型對飛行器的繞流流場進行數學建模；其次，用Gambit軟件對飛行器的外形進行幾何建模，并進行網格劃分和邊界條件設定；最后，在Fluent中進行相關參數設置和氣動力數值計算及仿真，得出了該型號飛行器升力系數、阻力系數和俯仰力矩系以及飛行器表面的壓力分布和來流速度分布，并對結果進行分析。該研究表明單方程模型能夠很好地解決具有壁面限制的飛行器氣動力數值計算問題。

關鍵詞： 飛行器； 氣動力； 湍流模型； Spalart?Allmaras

中圖分類號： TN911?34 文獻標識碼： A 文章編號： 1004?373X（2014）16?0068?03

Fluent?based calculation method of aircraft aerodynamic parameters

LI Nan1， NI Yuan2， LI Ju?feng2， NIU Jia?hui2， TIAN Hua2

（1. Xian Institute of Electromechanical Information Technology， Xian 710065， China；

2. College of Electronics and Information Engineering， Xian Technological University， Xian 710032， China）

Abstract： The aerodynamics simulation of an aircraft was performed in a given attack angle and Mach number by means of the single?equation （Spalart?Allmaras） model in turbulence model and CFD software Fluent. First， Spalart?Allmaras model is used to achieve the mathematical modeling for the flow field of the aircraft， the Fluent software package， Gambit is adopted to conduct the geometric modeling of the aircrafts outline， and then meshes are generated and the boundary conditions are set. Finally， the relevant parameters are set in Fluent， and lift coefficient， drag coefficient， pitching moment coefficient， and pressure and flow velocity distribution on the aircraft surface are obtained by the aerodynamic numerical calculation and simulation. The simulation results also need to be analyzed. The research result demonstrates the Spalart?Allmaras model can achieve the aerodynamic numerical calculation of the aircraft with wall?surface restriction.

Keywords： aircraft； aerodynamics； turbulence model； Spalart?Allmaras

飛行器在飛行過程中受到周圍氣流的作用，其表面的壓強分布不均勻，因此引起的壓力差和氣流對飛行器表面產生的粘性摩擦力合共同作用，形成了飛行器上的空氣動力[1]，而阻力、升力和俯仰力矩等是研究計算空氣動力的重要參數。在以往的研究中，氣動力參數的獲得多數依靠試驗或根據飛行器的外形進行大量的計算，例如風洞試驗和實彈試驗等。但這樣的研制周期長、試驗耗資大、成本高，并且傳統試驗方法的精準度有限。

隨著計算機和計算流體力學（Computational Fluid Dynamics，CFD）軟件的發(fā)展，仿真可以提高系統性能和研制質量，縮短研制周期、減少實驗成本。文中采用更先進的Fluent軟件對空氣動力進行計算，使所需的氣動力參數的獲得變得更加簡便，精度更高。本文采用Fluent軟件對研究的飛行器進行建模和仿真，并計算得出所需的阻力，升力和俯仰力矩等參數。

1 數學模型

在CFD軟件出現之前，計算空氣動力常采用的風洞實驗，雖然它能夠比較準確地控制實驗條件和實驗項目內容的多樣性等優(yōu)點，但是邊界效應、支架干擾和相似準則不能滿足能缺點卻無法避免[2]。隨著CFD的出現，尤其是Fluent軟件的產生和發(fā)展，使空氣動力的計算變得越發(fā)便捷和準確，且克服了風洞實驗的不足之處。Fluent用來模擬從不可壓縮到高度可壓縮范圍內的復雜流動。由于采用了多重網格加速收斂技術和多種求解方法，因此Fluent能達到最佳的求解精度和收斂速度 [3]。

Fluent包含豐富的物理模型，例如計算流體流動和熱傳導模型（包括層流，定常和非定常流動，自然對流、湍流，不可壓縮和可壓縮流動，紊流，周期流，旋轉流，時間相關流等）；相變模型，輻射模型，離散相變模型，多相流模型，化學組分輸運和反應流模型等。對每一種物理問題的流動特點都有適合它的數值解法，可對顯式或隱式差分格式進行選擇，可以在計算速度、穩(wěn)定性、精度等方面達到最佳[4]。

Fluent提供的湍流模型包括：雙方程模型（重整化群κ?ε模型、標準κ?ε模型、可實現（Realizable）κ?ε模型），單方程（Spalart?Allmaras）模型、雷諾應力模型和大渦模擬[5]。其中單方程模型相對簡單，對網格的質量較為不敏感，只求解一個有關渦粘性的輸運方程，計算量相對較?。徊⑶以撏牧髂Ｐ捅容^適合于具有壁面限制的流動問題，對有逆壓梯度的邊界層問題能夠給出很好的計算結果[6?7]。Spalart?Allmaras模型的求解變量是[ν]，表征出了近壁（粘性影響）區(qū)域以外的湍流運動粘性系數。[ν]的輸運方程為：

[ρdνdt=Gν+1σν??xj（μ+ρν）?ν?xj+Cb2ρ?ν?xj-Yν]

式中：[Gν]是湍流粘性產生項；[Yν]是由于壁面阻擋與粘性阻尼引起的湍流粘性的減少；[σν]和[Cb2]是常數；ν是分子運動粘性系數。對于本文所研究的飛行器模型，屬于雷諾數很大，黏度很小，有邊界層的湍流，并且對網格劃分要求不高，因此選用單方程模型進行模擬運算。

2 計算方法與過程

用Fluent軟件包的前置處理器GAMBIT軟件，按照由點到線，由線到面，由面到體的原則對飛行器進行建模，并根據尺寸比例在外圍設置一個繞流流場，繞流流場設置過程為：外圍是一個半徑為7 200 mm，高為3 600 mm的大圓柱形，用來作為飛行器的繞流流場。大圓柱內部設了一個小圓柱體，其半徑為1 800 mm，高為3 600 mm（飛行器位于計算域的正中），并進行較密的網格劃分，兩個圓柱之間的區(qū)域可將網格劃分得稀疏些，這樣既保證了計算的準確性，又減少了計算量。飛行器的舵翼部位采用密集的網格劃分，以反映氣流的劇烈變化，頭部和彈體采用兩頭密、中間稀的網格劃分；由于飛行器是軸對稱的，為了提高計算效率，減少計算時間，取整個模型的一半進行計算。整個區(qū)域的網格有2 018 432個。

對飛行器的相關區(qū)域進行網格劃分，如圖1所示。

把生成的網格文件導入FLUENT求解器里進行求解，具體的求解過程如下：

（1） 粘性模型采用單方程模型，使用運輸方程求解；

（2） 設置流體材料屬性：密度為ideal?gas，材料設為air，在“Viscosity”一項中選Sutherland，采用三系數方法[8]，對于高速可壓縮流動氣體采用描述氣體粘度的Suther?land 定律較比較合適；

（3） 壁面條件：無滑移，壁面粗糙度選為0.5，所有其他標量選擇不可滲透壁面條件；

（4） 數值計算過程中的差分格式選擇：動量、湍流動能、湍流耗散率均選用second order upwind scheme，即二階迎風格式；壓力插值選用默認的standard 方法[9]；

（5） 邊界條件設置：設置大氣壓為101 325 Pa，馬赫數為0.8，攻角為4°，粘性比為10；

（6） 松弛因子的設置：求解時設置壓力項松弛因子為0.8，密度、質量力項為1，動量項為0.5，湍動能項為0.6，耗散率項為0.6，湍流粘性項為0.6；

（7） 設置迭代次數6 500次，求解。

圖1 計算域網格劃分

3 計算實例

在給定的攻角和馬赫數下，阻力系數、升力系數和俯仰力矩系數隨著迭代次數的增加而也跟著振蕩，一直到迭代6 000次以后，曲線的變化逐漸趨于穩(wěn)定，待曲線穩(wěn)定以后，就可以在圖中讀出阻力系數、升力系數、俯仰力矩系數等相關參數的值。圖2為阻力系數隨迭代過程變化曲線，圖3為升力系數隨迭代過程變化曲線，圖4為俯仰力矩系數隨迭代過程變化曲線。從圖2～圖4可以看出，其氣動特性是收斂的，表明所設計的系統是穩(wěn)定的。仿真出的氣動參數和在該攻角和馬赫數下的試驗參數表對比之后有誤差，誤差為0.05%，在誤差允許范圍內，因此可以驗證出仿真數據的準確性。這些氣動參數將為下一步實驗提供數據。

圖2 阻力系數隨迭代過程變化曲線

圖3 升力系數隨迭代過程變化曲線

圖4 俯仰力矩系數隨迭代過程變化曲線

如果采用湍流模型中的其他模型，比如雙方程模型，計算過程相對復雜，對計算網格的要求更加精確，即壁面條件要求更高，一旦網格劃分不符合要求，則計算時很難得到收斂曲線，并且在網格劃分準確的情況下計算次數要達到10 000萬次，而在實際項目中，不需要對飛行器的網格進行精確的劃分，而單方程模型簡單、計算次數少，對網格的要求不敏感，有效計算出了實驗所需的相關參數，大大地節(jié)省了實驗時間。

4 結 論

在給定攻角和馬赫數的前提下，采用湍流模型中的單方程模型，借助CFD軟件Fluent對某飛行器進行氣動力計算，可以得出飛行器的升力系數、阻力系數、俯仰力矩系數。這些參數將為下一步飛行器控制系統設計提供依據。

參考文獻

[1] 趙洪章，岳春國，李進賢.基于Fluent的飛行器氣動特性計算[J].彈箭與制導學報，2007，27（2）：203?205.

[2] 于勇.FLUENT入門與進階教程[M].北京：北京理工大學出版社，2008.

[3] 劉明侯.計算流體和傳熱傳質[D].合肥：中國科技大學，2005.

[4] 龐英良，宋衛(wèi)東，楊曉霖.基于Fluent 的末制導炮彈初始段氣動仿真[J].兵工自動化，2009，28（2）：13?15.

[5] 錢翼稷.空氣動力學[M].北京：北京航空航天大學出版社，2009.

[6] 尹志林.某隱身巡航導彈氣動及雷達目標特性分析[D].南京：南京航空航天大學，2009.

[7] 張偉，張安堂，肖宇.基于坐標旋轉數字計算方法的三維坐標變換[J].探測與控制學報，2011，33（2）：73?76.

[8] 朱銳，董二寶，張杰，等.頭部可偏轉飛行器氣動仿真與外形優(yōu)化[J].機械與電子，2008（8）：6?8.

[9] LANDERS M G， HALL L H. Deflectable nose and canard controls for a fin?stabilized projectile at supersonic and hypersonic speeds， AIAA 2003?3805 [R]. Orlando， Florida， USA： AIAA， 2003.

Fluent提供的湍流模型包括：雙方程模型（重整化群κ?ε模型、標準κ?ε模型、可實現（Realizable）κ?ε模型），單方程（Spalart?Allmaras）模型、雷諾應力模型和大渦模擬[5]。其中單方程模型相對簡單，對網格的質量較為不敏感，只求解一個有關渦粘性的輸運方程，計算量相對較??；并且該湍流模型比較適合于具有壁面限制的流動問題，對有逆壓梯度的邊界層問題能夠給出很好的計算結果[6?7]。Spalart?Allmaras模型的求解變量是[ν]，表征出了近壁（粘性影響）區(qū)域以外的湍流運動粘性系數。[ν]的輸運方程為：

[ρdνdt=Gν+1σν??xj（μ+ρν）?ν?xj+Cb2ρ?ν?xj-Yν]

式中：[Gν]是湍流粘性產生項；[Yν]是由于壁面阻擋與粘性阻尼引起的湍流粘性的減少；[σν]和[Cb2]是常數；ν是分子運動粘性系數。對于本文所研究的飛行器模型，屬于雷諾數很大，黏度很小，有邊界層的湍流，并且對網格劃分要求不高，因此選用單方程模型進行模擬運算。

2 計算方法與過程

用Fluent軟件包的前置處理器GAMBIT軟件，按照由點到線，由線到面，由面到體的原則對飛行器進行建模，并根據尺寸比例在外圍設置一個繞流流場，繞流流場設置過程為：外圍是一個半徑為7 200 mm，高為3 600 mm的大圓柱形，用來作為飛行器的繞流流場。大圓柱內部設了一個小圓柱體，其半徑為1 800 mm，高為3 600 mm（飛行器位于計算域的正中），并進行較密的網格劃分，兩個圓柱之間的區(qū)域可將網格劃分得稀疏些，這樣既保證了計算的準確性，又減少了計算量。飛行器的舵翼部位采用密集的網格劃分，以反映氣流的劇烈變化，頭部和彈體采用兩頭密、中間稀的網格劃分；由于飛行器是軸對稱的，為了提高計算效率，減少計算時間，取整個模型的一半進行計算。整個區(qū)域的網格有2 018 432個。

對飛行器的相關區(qū)域進行網格劃分，如圖1所示。

把生成的網格文件導入FLUENT求解器里進行求解，具體的求解過程如下：

（1） 粘性模型采用單方程模型，使用運輸方程求解；

（2） 設置流體材料屬性：密度為ideal?gas，材料設為air，在“Viscosity”一項中選Sutherland，采用三系數方法[8]，對于高速可壓縮流動氣體采用描述氣體粘度的Suther?land 定律較比較合適；

（3） 壁面條件：無滑移，壁面粗糙度選為0.5，所有其他標量選擇不可滲透壁面條件；

（4） 數值計算過程中的差分格式選擇：動量、湍流動能、湍流耗散率均選用second order upwind scheme，即二階迎風格式；壓力插值選用默認的standard 方法[9]；

（5） 邊界條件設置：設置大氣壓為101 325 Pa，馬赫數為0.8，攻角為4°，粘性比為10；

（6） 松弛因子的設置：求解時設置壓力項松弛因子為0.8，密度、質量力項為1，動量項為0.5，湍動能項為0.6，耗散率項為0.6，湍流粘性項為0.6；

（7） 設置迭代次數6 500次，求解。

圖1 計算域網格劃分

3 計算實例

在給定的攻角和馬赫數下，阻力系數、升力系數和俯仰力矩系數隨著迭代次數的增加而也跟著振蕩，一直到迭代6 000次以后，曲線的變化逐漸趨于穩(wěn)定，待曲線穩(wěn)定以后，就可以在圖中讀出阻力系數、升力系數、俯仰力矩系數等相關參數的值。圖2為阻力系數隨迭代過程變化曲線，圖3為升力系數隨迭代過程變化曲線，圖4為俯仰力矩系數隨迭代過程變化曲線。從圖2～圖4可以看出，其氣動特性是收斂的，表明所設計的系統是穩(wěn)定的。仿真出的氣動參數和在該攻角和馬赫數下的試驗參數表對比之后有誤差，誤差為0.05%，在誤差允許范圍內，因此可以驗證出仿真數據的準確性。這些氣動參數將為下一步實驗提供數據。

圖2 阻力系數隨迭代過程變化曲線

圖3 升力系數隨迭代過程變化曲線

圖4 俯仰力矩系數隨迭代過程變化曲線

如果采用湍流模型中的其他模型，比如雙方程模型，計算過程相對復雜，對計算網格的要求更加精確，即壁面條件要求更高，一旦網格劃分不符合要求，則計算時很難得到收斂曲線，并且在網格劃分準確的情況下計算次數要達到10 000萬次，而在實際項目中，不需要對飛行器的網格進行精確的劃分，而單方程模型簡單、計算次數少，對網格的要求不敏感，有效計算出了實驗所需的相關參數，大大地節(jié)省了實驗時間。

4 結 論

在給定攻角和馬赫數的前提下，采用湍流模型中的單方程模型，借助CFD軟件Fluent對某飛行器進行氣動力計算，可以得出飛行器的升力系數、阻力系數、俯仰力矩系數。這些參數將為下一步飛行器控制系統設計提供依據。

參考文獻

[1] 趙洪章，岳春國，李進賢.基于Fluent的飛行器氣動特性計算[J].彈箭與制導學報，2007，27（2）：203?205.

[2] 于勇.FLUENT入門與進階教程[M].北京：北京理工大學出版社，2008.

[3] 劉明侯.計算流體和傳熱傳質[D].合肥：中國科技大學，2005.

[4] 龐英良，宋衛(wèi)東，楊曉霖.基于Fluent 的末制導炮彈初始段氣動仿真[J].兵工自動化，2009，28（2）：13?15.

[5] 錢翼稷.空氣動力學[M].北京：北京航空航天大學出版社，2009.

[6] 尹志林.某隱身巡航導彈氣動及雷達目標特性分析[D].南京：南京航空航天大學，2009.

[7] 張偉，張安堂，肖宇.基于坐標旋轉數字計算方法的三維坐標變換[J].探測與控制學報，2011，33（2）：73?76.

[8] 朱銳，董二寶，張杰，等.頭部可偏轉飛行器氣動仿真與外形優(yōu)化[J].機械與電子，2008（8）：6?8.

[9] LANDERS M G， HALL L H. Deflectable nose and canard controls for a fin?stabilized projectile at supersonic and hypersonic speeds， AIAA 2003?3805 [R]. Orlando， Florida， USA： AIAA， 2003.

Fluent提供的湍流模型包括：雙方程模型（重整化群κ?ε模型、標準κ?ε模型、可實現（Realizable）κ?ε模型），單方程（Spalart?Allmaras）模型、雷諾應力模型和大渦模擬[5]。其中單方程模型相對簡單，對網格的質量較為不敏感，只求解一個有關渦粘性的輸運方程，計算量相對較?。徊⑶以撏牧髂Ｐ捅容^適合于具有壁面限制的流動問題，對有逆壓梯度的邊界層問題能夠給出很好的計算結果[6?7]。Spalart?Allmaras模型的求解變量是[ν]，表征出了近壁（粘性影響）區(qū)域以外的湍流運動粘性系數。[ν]的輸運方程為：

[ρdνdt=Gν+1σν??xj（μ+ρν）?ν?xj+Cb2ρ?ν?xj-Yν]

式中：[Gν]是湍流粘性產生項；[Yν]是由于壁面阻擋與粘性阻尼引起的湍流粘性的減少；[σν]和[Cb2]是常數；ν是分子運動粘性系數。對于本文所研究的飛行器模型，屬于雷諾數很大，黏度很小，有邊界層的湍流，并且對網格劃分要求不高，因此選用單方程模型進行模擬運算。

2 計算方法與過程

用Fluent軟件包的前置處理器GAMBIT軟件，按照由點到線，由線到面，由面到體的原則對飛行器進行建模，并根據尺寸比例在外圍設置一個繞流流場，繞流流場設置過程為：外圍是一個半徑為7 200 mm，高為3 600 mm的大圓柱形，用來作為飛行器的繞流流場。大圓柱內部設了一個小圓柱體，其半徑為1 800 mm，高為3 600 mm（飛行器位于計算域的正中），并進行較密的網格劃分，兩個圓柱之間的區(qū)域可將網格劃分得稀疏些，這樣既保證了計算的準確性，又減少了計算量。飛行器的舵翼部位采用密集的網格劃分，以反映氣流的劇烈變化，頭部和彈體采用兩頭密、中間稀的網格劃分；由于飛行器是軸對稱的，為了提高計算效率，減少計算時間，取整個模型的一半進行計算。整個區(qū)域的網格有2 018 432個。

對飛行器的相關區(qū)域進行網格劃分，如圖1所示。

把生成的網格文件導入FLUENT求解器里進行求解，具體的求解過程如下：

（1） 粘性模型采用單方程模型，使用運輸方程求解；

（2） 設置流體材料屬性：密度為ideal?gas，材料設為air，在“Viscosity”一項中選Sutherland，采用三系數方法[8]，對于高速可壓縮流動氣體采用描述氣體粘度的Suther?land 定律較比較合適；

（3） 壁面條件：無滑移，壁面粗糙度選為0.5，所有其他標量選擇不可滲透壁面條件；

（4） 數值計算過程中的差分格式選擇：動量、湍流動能、湍流耗散率均選用second order upwind scheme，即二階迎風格式；壓力插值選用默認的standard 方法[9]；

（5） 邊界條件設置：設置大氣壓為101 325 Pa，馬赫數為0.8，攻角為4°，粘性比為10；

（6） 松弛因子的設置：求解時設置壓力項松弛因子為0.8，密度、質量力項為1，動量項為0.5，湍動能項為0.6，耗散率項為0.6，湍流粘性項為0.6；

（7） 設置迭代次數6 500次，求解。

圖1 計算域網格劃分

3 計算實例

在給定的攻角和馬赫數下，阻力系數、升力系數和俯仰力矩系數隨著迭代次數的增加而也跟著振蕩，一直到迭代6 000次以后，曲線的變化逐漸趨于穩(wěn)定，待曲線穩(wěn)定以后，就可以在圖中讀出阻力系數、升力系數、俯仰力矩系數等相關參數的值。圖2為阻力系數隨迭代過程變化曲線，圖3為升力系數隨迭代過程變化曲線，圖4為俯仰力矩系數隨迭代過程變化曲線。從圖2～圖4可以看出，其氣動特性是收斂的，表明所設計的系統是穩(wěn)定的。仿真出的氣動參數和在該攻角和馬赫數下的試驗參數表對比之后有誤差，誤差為0.05%，在誤差允許范圍內，因此可以驗證出仿真數據的準確性。這些氣動參數將為下一步實驗提供數據。

圖2 阻力系數隨迭代過程變化曲線

圖3 升力系數隨迭代過程變化曲線

圖4 俯仰力矩系數隨迭代過程變化曲線

如果采用湍流模型中的其他模型，比如雙方程模型，計算過程相對復雜，對計算網格的要求更加精確，即壁面條件要求更高，一旦網格劃分不符合要求，則計算時很難得到收斂曲線，并且在網格劃分準確的情況下計算次數要達到10 000萬次，而在實際項目中，不需要對飛行器的網格進行精確的劃分，而單方程模型簡單、計算次數少，對網格的要求不敏感，有效計算出了實驗所需的相關參數，大大地節(jié)省了實驗時間。

4 結 論

在給定攻角和馬赫數的前提下，采用湍流模型中的單方程模型，借助CFD軟件Fluent對某飛行器進行氣動力計算，可以得出飛行器的升力系數、阻力系數、俯仰力矩系數。這些參數將為下一步飛行器控制系統設計提供依據。

參考文獻

[1] 趙洪章，岳春國，李進賢.基于Fluent的飛行器氣動特性計算[J].彈箭與制導學報，2007，27（2）：203?205.

[2] 于勇.FLUENT入門與進階教程[M].北京：北京理工大學出版社，2008.

[3] 劉明侯.計算流體和傳熱傳質[D].合肥：中國科技大學，2005.

[4] 龐英良，宋衛(wèi)東，楊曉霖.基于Fluent 的末制導炮彈初始段氣動仿真[J].兵工自動化，2009，28（2）：13?15.

[5] 錢翼稷.空氣動力學[M].北京：北京航空航天大學出版社，2009.

[6] 尹志林.某隱身巡航導彈氣動及雷達目標特性分析[D].南京：南京航空航天大學，2009.

[7] 張偉，張安堂，肖宇.基于坐標旋轉數字計算方法的三維坐標變換[J].探測與控制學報，2011，33（2）：73?76.

[8] 朱銳，董二寶，張杰，等.頭部可偏轉飛行器氣動仿真與外形優(yōu)化[J].機械與電子，2008（8）：6?8.

[9] LANDERS M G， HALL L H. Deflectable nose and canard controls for a fin?stabilized projectile at supersonic and hypersonic speeds， AIAA 2003?3805 [R]. Orlando， Florida， USA： AIAA， 2003.