不同流動狀態下舵面氣動熱彈性分析

范賽鋒，員海瑋

南京航空航天大學，江蘇 南京 210016

新一代高超聲速飛行器需要長時間保持高速飛行，受到的強烈氣動加熱使材料屬性和結構應力分布發生改變，進而導致結構的模態等特性發生重大改變。同時，高超聲速流場周圍存在強激波間斷、流動轉捩等復雜的流動。以上因素使得高超聲速飛行器的氣動彈性問題與常規飛行器存在明顯不同，其結構上承受非線性氣動力和強烈氣動熱的作用，多場耦合問題變得異常復雜，這些問題極大地增加了高超聲速飛行器的設計難度，因此高超聲速飛行器氣動熱彈性的研究具有非常重要的工程意義。

分析高超聲速氣動熱彈性問題時，對氣動熱流的處理，一般采用計算特定條件或者定常流動下的溫度場，計算結構在此溫度場下的傳熱，得到穩態的結構內部溫度分布，以此熱狀態作為新的結構屬性，進行氣動彈性計算。

張偉偉等[1]通過對氣動熱、結構溫度場和熱結構仿真的分步計算方法，在時域內研究了確定溫度分布條件下的氣動熱彈性分析。V.J.Shinde 等[2]通過直接數值計算方法（DNS）求解Navier-Stokes 方程，對柔性面板的過渡激波邊界層內的相互作用進行了數值模擬，結果表明柔性板相對于剛性板有更強的流固耦合作用。沈恩楠等[3]通過建立全時域的多物理場計算方法，采用有限體積法求解高超聲速條件下的流場和結構的溫度場，結果表明全時域耦合分析方法在模擬結構振動對流場和溫度場的影響方面，要明顯優于同步計算方法。高超聲速氣動熱彈性問題，區別于經典氣動彈性問題的關鍵就是需要考慮氣動熱和結構傳熱之間的耦合計算，詳細的介紹可以參見參考文獻[4]，本文限于篇幅不再贅述。

眾多學者在氣動彈性問題研究上已取得較多進展。劉成等[5]在假設舵面不同位置發生轉捩的條件下，構造出不同的熱分布模型，考慮熱應力和材料在不同溫度下的屬性。通過計算歐拉方程得到無黏流場，應用基于計算流體力學（CFD）技術的當地流活塞理論，計算高超聲速翼型的非定常氣動力進行氣動熱彈性分析，研究了轉捩位置的不同對舵面氣動熱彈性的影響。郭同慶等[6]采用單向氣動彈性—氣動熱耦合方法，開發了高超飛行器沿著預定的彈道飛行時的氣動熱彈性分析方法，通過CFD 和計算結構力學（CSD）耦合的方法，對每個軌跡點上瞬態加熱過程進行了顫振計算，從而給出了各狀態點的顫振邊界。N.Lamorte等[7]利用CFD 技術結合徑向基函數，建立了計算高超聲速下結構氣動熱彈性的框架，論證了從層流到湍流的過渡狀況，以及熱應力和高超聲速飛行器穩定性裕度之間的重要關系。Z.B.Riley 等[8]使用一個基本的二維氣動熱彈性模型，在簡支撐板假設下，研究了邊界層穩定性、氣動加熱和表面條件之間的關系，綜合研究結果表明，完全湍流條件下的計算結果，對預測表面板的熱結構相應并總趨于保守，轉捩條件下的熱彈性響應比全湍流條件下更加突出。趙仕偉等[9]提出一種基于氣動力降階模型的柔性后緣可變形機翼氣動彈性分析方法，在保證計算精度的同時，大幅提高氣動彈性分析效率。胡家亮等[10]針對某機翼變剛度低速顫振風洞試驗要求，采用3D 打印結構相似顫振模型方案，通過可拆卸蒙皮設計實現模型局部剛度可變，完成了模型結構設計、綜合優化、制造和風洞試驗，研究了機翼顫振耦合機理和局部剛度對顫振特性的影響，驗證了3D打印顫振模型設計方法的可行性和有效性。J.D.Thayer等[11]對按照規定運動的懸臂板進行了Ma2.0 流動下三個模擬分辨率上的研究，分別是非定常雷諾平均（URANS）、延遲渦模擬（DDES）和隱式大渦模擬（ILES），研究表明誘導載荷在各個分辨率的模擬方法上，預測是一致的，也就是分辨剪切層的中高頻動力學，可能并不是識別流體—結構耦合的必備條件。

