基于速度分裂法的翼型陣風響應及減緩數值模擬

高 坤，郭同慶，紀哲翰，周 迪，陸志良

(南京航空航天大學 非定常空氣動力學與流動控制工業和信息化部重點實驗室，南京 210016)

0 引 言

陣風載荷響應作為一種氣動彈性動力響應問題，是指飛行器在大氣中飛行時遇到突然的陣風擾動，這一擾動氣流將改變飛機氣動升力面的有效迎角，引起升力等的突然改變，使得飛行器產生涉及彈性變形的陣風載荷響應。當飛行器處于陣風環境中時，飛行器結構將承載較大的陣風載荷，這將導致飛行器操縱困難且容易發生結構疲勞損壞。隨著航空科學技術的發展，大型飛行器的結構設計朝著低重量和大柔性方向發展，這導致陣風對大型飛行器的飛行影響更為復雜[1]。因此，對于現代飛行器設計，其陣風響應分析的非定常氣動模型由傳統線化理論逐漸轉向更為準確的CFD技術，其中陣風效應的嵌入主要分為遠場邊界條件法和給定流場陣風速度法兩類。

將陣風加入CFD程序最直接的方法是修改計算域遠場邊界條件，稱為遠場邊界條件（far-field boundary condition, FBC）法[2]。該方法為防止因計算域網格分辨率不足引起數值耗散導致陣風響應計算精度下降，要求計算域網格從物體表面到遠場邊界均需要非常精細，這將導致高昂的計算成本。為此，Tang和Bader[3]提出根據陣風在計算域中所處位置的變化對計算網格進行局部細化，該方法采用高階多項式插值來重新分布網格點以獲得陣風區域更新后的精細網格[4]。

給定速度法根據每一時刻陣風場在計算域中的空間位置來計算陣風速度，因此不存在遠場邊界條件法中因遠離物面區域網格分辨率不足而引起的數值耗散問題，被廣泛應用于陣風響應數值計算。給定速度法包括網格速度法（field velocity method, FVM）和速度分裂法（split velocity method, SVM）。其中，FVM由Parameswaran和Bader提出[5]，該方法借用動網格系統流動控制Euler/N-S方程，將陣風速度等效為流場當地網格速度來模擬陣風運動，但實際并不會引起網格運動。他們采用FVM模擬翼型穿越階躍陣風時的氣動響應，計算結果與理論分析結果十分接近。然而，FVM并非通過嚴格理論推導得來[6]，也沒有考慮翼型對陣風的影響。顧寧等[7]指出當FVM應用于N-S方程時需要加入黏性通量修正。相比而言，由Wales等[8]提出的SVM則是以非定常Euler/N-S方程為基礎，將陣風條件下的流場速度分解為陣風速度與無陣風作用下的流場背景速度兩部分，嚴格理論推導出SVM陣風模擬控制方程。與FVM相比，除陣風速度外，SVM還包括了陣風對黏性通量的影響以及與陣風時空導數相關的源項。此方法已被證明在較短陣風尺度、跨聲速等條件下的陣風響應模擬中能夠獲得更為準確的計算結果[8]。

目前，國內還尚未開展基于SVM的飛行器陣風響應模擬和減緩研究。本文針對One-Minus-Cosine離散陣風，首先基于動態網格系統非定常N-S方程，理論推導SVM陣風響應模擬的控制方程。然后耦合翼型結構運動方程，發展基于SVM的彈性翼型陣風響應CFD/CSD時域耦合算法。進而采用控制翼型旋轉，發展彈性翼型陣風載荷減緩的數值模擬方法。最后采用上述方法，開展剛性及彈性翼型陣風響應和減緩分析。

1 陣風模型

本文針對One-Minus-Cosine型離散陣風開展翼型陣風響應分析。飛行器以速度V水平飛行，突然受到速度為wg的垂向陣風作用。wg為陣風速度Vg在直角坐標系下的垂向分量，定義如下[9]：

其中：Lg為陣風尺度；wg=Vtan(Δα)為陣風速度幅值；Δα為陣風速度幅值等效迎角。

2 陣風響應CFD模擬方法

2.1 網格速度法（FVM）

本節簡要給出FVM[10-11]的控制方程，便于與SVM方法進行比較。動態網格系統下非定常N-S方程為：

其中：? 、為背景流場速度的速度分量；xt、zt為網格速度Vt的分量；ρ為流體密度；E為單位質量流體總能；p為流體壓強；γ為比熱比。

根據相對運動原理，如果計算域網格的運動速度為Vt，則等價于在計算域靜止情況下受到速度為-Vt的氣流作用[12]。因此，基于動態網格系統非定常N-S方程（2），FVM將陣風速度表示為相反方向的網格速度來模擬陣風作用[13]，即：

