基于重疊場源法的非定常氣動力計算研究

張輝, 李杰, 韓杰

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

?

基于重疊場源法的非定常氣動力計算研究

張輝, 李杰, 韓杰

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

摘要：面向氣動彈性工程應用，將基于雷諾平均N-S方程的定常流動數值解引入跨聲速小擾動方程，在小擾動條件下非定常激波效應由定常激波效應確定，從而發展了一種可用于跨聲速流動的非定常氣動力計算方法。數值算法利用塊三對角近似技術大大提高了算法的計算效率、節省了計算所需內存空間，并采用了重疊場源策略為復雜構型的非定常氣動力計算提供了有利保障。文中以F5機翼為計算算例，研究了場源模型參數對非定常氣動力計算結果的影響，并驗證了塊三對角近似算法的可靠性；以CRM WBH翼身尾構型為算例，檢驗了針對復雜構型重疊場源策略的可行性。

關鍵詞：非定常流動；跨聲速流動；馬赫數；氣動力；壓力分布；計算效率

非定常氣動力的計算對氣動彈性分析、動導數計算及機動、突風載荷計算等都具有特別重要的意義[1-3]。長期以來，廣泛應用于氣動彈性分析的氣動力計算方法是偶極子格網法[4-5](doublet lattice method,DLM)，多個商用軟件均采用了DLM進行氣動彈性分析，如MSC.NASTRAN的氣動彈性模塊。DLM方法能夠提供亞聲速及超聲速范圍內較為精確的非定常氣動力，其最大的優點還在于能夠生成氣動力影響系數(AIC)矩陣。AIC矩陣不依賴于飛機結構參數，因此，在結構設計循環中只需計算1次而重復使用。但是DLM是基于線化勢流理論的方法，無法解決非線性較強的流場，不適用于跨聲速非定常氣動力計算。計算流體力學(CFD)方法通過求解Euler或Navier-Stokes方程能夠給出精確的跨聲速流場解，不過基于CFD的流場求解不能生成AIC矩陣，因而不能利用目前已經發展較完備的基于AIC矩陣的氣動彈性分析方法。而目前用于氣動彈性分析的CFD/CSD耦合計算方法[6-7]針對復雜飛機構型的適用性還有待發展，另外耦合算法需要巨大的計算資源和時間耗費，這樣均不利于工程應用。因此，發展一種能夠生成AIC矩陣的高效非定常跨聲速氣動力計算方法對跨聲速范圍內飛行器的氣動彈性特性研究意義重大。

非定常場源方法即是一種能夠生成AIC矩陣的方法，并且通過引入CFD定常解考慮了跨聲速激波效應，從而適用于非定?？缏曀贇鈩恿Φ挠嬎?。早在20世紀80年代，國外一些研究者就已經開始了對非定常場源法的研究，Larry L.Erickson和Shawn M.Strande[8]給出了利用場源法將面元方法推廣以求解跨聲速流動的理論基礎；M.H.L Hounjet[9]研究了基于場源法和有限差分的混合方法求解非定常跨聲速流動問題，該方法綜合了有限差分和場源積分方法的優點，從而大幅地減少了計算時間，使得其可用于常規的顫振分析；Lutz Gebhardt和Dmitri Fokin等[10]研究了用于跨聲速飛行器氣動設計的場源方法，改進了自適應場源網格生成方法使得場源法的實際應用更加現實；Chen和Gao等[11]基于重疊場源策略通過求解關于時間線化的跨聲速小擾動方程的積分方程生成非定?？缏曀贇鈩恿τ绊懴禂稻仃?，從而實現復雜構型的跨聲速氣彈分析，并利用塊三對角近似技術求解大型的體單元影響系數矩陣，大大地提高了計算效率。然而，國內對基于場源法的跨聲速非定常氣動力計算的研究目前尚無公開發表的文獻。

本文深入研究了非定常場源法的算法實現，將基于雷諾平均N-S方程的定常流動數值解引入跨聲速小擾動方程，在小擾動條件下非定常激波效應由定常激波效應確定，從而發展了一種可用于跨聲速流動的非定常氣動力計算方法。

1非定常場源法的求解方程

關于結構振蕩幅值線化非線性跨聲速小擾動方程得到時間線化的跨聲速小擾動速勢方程如下：

(1)

(2)

φs表示表面面元強度對速度勢的影響，可寫為

(3)

(4)

(5)

式中，Δσv表示穿過激波面體源強度的跳躍。當不存在跨聲速激波時，Δσv=0,φshock自動消失。另外，當激波出現時，也可通過對φv完成如下分部積分運算以消掉φshock。

(6)

得到

(7)

式中，xs表示激波位置，ε表示激波面無限小厚度。由(5)式和(7)式可得

(8)

2定常流動數值解的引入

本文用于求解定常流動數值解的控制方程如下