在航空航天實際工程中，耦合求解URANS方程和非線性有限元的CFD/CSD方法使用廣泛程度并不高，國內的很多設計場所仍然依賴于基于勢流模型和結構模態的線化方法。本文采用成熟的Fluent 和ANSYS 求解舵面模型的非定常流場和非線性結構場的耦合問題，研究了多種條件對舵面氣動彈性的影響，同時探索了適用于多場耦合分析的高精度通用計算方法，為解決流固耦合分析中問題復雜性和易分析性之間的矛盾提供了參考。

1 流固熱耦合基本理論和控制方程

高超聲速氣動熱彈性問題是典型的多場耦合計算問題。本文采用載荷轉移法進行氣動熱—結構傳熱耦合計算，流體部分的實質是等溫壁面的流場計算，而結構部分是第二類熱傳導邊界條件下的熱平衡計算，同時需要進行壁面輻射模擬。輻射方程是四次非線性方程，因此需要進行反復迭代以獲得收斂的壁面溫度，作為邊界條件反饋給流場計算。載荷轉移法對初始值非常敏感，為此需要采用松弛系數控制每個迭代步加載到結構上的邊界熱流。

1.1 流體控制方程

舵面在流場中，遵循質量守恒、動量守恒、能量守恒三大物理守恒定律[12]，推導可以得到流體力學的基本控制方程，即連續性方程、動量方程和能量方程，具體方程可以很方便在相關文獻和計算流體力學書籍中找到，本文不再列出。

1.2 結構控制方程及耦合邊界條件

結構基本控制方程如下

式中，t為時間；M、C、K分別為結構的質量矩陣、阻尼矩陣和剛度矩陣；u為結構位移矢量；F為結構所受的載荷矢量，在這里為氣動力載荷。

為實現耦合計算，需要設置流場和結構的交界面，交界面上的邊界條件是

式中，σs為交界面上結構應力矢量；n是交界面上法矢量；p為交界面壓力矢量；σv為氣動力矢量；us、uf分別為結構位移矢量和流體位移矢量。即式（2）左式表示交界面上力的平衡，右式表示流體和結構的位移一致。

1.3 湍流模型

目前應用廣泛的湍流模型中，參照各湍流模型的適用范圍，選取Menter[13]提出的SST 二方程湍流模型。Menter在2004年對原SSTk-ω模型進行了改進，改進的SSTk-ω模型k和ω輸運方程如式（3）所示

式中，ρ為流體密度；Ui(i= 1,2,3)表示三個方向位移；xi(i=1,2,3)表示三個方向，同Ui對應；μ為流體黏度；μt為湍流動黏度；混合函數F1的定義如下。

湍流渦黏性定義如下

式中，S為應變率不變量；F2是第二混合函數，它的定義如下

SST 模型中同時添加了產生項限制器，以避免在停滯區域產生湍流堆積，具體如下

以上各式中常數項數值為

F.R.Menter用應變率不變量S替換了原SSTk-ω模型式（5）中的渦量，同時使用常數10替換了式（7）中原來建議的20。

1.4 動網格技術

本文將采用彈簧類比—超限插值方法（TFI）混合的方法實現動網格，其主要過程是：（1）將整個流場網格分為多個塊，將不同塊的角和點之間用彈簧連接，通過求解靜力平衡方程得到各個節點的位移；（2）對于內部的節點，采用TFI方法將邊界節點通過位移插值的方法分布到網格塊的內部節點上。下面簡要介紹彈簧類比法和TFI方法。

彈簧類比法是基于有限元思想的一種動網格方法，在網格節點之間通過“彈簧” 連接，使用迭代法求解節點的靜力平衡方程實現網格的變形，其基本方程就是如下的靜力平衡方程

式中，kij是邊i-j的彈簧剛度；kt,jj是為了避免單元過度變形而附加的扭轉剛度； Δxj是與有彈簧連接的第j個節點的位移； Δxi為待求的位移；N為與待求節點有彈簧連接的總節點數量。彈簧類比法的優點是可以適應結構網格和非結構網格，且其魯棒性良好。

超限插值方法（TFI）是一種適用于結構網格的動網格方法，通過采用代數插值將邊界上位移按距離邊界的遠近插值到內部網格。TFI方法的優勢是在小變形條件下能保證原有網格質量，因此非常適用于邊界層內計算，同時，代數插值的計算量小，速度很快；不足之處是僅能適用于結構網格，因此流體計算網格將采用結構網格劃分。

2 氣動熱彈性問題分析流程

動氣動彈性的分析，基本流程與靜氣動彈性的分析類似，不同的是，動氣動熱彈性分析增加了結構慣性力，關注點不再是結構的溫度和位移的平衡，而是結構在擾動下的動穩定性問題，通常采用時域方法進行結構動響應的計算，具體流程如圖1所示。