需要注意的是，FVM只是利用了網格速度的概念，在實施過程中并不引起實際的網格運動。盡管FVM是目前陣風響應CFD分析的主要方法，但該方法并未得到嚴格的理論闡明。

2.2 速度分裂法（SVM）

SVM將陣風條件下的流場速度分解為陣風速度與流場背景速度兩部分，即

其中，u、w為流場總速度V的速度分量。為了開展考慮氣動彈性的陣風響應分析，接下來將現有的基于固定網格系統的SVM方法[14]拓展至動態網格系統。包含陣風作用的動態網格系統下非定常N-S方程可寫為：

聯合式（5）～式（7），并利用背景流場的連續方程和動量方程，可得SVM的控制方程為：

其中黏性應力張量的表達式為：

其中，μ為黏性系數。式（8）中附加源項為：

若采用Euler方程，黏性通量為0，則附加源項為：

Sm、Se分別定義如下：

不難看出，相比于SVM，傳統FVM的動量方程和能量方程缺少了上述附加源項。該源項與陣風場的時間和空間導數相關，反映了背景流場對陣風場的影響。對于小尺度陣風或存在激波情形，附加源項的作用將會凸顯。

本文離散陣風形式由式（1）給定，因此源項中陣風相關的導數可以解析求解。采用有限體積法、雙時間推進法[15]求解包含陣風作用的動態網格系統下的非定常N-S方程（8），采用并行算法提高計算效率；基于RANS模擬湍流，湍流模型選用SA一方程模型[16]。遠場采用無反射邊界條件，物面采用黏性無滑移邊界條件，即本文后續彈性翼型陣風響應及載荷減緩計算均采用無限插值法（TFI）[17]生成動態網格。

3 彈性翼型陣風響應CFD/CSD時域模擬方法

3.1 三自由度結構運動方程

圖1給出了一個具有沉浮、俯仰和控制面偏轉的三自由度典型二元翼段力學模型。圖中三個自由度為：剛心E的沉浮位移h，向下為正；翼段繞剛心的俯仰角位移α，抬頭為正；控制面偏轉角 β ，向下偏轉為正。彈性軸E距離翼弦中點的距離為ahb，ah為基于半弦長b的無量綱距離，彈性軸在翼弦中點之后時ah＞0。Cβb表示控制面鉸鏈軸距翼型弦線中點的距離，其氣動彈性運動方程為[18-21]：

圖1 三自由度典型翼段力學模型[21]Fig.1 Sketch of an airfoil with three degrees of freedom[21]

其中：m為翼段質量；Iα、Iβ分別為翼段繞剛心和控制面繞鉸鏈軸的轉動慣量；Sα、Sβ分別為翼段和控制面的質量靜矩；Th、Tα和Tβ分別為翼段的沉浮、俯仰和控制面轉動阻尼系數；Kh、Kα和Kβ分別為翼段的沉浮、俯仰和控制面轉動剛度；L、Mα和Hβ分別為升力、俯仰力矩和鉸鏈力矩；CL、Cm和CH分別為升力系數、俯仰力矩系數和鉸鏈力矩系數； ρ∞為自由來流密度；c為翼型弦長。

式（16）的無量綱形式為：

為了便于在時域中求解結構運動方程（18），將其轉換為狀態空間形式：

其中，I為單位矩陣。

本文采用隱式歐拉法求解方程（20）：

得到n+ 1時刻的狀態變量：

3.2 CFD/CSD時域耦合計算流程

給定計算狀態，首先不考慮陣風作用，CFD計算獲得該狀態下的定常流場，然后基于定常流場開展陣風響應時域計算。基于定常流場，引入陣風場，在時域內耦合求解氣動方程和結構運動方程，從而求解出陣風作用下翼型氣動響應和結構位移響應。

4 基于翼型俯仰運動的陣風減緩

類似旋轉翼尖，本文嘗試通過控制翼型俯仰運動來進行陣風減緩。將翼型的俯仰控制角度αc與其位移響應、速度響應和加速度響應聯系起來，

其中，G1、G2、G3、G4、G5、G6為控制系統增益。

則 式 （18） 右 側 的 σa應 改 寫 為 σ =σa+σc， 其 中為控制力向量。

5 算例分析

5.1 剛性翼型陣風響應計算

分別采用FVM和SVM方法，開展NACA0012剛性翼型的One-Minus-Cosine陣風響應計算分析。計算條件為：Ma= 0.3，迎角 0°，Re= 7.0 × 106；力矩參考點位于0.4c；定義無量綱時間S=2V∞t/c。采用C型結構網格，陣風從遠場邊界外進入計算域。根據文獻[8]中不同網格數對比計算結果，翼型周向和法向網格數取為257 × 129。翼型弦長c= 7.5m，遠場邊界位于15倍弦長處，第一層網格高度為1×10-3c，如圖2所示。本節除黏性對陣風響應影響分析采用N-S方程計算之外，其余算例均采用Euler方程。

