超聲速及高超聲速壁板顫振中的湍流邊界層效應

張 兵，韓景龍，錢 凱

（南京航空航天大學航空宇航學院，江蘇 南京210016）

引 言

壁板顫振是壁板結構在高速氣流中產生的一種自激振蕩現象，屬于典型的流固耦合問題。幾乎所有的高速飛行器的蒙皮都屬于壁板結構，發生壁板顫振會引發結構疲勞破壞，給飛行器安全帶來嚴重威脅。在過去的幾十年中，國內外學者已經就壁板顫振問題開展了大量研究［1～3］，但到目前為止，和試驗結果完全一致的分析結果還很少見［4］。其原因在于影響壁板顫振的因素較復雜，使得試驗條件難以準確模擬。

早在1963年，Fung就指出在壁板顫振的諸多影響因素中，湍流邊界層厚度效應是計算結果和試驗不一致的主要原因［5］。后來Muhlstein等專門進行了低超聲速（馬赫數1.0～1.4）條件下不同湍流邊界層厚度壁板顫振的試驗研究，結果表明隨著馬赫數增加，湍流邊界層厚度對顫振邊界影響逐漸增大，并在馬赫數為1.2時達到最大，而后隨馬赫數增加而迅速減小［6］。Dowell通過求解線性擾動方程來計算表面壓力場，并通過假定法向速度分布滿足1／7次方指數率來模擬邊界層，其計算結果和試驗結果也存在較大誤差［7］。近年來，隨著CFD技術的發展，CFD／CSD耦合分析方法逐漸成為解決流固耦合問題的主要手段。Friedmann和Zhong Xaolin等曾分別采用三階活塞理論、Euler以及N－S方程對給定運動的二維壁板進行對比計算，計算馬赫數為10.05［8］。結果發現Euler方程和三階活塞理論的氣動力十分吻合，但采用N－S的結果與二者有著本質的不同，并給出了高超條件下氣動力模型需要進一步考核的結論。隨后，Bendisken，Gordnier，Selvam先后采用了CFD方法對跨聲速、超聲速壁板顫振問題進行研究，對比了無粘和粘性結算結果，但缺少試驗驗證，且沒有針對邊界層效應進行討論［9～12］。最近，Atsushi等采用 CFD方法以及基于von Kármán板理論的有限元方法對Muhlstein等的試驗模型進行了計算［13］，得到了和試驗一致的結果，證明了采用 RANS（Reynolds－Averaged Navier－Stokes equations）方程可以獲得較準確的氣動力，但該文計算最大馬赫數為2.4，并未考慮高超聲速情形。因此，對于高超聲速壁板顫振問題，有必要采用能反映流場粘性特征的N－S方程來開展研究。

本文通過采用求解RANS方程實現非定常氣動力計算，并和結構有限元分析軟件耦合實現壁板顫振時域計算。建立和文獻［6］相同的壁板顫振模型，首先計算低超聲速下不同湍流邊界層厚度的顫振邊界，并和試驗結果進行比較驗證；然后將計算方法推廣到高超聲速，分析壁板顫振中的湍流邊界層效應，探討其作用機理。

1 計算方法

基于CFD／CSD的流固耦合計算分主要包括3部分：結構分析、氣動力計算、載荷及位移插值。下面介紹本文采用的方法。

1.1 基于有限元方法的結構動力學計算

以前研究較多采用von Kármán板理論對結構進行建模，考慮到目前結構動力學分析的有限元方法已經較為完善，以ANSYS，NASTRAN為代表的通用有限元分析軟件可以完成包含材料非線性、幾何大變形等非線性因素在內的結構分析。本文采用ANSYS軟件來實現結構計算。

采用有限元方法對壁板結構進行離散，其動力學方程為

式中x為節點自由度向量，M為質量矩陣，C為結構阻尼矩陣，K為結構剛度矩陣，F為載荷向量。當存在幾何非線性效應時，剛度矩陣K不再是常數矩陣，而是自由度x的函數，需要采用迭代法求解。在ANSYS中采用全瞬態分析，開啟大變形非線性效應，計算過程采用Newmark方法進行時間推進。

1.2 基于CFD方法的氣動力計算

采用迎風格式的有限體積方法是目前CFD計算的主流方法，本文采用自主開發的基于上述方法的三維N－S方程求解程序來實現氣動力計算。采用有限體積方法離散后的N－S方程為