圖1 氣動熱彈性分析流程Fig.1 Aerothermoelastic analysis process

高超聲速舵面一般采用薄翼型設計，截面形狀主要有對稱菱形、梯形以及雙楔形等，如美國國家航空航天局（NASA）的典型高超聲速驗證飛行器X-43 的全動平尾截面為薄雙楔形，而其垂尾采用的是薄梯形。

本文研究舵面結構在經歷了長時間高超聲速飛行后的氣動熱彈性問題，建立高超聲速舵面模型平面尺寸如圖2所示。從圖2 可知，本文高超聲速舵面后掠角為32°，舵面底部中間處為固定面，如圖2（b）中A—A陰影部分所示，考慮到加工誤差和實際情況，舵面的前后緣均有厚度。舵面采用耐高溫鈦合金制造，考慮材料隨溫度的屬性變化，主要屬性變化如圖3所示。

圖2 舵面模型尺寸Fig.2 Model dimensions of rudder

圖3 材料屬性隨溫度變化曲線Fig.3 Curve of material properties with temperature

當前針對高超聲速氣動彈性中非定常氣動力的計算，基本會采用活塞理論，其他的一些方法（如非定常激波膨脹波理論和非定常牛頓碰撞理論）有時也被采用。這些方法都假設高超聲速氣流無黏性并且忽略真實氣體效應。盡管作了簡化，這些近似手段還是在一些特定情況下得到了足夠精確的結果[14]。考慮到本文模型和流場特性，也可以在保證足夠精度的情況下，忽略高超聲速氣流的黏性作用，采用工程算法計算非定常氣動力。

為了減小網格劃分對計算結果的影響，對流場網格進行無關性驗證，此算例的網格無關性驗證信息和評價結果分別見表1和表2。

表1 網格無關性驗證信息Table 1 Grid independence verification information

表2 網格收斂性評價Table 2 Grid convergence evaluation

采用舵面5°迎角時的升力系數作為驗證參數，不同網格密度下的升力系數計算結果如圖4所示。

圖4 不同網格密度下升力系數Fig.4 Lift coefficient at different mesh densities

從圖4 可知，不同網格密度下舵面升力系數差異隨著網格密度的增加趨于平穩，較密網格和極密網格之間的差異較小，為了確定需要采用的網格劃分形式，對不同網格密度的升力系數差異進行定量分析，見表2。

由表2可知，采用較密網格劃分最符合計算精度和所需計算資源平衡，因此流場計算網格采用較密網格方式劃分。

圖5 是舵面的流場計算網格，在產生激波和膨脹波位置均對網格進行了加密，以更加精確地計算流場特性。壁面第一層網格垂直高度根據網格無關性的驗證結果，采用5×10-4m，水平第一層網格寬度采用1×10-3m，網格垂直和水平增長率均為1.05，網格總數160萬個，來流數據采用高度為10km 處的標準大氣數據，其中靜壓為26499Pa，靜溫為223K，舵面迎角為0°，對流通量格式采用相對穩定的Roe格式，時間積分格式采用基于隱式的LU-SGS。

圖5 流場網格Fig.5 Flow field grid

使用Fluent/ANSYS 耦合的方法計算舵面的氣動彈性時，為了保證計算準確，采取緊耦合的方法，時間步長1×10-4s，選取翼梢上的一點作為監測點，根據其位移響應曲線情況，來判斷舵面是否發生顫振，監測點不同馬赫數下Z方向的位移響應如圖6所示。

圖6 舵面不同馬赫數下位移響應曲線Fig.6 Displacement response curve of rudder under different Mach numbers

從圖6（a）可知，監測點Z方向的位移響應是收斂的，但是收斂的速度緩慢，可見Ma5.75 下，已經接近舵面的臨界顫振速度。由圖6（b）可見，監測點Z方向的位移響應發散。綜上，此舵面常溫下的顫振馬赫數區間為5.75～6。

接著進行氣動熱作用下的氣動熱彈性分析，首先對舵面進行靜熱配平計算。熱配平計算時為了防止因加載的初始熱流以及氣動力過大，導致動網格計算出錯，所以在氣動熱流和氣動力加載時引入松弛系數0.1。同時，為了研究不同流動狀態、馬赫數和動壓對舵面氣動彈性的影響，根據先期計算，后續將分別計算層流和湍流條件下Ma4.5和Ma5下的舵面顫振特性，以及馬赫數相同時，層流和湍流同馬赫數不同動壓下的舵面顫振特性。下面將以湍流Ma4.5為例進行分析，圖7 是進行靜熱配平時翼梢上前、中、后三個監測點的溫度和位移變化曲線。