圖2 NACA0012翼型計算網格Fig.2 Computational mesh for NACA0012

5.1.1 時間步長的確定

陣風尺度Lg= 1c，等效迎角 Δα= 2.0°。分別取時間步長 ΔS為 0.1、0.01和 0.001，計算結果如圖3所示。對于FVM和SVM， ΔS= 0.1時的升力系數最大值和力矩系數最大值均大于 ΔS= 0.01或 ΔS= 0.001時的最大值，而 ΔS= 0.01與 ΔS= 0.001的結果相吻合，因此本文時間步長取 ΔS= 0.01。

圖3 時間步長對FVM和SVM陣風響應的影響Fig.3 Effects of time step on gust response obtained by FVM and SVM

圖4比較了陣風尺度為1c時本文計算得到的升力系數陣風響應與文獻[8]的結果。對于小尺度陣風，陣風梯度增大、作用域更為集中，在陣風離開翼型后緣的時間段內，升力系數產生小幅波動；除波動幅值存在一定差異外，本文FVM和SVM兩種方法計算得到的升力系數峰值及其出現時刻均與文獻結果吻合良好。

圖4 FVM和SVM氣動力陣風響應與文獻結果對比Fig.4 Comparisons of aerodynamic responses between present and ref.[8]data

進一步研究表明，陣風離開翼型后緣時升力系數小幅波動的幅值與翼型后緣厚度的處理方式直接相關。圖5比較了尖后緣和鈍后緣NACA0012翼型計算得到的升力系數陣風響應，鈍后緣翼型后緣厚度為1.26‰c。相比于尖后緣，陣風離開鈍后緣翼型時升力系數響應的波動幅值更大，峰值略微降低。

5.1.2 陣風尺度的影響

分別采用FVM、SVM計算陣風尺度 1c、5c和10c，等效迎角 Δα= 2.0°的陣風響應，計算結果如圖6所示。在陣風尺度5c和10c時，FVM和SVM的計算結果幾乎一致。當陣風尺度減小到1c時，陣風場的時空導數增大，SVM方法中源項的作用凸顯出來，導致FVM和SVM的計算結果出現差異，SVM計算得到的氣動參數陣風響應峰值略高于FVM。SVM方法的推導過程表明FVM只是SVM忽略源項的一種近似方法，因此理論上SVM的計算結果更為準確。隨著陣風尺度的增加，陣風場的時空導數減小，陣風離開翼型后緣時氣動系數的小幅波動現象消失。

圖6 陣風尺度對FVM和SVM的陣風響應影響Fig.6 Gust length effect on gust responses obtained by FVM and SVM

5.1.3 黏性對陣風響應的影響

針對陣風尺度1c和10c、等效迎角 Δα= 2.0°，采用SVM方法開展Euler、N-S方程的陣風響應對比計算，其中N-S方程計算網格的第一層網格高度為1×10-5c，其他網格參數與Euler方程計算網格相同。計算結果對比如圖7所示，黏性效應導致陣風響應的升力系數峰值降低，陣風尺度越大，黏性對陣風場的影響越大，升力減小量越明顯；黏性作用下，1c陣風離開翼型后緣時的升力系數波動明顯減弱。

圖7 黏性對陣風響應的影響Fig.7 Viscosity effects on gust response

5.2 彈性翼型陣風響應計算

采用非定常Euler方程計算彈性NACA0012翼型的陣風響應。翼型結構模型參數[22]由表1給出。

表1 NACA0012翼型結構參數Table 1 NACA0012 parameters

5.2.1 陣風響應CFD/CSD耦合算法驗證

選取文獻[22]的計算參數，Ma= 0.3，V*= 0.2，翼型弦長c= 1.0 m。翼型結構位移的初始狀態 ξ (0) = 0，α(0)= 0°，陣風尺度Lg= 6.25 m，等效迎角 Δα= 0.57°。考慮到文獻計算方法為FVM，本算例也采用FVM。計算結果對比如圖8所示。本文計算結果與文獻[22]數據吻合，驗證了本文彈性翼型陣風響應CFD/CSD耦合算法的可靠性。

圖8 彈性翼型陣風響應與文獻結果對比Fig.8 Comparison of gust responses of an elastic airfoil between present and reference results

5.2.2 彈性對翼型陣風響應的影響

采用SVM計算陣風尺度1c和10c、等效迎角 Δα=2.0°的陣風響應，V*= 0.2，翼型結構初始狀態 ξ (0) =0，α(0)= 0°。圖9為剛性翼型和彈性翼型的氣動系數陣風響應對比。圖10為彈性翼型的結構響應結果。在One-Minus-Cosine陣風穿越彈性翼型期間，翼型的氣動參數和俯仰角均先增大然后減小至極小值，沉浮位移則一直增大至峰值；10c陣風穿越彈性翼型時，升力和俯仰角近乎同時達到峰值，導致升力峰值較剛性翼型增大；1c陣風穿越彈性翼型時，俯仰角峰值遠滯后于升力峰值，彈性對升力峰值幾乎沒有影響。