(9)

式中，Q=(ρ,ρu,ρv,ρw,ρe)T為守恒向量；ρ、e、(u,v,w)分別為密度、單位質量氣體的總能量和直角坐標系下的速度分量；v為控制體體積；s為控制體表面積；F為通過表面s的黏性通量和無黏通量的和；n為邊界外法向單位矢量。黏性項采用中心差分格式離散，無黏項采用Roe三階迎風偏置通量差分裂方法離散，采用隱式近似因子分解(AF)方法進行時間推進，選用SA湍流模型，通過多重網格技術來加速收斂。

采用場源法求解非定常氣動力時，需要建立場源模型，場源模型包括面單元和體單元。首先建立面元氣動模型，通常將飛機部件分為翼面類和機身類部件處理，翼面類部件由位于翼面類部件均平面的非定常渦面來模擬，機身類部件表面離散為體表面單元，每個體表面單元上布置非定常源，以模擬由于體的體積效應產生的氣動力分布；其次圍繞翼面類部件(升力面)或機身類部件的面元氣動模型定義體塊，然后將體塊分割為若干體單元。場源模型建立后，即可將上述定常流場數值結果(包括當地馬赫數及擾動速度分量)插值到每個體單元上。建立合適的場源模型直接決定著插值精度，進而影響最終非定常氣動力計算的精度。在建立場源模型時，確保場源模型中的物面邊界盡可能與CFD計算模型的物面邊界重合是建立合適的場源模型的必要條件。此外，體塊高度和體塊分割層數是場源模型的2個重要參數，這2個參數的選取對計算結果有重要影響，選擇的依據是確保場源模型包含非線性流動區域。算例驗證部分詳細研究了體塊高度和體塊分割層數對計算結果的影響，并歸納總結了體塊高度和體塊分割層數選取的一些原則。

3AIC矩陣的生成及算法策略

3.1AIC矩陣的生成

位于面單元控制點處的法向擾動速度為

(10)

體單元控制點處的速度勢為

(11)

式中，矩陣[A]和[C]分別為面奇點對面單元和體單元的影響系數矩陣，矩陣[B]和[D]分別為體源對面單元和體單元的影響系數矩陣。引入如下有限差分算子[T]

(12)

將方程(12)代入(10)式和(11)式，得到

(13)

3.2塊三對角近似技術

上節中矩陣[E]為滿系數復數矩陣，其階數與體單元數相同。對于簡單的三維問題需要較多體單元數，并且體單元數隨構型復雜程度的增加而急劇增加。這種大型矩陣求逆對常規的氣彈及載荷分析是不實用的。因此，有必要采用塊三對角近似技術解決大型滿系數矩陣[E]求逆的問題。塊三對角近似技術的思路[11]是：首先體單元被分為許多子塊，同一個子塊內的體單元組合在一起。這樣矩陣[E]可寫為

(14)

式中，[EB]為塊三角矩陣，其三角塊包含自身塊和相鄰子塊的影響系數；[Eε]包含三角塊處的零元素和非相鄰塊的影響系數。下一步比較矩陣[Eε]和矩陣[EB]內系數的量級，可以看出[Eε]內所有系數均是小量，原因是方程(4)的積分方程的積分函數包含1/R函數，當點(x0,y0,z0)遠離非相鄰子塊時，1/R快速衰減。因此，矩陣[E]的逆可由下式近似

(15)

因為[EB]為塊三角矩陣，所以[EB]-1可利用塊三角矩陣求解技術有效計算。

3.3針對復雜構型的重疊場源策略

針對復雜飛機構型，本文采用重疊場源策略以減小網格生成的難度，減少計算所需內存并節省計算時間。首先針對復雜構型的每一個部件獨立生成體單元網格，然后基于以下原則[11]構建這種重疊體單元模型的影響系數矩陣：①不同體塊內的體單元互不影響；②同一個體塊內的體單元僅影響與其相關聯的面元；③所有面元影響所有體單元。

4算例驗證

4.1F5機翼

