張輝, 李杰, 韓杰

(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安　710072)

?

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

張輝, 李杰, 韓杰

(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安710072)

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

關(guān)鍵詞：非定常流動；跨聲速流動；馬赫數(shù)；氣動力；壓力分布；計算效率

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

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

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

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

關(guān)于結(jié)構(gòu)振蕩幅值線化非線性跨聲速小擾動方程得到時間線化的跨聲速小擾動速勢方程如下：

(1)

(2)

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

(3)

(4)

(5)

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

(6)

得到

(7)

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

(8)

2定常流動數(shù)值解的引入

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

(9)

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

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

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

3.1AIC矩陣的生成

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

(10)

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

(11)

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

(12)

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

(13)

3.2塊三對角近似技術(shù)

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

(14)

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

(15)

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

3.3針對復(fù)雜構(gòu)型的重疊場源策略

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

4算例驗證

4.1F5機翼

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

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

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

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

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

4.2CRM WBH翼身尾構(gòu)型

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

圖7　CRM WBH構(gòu)型場源模型

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

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

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

5結(jié)論

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

參考文獻：

[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(xué)計算[J]. 北京：北京航空航天大學(xué)出版社，1991

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

[6]徐敏，安效民，陳士櫓．一種CFD/CSD耦合計算方法[J]．航空學(xué)報，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

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

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

中圖分類號：V215

文獻標(biāo)志碼：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