圖9 剛性和彈性翼型氣動力陣風響應對比Fig.9 Comparison of aerodynamic gust response between rigid and elastic airfoils

圖10 彈性翼型結構響應Fig.10 Structural response of an elastic airfoil

5.3 彈性翼型陣風載荷減緩

以NACA64A010翼型為研究對象，翼型弦長c=1.0 m，其他網格參數與第5.1節Euler方程計算網格相同，如圖11所示。翼型結構參數[23]如表2所示，取Ma= 0.3、V*= 0.337 4。采用SVM進行翼型陣風載荷減緩響應計算。陣風尺度Lg= 10c，等效迎角 Δα=2.0°。翼型結構位移 的 初 始狀態 ξ (0) = 0、α(0) = 0°、β(0)= 0°。本文以減緩翼型在通過陣風時產生的最大升力系數為目標。

表2 NACA64A010翼型結構參數Table 2 NACA64A010 parameters

圖11 NACA64A010翼型計算網格Fig.11 Computational mesh for NACA64A010

彈性翼型陣風響應分析表明，在陣風穿越翼型期間，升力系數先增大后減小至極小值，而翼型沉浮位移則一直增大至極值，因此本文選取沉浮速度值作為控制輸入量。圖12、圖13為增益G3分別取5.0和10.0時閉環控制系統與開環系統陣風響應的結果對比。采用G3閉環控制情況下，沉浮運動峰值顯著降低并迅速收斂；在升力系數峰值之前，由于翼型沉浮速度為負值，G3等效于給彈性翼型增加低頭俯仰控制力矩，導致俯仰角、升力和力矩系數峰值明顯降低，G3= 10時升力系數峰值降低49.758%；當陣風離開翼型后，翼型沉浮速度為正值，俯仰控制力矩使翼型俯仰角逐漸增大至正極大值點，并高于無控制情況，導致此時氣動系數響應極值也均略高于無控制情況，并接近或高于陣風穿越翼型時的峰值。

圖12 開環系統和閉環控制系統翼型氣動力陣風響應對比Fig.12 Comparison of aerodynamic gust responses between open- and closed-loop control systems

圖13 開環系統和閉環控制系統翼型結構響應對比Fig.13 Compsrison of structural response between open- and closed-loop control systems

接下來采用G2= -5.0、G3= 10.0同時控制翼型的沉浮和俯仰運動。圖14、圖15分別為G2、G3組合的閉環控制陣風響應結果和僅G3的響應結果對比。增加俯仰控制后，俯仰角位移幅值顯著減小，導致陣風作用期間的氣動參數峰值高于無俯仰控制情況，但陣風離開翼型后的氣動參數極值則顯著降低。

圖14 G2 = -5.0、G3 = 10.0和僅有G3 = 10.0氣動力響應對比Fig.14 Comparison of aerodynamic response between G2 = -5.0, G3 = 10.0 and only G3 = 10.0

圖15 G2 = -5.0、G3 = 10.0和僅有G3 = 10.0結構響應對比Fig.15 Comparison of structural response between G2 = -5.0, G3 = 10.0 and only G3 = 10.0

6 結 論

本文基于動態網格系統下非定常N-S方程，理論推導出SVM陣風模擬的控制方程，發展出基于SVM的彈性翼型陣風響應CFD/CSD時域耦合算法和基于翼型俯仰控制的陣風減緩模擬方法。算例分析了NACA0012剛性、彈性翼型的One-Minus-Cosine陣風響應，計算結果與文獻數據吻合；分析了陣風尺度、流體黏性和結構彈性對NACA0012翼型陣風響應的影響；實現了基于俯仰控制的NACA64A010彈性翼型陣風減緩。主要結論如下：

1） 理論推導表明，傳統FVM是忽略SVM中陣風相關附加源項的一種近似方法，對于大尺度陣風，兩種方法的計算結果一致，源項的作用隨著陣風尺度的減小而凸顯，理論上SVM更為準確；

2） 黏性效應導致NACA0012翼型升力系數陣風響應峰值降低，陣風尺度越大，升力減小越明顯，黏性同時減弱了1c陣風離開翼型后緣時的升力系數波動；

3） 結構彈性引起NACA0012翼型俯仰角陣風響應峰值增加，陣風尺度越大，升力系數增加越顯著；

4） 俯仰控制能夠有效減緩NACA64A010彈性翼型的陣風載荷，沉浮速度控制輸入量G3= 10.0時升力系數陣風響應峰值減緩49.758%，結合G2、G3則能夠同時減緩沉浮和俯仰響應。