本節以F5機翼[12]為算例，建立不同參數下的場源模型，取計算馬赫數為0.90，減縮頻率為0.275，進行計算分析并與實驗值對比以研究場源模型不同參數對非定常氣動力計算精度的影響。主要研究了場源模型的體塊高度H、體塊分割層數Nl、體塊前(后)緣向前(后)延伸量EXT及分割數Ne對非定常氣動力計算精度的影響。以機翼根弦長C為參考，分別取體塊高度0.25C、0.50C、1.00C和1.50C，保持其他參數不變，其中體塊層數為7。計算所得展向51.5%位置處非定常壓力系數結果與實驗值對比如圖1所示，當體塊高度較小，如H=0.25C和0.50C時，計算結果明顯小于實驗值，計算對激波強度的預測明顯偏??；隨著體塊高度的增加，在H=1.00C時，計算結果與實驗值比較接近，計算精度較好；隨著體塊高度的繼續增加，在H=1.50C時，計算結果則高于實驗值，顯然這不是一個理想的結果，因為這與之前的理論分析相悖，按照理論分析當H=1.50C時其場源模型比其他場源模型大，所以合理的結果應該是此時的計算結果應更加接近實驗值，但圖2結果并非如此。進一步，保持體塊高度為H=1.50C，選取體塊層數分別為5、7、9、11完成計算，同樣給出展向51.5%位置處非定常壓力系數的計算結果與實驗值對比如圖2所示，對比結果表明，隨著體塊層數的增加，計算結果不斷逼近實驗值，并且Nl=11對應的計算結果略好于圖1中體塊高度H=1.00C的計算結果。因此，可以得到如下結論：①增大場源模型以確保場源模型包含非線性流動區域可以提高非定常氣動力的計算精度；②若保持其他參數不變，只是增大體塊高度會使非定常氣動力的計算精度變差，原因是隨著體塊高度的不斷增加，使得體單元分布越來越稀疏，過于稀疏的體單元分布降低了其計算精度；③保持其他參數不變，增加體塊層數，計算結果趨于收斂；另外，圖3給出了一個狀態所需的CPU計算時間隨體單元個數的變化曲線，可見CPU計算時間隨著體單元數的增加呈指數式增加。因此，為了兼顧精度與效率，對于該算例，選取體塊高度為1.0C，體塊層數為7較合適。

圖1　體塊高度的影響　　　　　　　圖2　體塊分割層數的影響圖3　計算時間隨體單元數的變化

取體塊高度為1.0C，體塊層數為7，分別選取體塊前(后)緣向前(后)延伸量EXT為0.00C、0.10C、0.25C、0.50C的場源模型進行計算，計算結果與實驗值對比如圖4所示，可見隨著延伸量EXT的增加，計算結果向實驗值靠近，但變化量很小。圖5給出了延伸段不同分割數的場源模型計算結果對比，可以看到分割數分別為2、4、6時，計算結果基本重合，表明該處分割數對計算結果幾乎沒有影響。

針對F5機翼，分別建立單塊場源模型和多塊場源模型，對于多塊場源模型，利用塊三對角近似技術進行求解，對比2種場源模型計算結果，驗證塊三對角近似的可靠性。圖6給出了機翼展向51.5%位置處壓力分布對比結果，可見三對角近似求解結果與單獨體塊求解結果基本重合，這表明了三對角近似求解是可靠的。此外，與單獨體塊求解相比，塊三對角近似求解節省了大約40%的CPU計算時間。

圖4　場源模型前后延伸量的影響　　　圖5　延伸段分割數的影響　　　圖6　塊三對角近似技術的驗證，k=0.275

4.2CRM WBH翼身尾構型

采用重疊場源策略建立CRM WBH飛機構型[13]場源模型，分別針對機翼、平尾、機身獨立生成體塊，每個體塊被分割為若干子塊生成體單元網格以應用塊三對角近似技術求解。根據上節研究結論，機翼和平尾對應體塊高度取各自的根弦長，體塊層數為7，機身對應體塊高度取為機身中段半徑長，體塊層數為8，建立場源模型如圖7所示。

圖7　CRM WBH構型場源模型

取計算馬赫數0.85，減縮頻率0.15進行計算。圖8給出了CRM WBH 構型俯仰振蕩模態時機翼展向50.2%處翼面非定常壓力分布對比結果。由于沒有實驗結果，本文以ZAERO軟件[14]的跨聲速方法計算結果作為參考驗證本文計算方法的可靠性，對比可見本文計算結果與ZAERO計算結果基本一致，僅在對激波強度的預測方面稍有差異，這說明了針對復雜飛機構型本文采用的重疊場源策略是可靠的。同時圖中給出了偶極子格網法的計算結果，顯然，由于偶極子格網法不能考慮跨聲速激波效應的影響，而對于當前計算狀態下機翼翼面展向50.2%位置處存在較強激波的狀況，偶極子格網法計算結果必然與考慮了跨聲速激波效應的計算方法所得結果存在較大差異。換言之，如果翼面不存在激波，也不存在流動分離，那么偶極子格網法計算結果就應該與本文氣動力方法計算結果接近。平尾翼面壓力分布對比結果有利地證實了這一點。圖9給出了俯仰振蕩模態時平尾展向50%處翼面非定常壓力分布對比結果。圖9中本文氣動力方法、偶極子格網法及ZAERO跨聲速方法計算結果可以看到，3種計算方法所得結果基本一致。

圖8　機翼展向50.2%處非定常壓力系數

圖9　平尾展向50%處非定常壓力系數

5結論