圖7 舵面監測點溫度和位移變化曲線Fig.7 Temperature and displacement curve of rudder surface monitoring point

從圖7 可知，由于在氣動熱流和氣動力加載時使用了松弛系數0.1，所以溫度有從低到高最后趨于穩定，Y負方向位移也有從小到大最后趨于穩定的過程，在迭代40步之后，溫度和位移值均趨于穩定，不再變化，認為此時靜熱配平完成。

圖8是Ma4.5湍流下熱配平后舵面的溫度分布云圖。從圖8可知，舵面溫度分布隨著舵面形狀變化，舵面前緣由于直接受到熱流沖擊，所以溫度最高，在舵面兩個轉折處產生的膨脹波等影響下，沿著X方向舵面溫度逐漸降低，后緣溫度最低。同時，由于迎角為0°，舵面是對稱雙楔形翼型，所以前緣表面溫度邊界與物面邊界有近似平行的分布。

圖8 熱配平后舵面溫度分布云圖Fig.8 Temperature distribution nephogram of rudder after hot trim

舵面在高超聲速氣流中，同時受到氣動力和氣動熱流的作用，舵面材料受熱會產生熱膨脹，同時氣動力也會使舵面因產生彈性變形而產生位移，圖9 是舵面熱配平后的總位移云圖，舵面位移從約束處沿著舵面展向逐漸增加，在翼梢處位移最大為8.4281mm。

圖9 熱配平后舵面位移云圖Fig.9 Displacement nephogram of rudder after hot trim

為了研究舵面的位移主要是由氣動加熱引起的，還是因為受到氣動力的作用，分別分析了只受到氣動熱和氣動力作用下的舵面位移，結果如圖10所示。

圖10 舵面不同載荷下的位移云圖Fig.10 Displacement nephogram of rudder under different loads

由圖10可知，舵面的熱位移主要是由于舵面受氣動熱而產生的，因氣動力產生的位移非常微小，最大位移只有0.1375mm，這也符合在迎角0°時，對稱舵面所受氣動力很小。

按照同樣的分析流程和方法分別計算層流和湍流時Ma4.5和Ma5下的舵面模態信息，將不同流動狀態下熱配平后計算得到前五階模態信息，表3中常溫為22℃，T表示湍流流動，L表示層流流動。

表3 不同流動下舵面前5階模態頻率Table 3 1～5 modal frequency of rudder under different flows

由表3可知，舵面受到氣動力和氣動加熱的作用，模態頻率有了明顯的下降，這是由于強烈的氣動加熱使材料的屬性發生改變。同時，在高超聲速流動中，湍流的傳熱系數遠大于層流，所以湍流狀態下的熱流密度遠高于層流，這很明顯會影響舵面的溫度分布，進而影響舵面模態頻率，表3中的各階頻率變化也體現了這一點，同馬赫數下，層流狀態下的頻率下降小于湍流下。相比于常溫狀態，Ma4.5 湍流下，1～5階的模態頻率分別下降了22.32%、16.99%、21.64%、21.07%和19.53% 。當湍流馬赫數從Ma4.5 提高到Ma5時，第一階模態頻率的下降幅度開始小于第二階(△f1=12.36%， △f2= 14.5%)，這就意味著隨著舵面溫度繼續升高，一階模態將會和二階模態相互靠近，根據經典顫振發生理論，彎—扭模態頻率越接近，將會越容易發生耦合，結構也就更容易顫振。

為了進一步分析舵面頻率的下降主要是由氣動加熱引起，還是由氣動力引起的，在Ma4.5湍流下，分別對只考慮氣動力作用，以及只考慮氣動熱作用的舵面進行分析，結果見表4。表4 中常溫為22℃，AF 表示只受氣動力作用，AH表示只受氣動熱作用。

表4 舵面不同載荷下前5階模態頻率Table 4 1～5 modal frequency of rudder under different loads

從表4可知，舵面模態頻率的下降由氣動加熱引起，當只受到氣動力作用時，模態頻率反而有所提高，這是由于舵面在氣動力作用下內部產生了應力，使舵面的剛度增加。可見，在高超聲速下，氣動熱流對結構的強烈加熱作用會極大影響飛行器結構特性，原因是隨著結構溫度的升高，舵面材料楊氏模量等參數反而降低，使得舵面的剛度降低，各階模態頻率下降。

在舵面的熱配平基礎上，按照動氣動熱彈性響應分析方法和流程，對舵面進行熱配平基礎上的結構響應分析。取翼梢上的一點作為監測點，耦合計算時間為1s，時間步長4×10-4s，通過監測點的位移響應是收斂還是發散，來判斷對應來流下舵面是否發生了顫振。圖11 是不同流動狀態和馬赫數下翼梢監測點Z方向的位移響應曲線。