式中Ω為控制體體積，W為控制體守恒變量，Fc和Fv分別為控制體第i個面上的對流通量和粘性通量，Q為控制體源項，Si為控制體第i個面的面積向量，nf為控制體面的個數。

為了精確捕捉激波，方程（2）中的對流項采用Roe格式，并通過MUSCL插值使空間離散精度達到二階。粘性通量采用二階中心差分格式計算。

時間離散方法采用二階Euler向后差分格式，如下式

式中 上標表示時間步序號，Δt為時間步長，R為右端殘值項，形式如下

采用Newton－Raphon方法求解非線性方程（3），通過將Wn＋1在Wn處線性化，并引入偽時間步得到隱式時間格式

式中 上標＊表示偽時間步。每一個偽時間步采用LUSGS方法求解。

由于流體方程是采用Euler描述法，而結構則采用Lagrangian描述法，因此對于包含網格運動的N－S方程其對流通量需要采用Arbitary Lagrangian Eulerian（ALE）公式計算

式中Vt為控制體面的運動速度。

此外為了減小時間離散誤差，控制體體積、面積矢量、面運動速度還需滿足幾何守恒率［14］。

為準確計算湍流邊界層，湍流模型選取Menter’s SST兩方程剪切應力模型［15］，考慮到高超聲速高雷諾數的特點，能量方程中需要包含湍流動能部分，此外湍流方程和質量、動量及能量方程全耦合求解。為避免駐點處湍流動能的累計誤差，湍流動能方程求解過程中的需要對產生項進行限制，本文取小于10倍的耗散項。

1.3 耦合計算平臺及計算流程

數據傳遞、載荷插值及迭代計算方法是耦合計算的3個關鍵問題。為了避免傳統I／O數據傳遞方式的效率瓶頸，本文在自主開發的基于共享內存的多場耦合計算平臺上進行多個軟件之間的數據交互［16］。

載荷插值精度直接決定了計算結果的精度，常用的載荷插值方法有樣條插值、形函數插值、守恒插值等。考慮到本文計算模型較簡單，結構網格節點和流體網格節點采用一對一形式，這樣可以避免由邊界載荷及位移差值帶來的計算誤差。

耦合計算流程依據耦合特點可分為松耦合、緊耦合兩種，松耦合方法各個物理場邊界條件采用另一側上個時間步的值，Gordnier曾指出這種兩個物理場間的不同步，在長物理時間計算中會引起較大累積誤差，建議采用緊耦合方式計算［10］。緊耦合方式需要在一個時間步對結構和流場進行內迭代求解來滿足時間步上的一致，其缺點是計算量較大，本文采用緊耦合方法計算。

2 計算模型

2.1 結構模型

計算結構模型如圖1所示，為周邊固支的矩形板。其中，a＝0.228 6m 為壁板沿流向尺寸，b＝2a為壁板展向尺寸，h＝0.001m為板厚度。結構材料為AZ31B－H24鎂合金，密度為1 762kg／m3，彈性模量為40.5GPa，泊松比為0.35。

圖1 壁板結構示意圖Fig.1 The panel structure model

ANSYS中采用SHELL63單元，結構沿流向和展向分別等距劃分20×40單元。ANSYS模態分析結果見表1，計算結果的第4，5階模態分別和試驗的第5，4階對應，這是由于試驗邊界并非完全固支的原因［6］。

表1 壁板結構模態的仿真結果Tab.1 Results of panel structure modal analysis

圖2 壁板結構模型前5階陣型圖Fig.2 First five mode shapes of panel structure

板內壓取流場初始計算結果作用力的平均值，為了激發壁板振動，在結構所有節點施加初始z方向3m／s的速度擾動。

2.2 流場網格及邊界條件

流場計算區域如圖3所示，流體網格一共分為12個結構網格塊，流場沿x，y，z三個方向的網格節點數為111×81×81，壁面第一層網格高度為1.0×10－6m，法向網格增長率為1.1。板所在平面（z＝0）內各個面均為絕熱壁面，為了防止入口處奇異，在距離入口0.1a距離的面設置為對稱面。其余各個邊界設置為超聲速遠場。

圖3 流場計算區域示意圖Fig.3 Computational domain of flow－field

邊界層厚度通過改變壁板前方到對稱面距離l來調整，其值取板中心位置處的厚度值，計算公式采用速度邊界層厚度公式

式中U為流體速度，U∞為來流速度。

3 計算結果

3.1 超聲速壁板顫振驗證