本文通過引入CFD定常解考慮跨聲速激波效應，面向氣動彈性分析發展了一種跨聲速非定常氣動力計算方法。文中以F5機翼為算例研究了場源模型參數對非定常氣動力計算結果的影響，得到了選擇合適場源模型的一般規律，并驗證了塊三對角近似算法的可靠性；以CRM WBH翼身尾構型為計算算例，單塊和多塊場源模型計算結果的對比表明針對復雜構型重疊場源策略是可行性的。算例計算結果與實驗及參考文獻結果的對比分析表明本文發展的方法適用于跨聲速非定常氣動力的計算。

參考文獻：

[1]Rafael N, Joao L. Comparison of Viscous and Inviscid Unsteady Aerodynamic Loads for Aeroelastic Analysis[R]. AIAA-2015-0767

[2]Hyungro Lee, Beom-Soo Kim, et al, Computational Study of the Stability Derivatives for the Standard Dynamic Model[R]. AIAA-2013-2658

[3]Leonard M, Ilhan T. Time Simulations of the Response of Maneuvering Flexible Aircraft[J]. Journal of Guidance, and Dynamics, 2004, 27(5): 814-828

[4]Albano E, Rodden W P. A Doublet-Lattice Method for Calculating Lift Distributions on Oscillating Surfaces in Subsonic Flows[J]. AIAA Journal, 1969, 7(2): 279-285

[5]管德. 非定?？諝鈩恿W計算[J]. 北京：北京航空航天大學出版社，1991

Guan De. Calculation on Unsteady Aerodynamics[M]. Beijing, Beijing University of Aeronautics and Astronautics Press, 1991 (in Chinese)

[6]徐敏，安效民，陳士櫓．一種CFD/CSD耦合計算方法[J]．航空學報，2006, 27(1)：33-37

Xu Min, An Xiaomin, Chen Shilu. CFD/CSD Couping Numerical Computational Methodology[J]. Acta Aeronautica et Astronautica Sinica, 2006, 27(1): 33-37 (in Chinese)

[7]Giulio Romanelli, Michele Castellani, Paolo Mantegazza, et al, Coupled CSD/CFD Non-Linear Aeroelastic Trim of Free-Flying Flexible Aircraft[R]. AIAA-2012-1562

[8]Erickson Larry L, Strande Shawn M. A Theoretical Basis for Extending Surface Paneling Methods to Transonic Flow[J]. AIAA Journal, 1985, 23(12): 1860-1867

[9]Hounjet M H L. A Field Panel/Finite Difference Method for Potential Unsteady Transonic Flow[J]. AIAA Journal, 1985, 23(4): 537-545

[10] Gebhardt Lutz, Fokin Dmitri, Lutz Thorsten, et al. An Implicit-Explicit Dirichlet-Based Field Panel Method for Transonic Aircraft Design[R]. AIAA-2002-3145

[11] Chen P C, Gao X W, Tang L. Overset Field-Panel Method for Unsteady Transonic Aerodynamic Influence Coefficient Matrix Generation[J]. AIAA Journal, 2004, 42(9): 1775-1787

[12] Malone J B, Sankar N L, et al. Unsteady Aerodynamic Modeling of a Fighter Wing in Transonic Flow[R]. AIAA-1984-1566

[13] John C, Edward N. Summary of the Fourth AIAA CFD Drag Prediction Workshop[R]. AIAA-2010-4547

[14] ZAERO Theoretical Manual[M]. ZONA Technology Inc Scottsdale, 2008

收稿日期：2015-10-22

基金項目：國家自然科學基金(11172240)、國家重點基礎研究發展計劃(2015CB755800)與航空科學基金(2014ZA53002)資助

作者簡介：張輝(1986—)，西北工業大學博士研究生，主要從事理論與計算流體力學研究。

中圖分類號：V215

文獻標志碼：A

文章編號:1000-2758(2016)03-0443-06

Calculation of Unsteady Aerodynamics Using Overset Field-Panel Method

Zhang Hui, Li Jie, Han Jie

(College of Aeronautics, Northwestern Polytechnical University, Xi′an 710072, China)

Abstract:The paper present a method of calculating unsteady aerodynamics in transonic flow, in which steady flow solution based on RANS solver code are imported into transonic small disturbance equations, thus, the unsteady shock location is dominated by the steady shock location in the small amplitude sense. Block-tridiagonal approximation is used to improve the computational efficiency and save the memory, and the overset field panel scheme is developed for complex configuration to calculate unsteady aerodynamics. Taking F5 wing as an example, the effect of parameters of field-panel model on the calculated results of unsteady aerodynamics is investigated, and the reliability of block-tridiagonal approximation is validated. In addition, the validation of overset field-panel scheme is performed by comparing the results between the one block modeling and the overset modeling.

Keywords:unsteady flow; transonic flow; Mach number; aerodynamics; pressure distribution; computational efficiency