圖11 層流和湍流下監測點的位移響應曲線Fig.11 Displacement response curves of monitoring points under laminar and turbulent flow

圖11中，在層流下，不管來流是Ma4.5還是Ma5，翼梢監測點的位移響應均是收斂的，說明在此來流條件下，舵面沒有發生顫振；而在Ma4.5湍流下，舵面熱結構的位移響應曲線呈發散趨勢，說明舵面在湍流Ma4.5下發生了顫振，同時從湍流Ma5的位移響應圖中，可以進一步確定舵面的顫振，此時發散的幅度較Ma4.5 時增大。為了進一步確定層流和湍流下各自的顫振速度區間，分別計算了層流Ma5.25和湍流Ma4.25 下舵面熱配平后翼梢同一監測點Z方向位移響應，結果如圖12所示。

圖12 不同流動下位移響應Fig.12 Displacement response under different flows

由圖12（a）可見，層流Ma5.25下，翼梢監測點位移響應呈發散趨勢，說明此時舵面已經發生顫振，結合圖11（b）可以確定，在層流條件下此舵面的顫振馬赫數區間為5～5.25。圖12（b）中，在Ma4.25 湍流下，翼梢監測點的位移響應呈收斂趨勢，同樣結合上文相關信息可以確定，在湍流條件下，此舵面的顫振馬赫數區間為4.25～4.5。

舵面常溫下的顫振馬赫數區間為5.75～6，因此，本文舵面模型在氣動熱作用下，顫振速度在層流狀態時下降了13.4%，而在湍流狀態時下降了26.9%。可見在高超聲速下，氣動加熱對舵面的顫振速度有顯著影響，因此在實際設計高超聲速舵面時，必須考慮氣動加熱作用對舵面結構特性和顫振速度的影響。

通過分析監測點位移響應數據可知，舵面在湍流Ma4.5 下的顫振頻率約為26Hz，此頻率介于一階模態頻率（19.523Hz）和二階模態頻率（58.619Hz）之間，這說明舵面顫振主要由一、二階模態耦合所引起，與經典顫振發生理論相吻合。

接著研究動壓對舵面氣動彈性的影響，下面以10km處標準大氣在湍流Ma4.5 下的動壓為基準（376kPa），計算舵面在相同馬赫數（Ma4.5）、不同動壓下同一監測點的位移響應變化。為了保證同馬赫數下動壓不同，顯然不同動壓下的空氣密度將會不同，密度減小會降低氣動加熱作用，為了減小其他因素的影響，在舵面熱配平時，仍采用10km 處標準大氣密度數據，只在舵面熱配平后進行動響應計算時，才采用不同的空氣密度數據。

選取層流和湍流作為不同流動條件，馬赫數均取5，計算以Ma4.5 湍流條件熱配平后，舵面模型在不同動壓下的位移響應，具體如圖13所示。

圖13 不同流動狀態和動壓下舵面位移響應Fig.13 Displacement response of rudder under different flow conditions and dynamic pressures

從圖13 可知，無論流動狀態是層流還是湍流，同馬赫數下動壓的降低均使得舵面的位移響應幅度減小，且本文在舵面熱配平時，密度是一致的，如果在熱配平階段同時考慮密度減小帶來的氣動加熱作用減弱，舵面的位移響應收斂趨勢將更加明顯，因此增大高超聲速飛行器的飛行高度，在同馬赫數下，不僅可以有效降低舵面氣動熱效應，同時也將延遲舵面顫振的出現。

3 結論

本文首先介紹了高超聲速下的氣動熱彈性分析方法和采用的湍流模型。通過對典型高超聲速翼型進行不同條件下的氣動熱彈性分析，得到的主要結論如下：

（1） 預測了舵面模型在常溫22℃下的顫振馬赫數區間為5.75～6； 對不同流動狀態和馬赫數下的舵面特性進行了計算，結果表明，不同流動狀態下對結構特性的影響有較大差異。同時，馬赫數的提高對結構顫振特性也有影響。

（2） 在湍流Ma4.5下考慮氣動熱的舵面顫振馬赫數區間，相比不考慮氣動熱時降低了26.9%，層流Ma4.5 時有13.4%的下降。

（3） 相同馬赫數下動壓的降低均會使舵面的位移響應幅度減小，因此，增大高超聲速飛行器的飛行高度，在相同馬赫數和飛行工況下，可以有效延遲舵面顫振的發生。