計算來流條件采用試驗值［6］，馬赫數范圍為1.05～1.4，時間步長取0.000 1s，最大子迭代步取6，子迭代步收斂殘差為0.01。結果中的顫振無量綱動壓由下式計算

式中υ為材料泊松比，Es為彈性模量，ρ∞為來流密度，U∞為來流速度。

3.1.1 Euler方程計算結果

圖4 Ma＝1.2時板中心點垂直位移響應時間歷程Fig.4 Vertical displacement response history of the center point at Ma＝1.2

如圖4為Ma＝1.2，λ＝105時的板中心點垂直方向無量綱位移（Δz／h）的時間歷程，壁板顫振形態為極限環振蕩，中心點的無量綱幅值為1.0，振動頻率為116.3Hz。如圖5為3個典型相位（中心點位值、零、最大值的3個相位的邊界層無量綱速度云圖，其中x值坐標范圍0～0.228 6為壁板，0.98的等值線為邊界層厚度。3個狀態下邊界層厚度基本不變化，但要超過來流厚度的10%，中心點的無量綱顫振幅值Δz／h在1左右，而邊界層厚度為24倍的板厚，因此振幅和邊界層厚度相比很小。

考慮板的振動位移，板中心點的邊界層厚度在一個周期內的變化規律相當于幅值為Δz的正弦運動，此時，邊界層起到了振動抑制的作用。

圖8所示為馬赫數1.2時無量綱顫振動壓隨邊界層厚度的變化曲線，為了便于比較，將Euler方程移分別為最小值、零和最大值）壁板表面的無量綱位移分布圖，由位移圖及頻率值可知，顫振主模態發生在一階，這與文獻［7］的結論吻合較好。

圖6為顫振邊界隨馬赫數的變化曲線，本文計算結果和文獻［7］結果一致，和試驗值吻合較好［6］。馬赫數為1.4時的顫振邊界和試驗值差別較大，是由于試驗結構邊界并不是完全固支，可以從結構模態計算結果看出。

圖5 Ma＝1.2，λ＝105時3個相位的壁板表面位移分布Fig.5 Panel displacement contours of three typical phases at Ma＝1.2，λ＝105

圖6 采用Euler方程的顫振動壓計算結果Fig.6 Flutter dynamic pressure results using Euler euqations

圖7 Ma＝1.2，λ＝180時3個典型相位的無量綱速度云圖Fig.7 Dimensionless velocity contours of three typical phases at Ma＝1.2，λ＝180

圖8 Ma＝1.2時顫振動壓隨邊界層厚度變化曲線Fig.8 Flutter dynamic pressure varying with boundary layer thickness at Ma＝1.2

3.1.2 N－S方程計算結果

通過求解RANS實現氣動力計算，計算了馬赫數為1.2時不同邊界層厚度下的壁板顫振動壓。圖7所示為λ＝180，δ／h＝0.11時中心點位移為最小的計算結果（δ＝0）也標記在圖中。其中空心圓點數據點為按照試驗結果外插得到的無粘結果。結果表明本文結果和試驗吻合較好，馬赫數為1.2時，邊界層厚度對顫振動壓影響非常大。當邊界層的無量綱厚度為0.1時，考慮邊界層效應的計算顫振動壓超過無粘結果的160%。

3.2 高超聲速壁板顫振

試驗中的來流馬赫數最大為1.4［6］，并未考慮高馬赫數下的情況，文獻［13］也僅計算到馬赫數為2.4。為考察高馬赫數下的邊界層效應，計算馬赫數選取了2.0，4.0，6.0，8.0四個工況，不考慮真實氣體效應，壁面按照絕熱壁面考慮，來流溫度為288 K。

如圖9結果為邊界層無量綱厚度δ／h＝0.087，Ma＝8.0，Δz／h＝1，λ＝4 806時的3個典型相位的無量速度云圖，由圖可知邊界層厚度基本保持不變。中心點位移為最小值、零、最大值的3個典型相位的位移如圖10所示，其運動形態上和圖5所示的低超聲速有較大不同。

圖9 Ma＝8，λ＝4 806時3個典型相位的無量綱速度云圖Fig.9 Dimensionless velocity contours of three typical phases at Ma＝8，λ＝4 806

圖10 Ma＝8，λ＝4 806時3個典型相位的無量綱位移云圖Fig.10 Dimensionless displacement contours of three typical phases at Ma＝8，λ＝4 806

圖11 不同馬赫數下N－S方程和Euler方程顫振動壓結果對比Fig.11 Comparison of flutter dynamic pressure at different Mach number using N－S and Euler equations

圖12 N－S方程和Euler方程結果差別隨馬赫數變化曲線Fig.12 Result difference varying with Mach number using Euler and N－S equations

如圖11為馬赫數為所有工況的結果比較圖，馬赫數超過2.0時，Euler方程和N－S方程的顫振動壓結果在絕對值上差別逐漸增大。圖12為二者結果的相對差別，馬赫數超過2后，相對差別基本保持在20%左右，且這個差別基本不隨馬赫數變化。高馬赫數下，采用N－S方程的計算顫振動壓高于Euler結果，此時邊界層主要起到了振動抑制的作用，在一定程度上提高了顫振穩定性。

4 結 論

（1）針對文獻［6］壁板顫振試驗模型，通過求解RANS方程實現氣動計算，和ANSYS軟件耦合實現壁板顫振的時域計算，計算結果和試驗值吻合較好，表明本文計算模型正確，計算方法精度較高。

（2）計算結果表明低超聲速階段湍流邊界層厚度對壁板顫振影響顯著，邊界層無量綱厚度為0.1時，其顫振動壓和Euler結果相比，差別最高可達160%，和試驗結果吻合較好。

（3）高超聲速情況下，邊界層厚度對顫振動壓仍有較大影響，其差別維持在20%左右，基本不隨馬赫數變化。湍流邊界層總體上起到了提高系統穩定性作用。因此準確預測壁板顫振問題，邊界層效應是必須考慮的因素。

［1］ Dowell E H.A review of the aeroelastic stability of plate and shells［J］.AIAA Journal，1970，8（3）：385—399.

［2］ Mei C Abdel－Motagaly K，Chen R.Review of nonlinear panel flutter at supersonic and hypersonic speeds［J］.Applied Mechanics Reviews，1999，52（10）：321—332.

［3］ 楊智春，夏巍.壁板顫振分析模型、數值求解方法和研究進展［J］.力學進展，2010，40（1）：81—98.Yang Zhichun，Xia Wei.Analytical models，numerical solutions and advances in the study of panel flutter［J］.Advances in Mechanics，2010，40（1）：81—98.

［4］ Dowell E H.Theoretical and experimental panel flutter studies in the Mach number range 1.0to 5.0［J］.AIAA Journal，1965，3（12）：2 292—2 304.

［5］ Fung Y C.Some recent contribution to panel flutter research［J］.AIAA Journal，1963，1（4）：898—909.

［6］ Muhlstein L，Gaspers P A，Riddle D W.An experimental study of the influence of the turbulent boundary layer on panel flutter［R］.NASA TND－4486，1968.

［7］ Dowell E H.Generalized aerodynamic forces on a flexible plate undergoing transient motion in a shear flow with an application to panel flutter［J］.AIAA Journal，1971，9（5）：834—841.

［8］ Nydick I，Friedmann P P，Zhong Xiaolin.Hypersonic panel flutter studies on curved panels［R］.AIAA－95－1458，1995.

［9］ Davis G A，Bendiksen O O.Transonic panel flutter［R］.AIAA－93－1476，1993.

［10］Gordnier R E，Visbal M R.Development of a three－dimensional viscous aeroelastic solver for nonlinear panel flutter［J］.Journal of Fluids and Structures，2002，16（4）：497—527.

［11］Gordnier R E，Visbal M R.High－fidelity computational simulation of nonlinear fluid－structure interaction problems structure interaction problems［J］.The Aeronautical Journal，2005，109（7）：301—331.

［12］Selvam R P，Visbal M R，Morton S A.Computation of nonlinear viscous panel flutter using a fully－implicit aeroelastic solver［R］.AIAA－98－1844，1998.

［13］Atsushi H，Takashi A.Effects of turbulent boundary layer on panel flutter［J］.AIAA Journal，2009，47（12）：2 785—2 791.

［14］Thomas P D，Lombard C K.Geometric conservation law and its application to flow computations on moving grids［J］. AIAA Journal，1979，17 （10）：1 030—1 037.

［15］ Menter F R.Two－equation eddy－viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32（8）：1 589—1 605.

［16］張兵，韓景龍.多場耦合計算平臺與高超聲速熱防護結構傳熱問題研究［J］.航空學報，2011，32（3）：400—409.Zhang Bing，Han Jinglong.Multi－field coupled computing platform and thermal transfer of hypersonic thermal protection structures［J］.Acta Aeronautica et Astronautica Sinica，2011，32（3）：400—409.