超聲速及高超聲速壁板顫振中的湍流邊界層效應(yīng)

張 兵，韓景龍，錢 凱

（南京航空航天大學(xué)航空宇航學(xué)院，江蘇 南京210016）

引 言

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

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

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

1 計算方法

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

1.1 基于有限元方法的結(jié)構(gòu)動力學(xué)計算

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

采用有限元方法對壁板結(jié)構(gòu)進行離散，其動力學(xué)方程為

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

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

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

式中Ω為控制體體積，W為控制體守恒變量，F(xiàn)c和Fv分別為控制體第i個面上的對流通量和粘性通量，Q為控制體源項，Si為控制體第i個面的面積向量，nf為控制體面的個數(shù)。

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

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

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

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

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

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

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

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

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

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

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

載荷插值精度直接決定了計算結(jié)果的精度，常用的載荷插值方法有樣條插值、形函數(shù)插值、守恒插值等?？紤]到本文計算模型較簡單，結(jié)構(gòu)網(wǎng)格節(jié)點和流體網(wǎng)格節(jié)點采用一對一形式，這樣可以避免由邊界載荷及位移差值帶來的計算誤差。

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

2 計算模型

2.1 結(jié)構(gòu)模型

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

圖1 壁板結(jié)構(gòu)示意圖Fig.1 The panel structure model

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

表1 壁板結(jié)構(gòu)模態(tài)的仿真結(jié)果Tab.1 Results of panel structure modal analysis

圖2 壁板結(jié)構(gòu)模型前5階陣型圖Fig.2 First five mode shapes of panel structure

板內(nèi)壓取流場初始計算結(jié)果作用力的平均值，為了激發(fā)壁板振動，在結(jié)構(gòu)所有節(jié)點施加初始z方向3m／s的速度擾動。

2.2 流場網(wǎng)格及邊界條件

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

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

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

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

3 計算結(jié)果

3.1 超聲速壁板顫振驗證

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

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

3.1.1 Euler方程計算結(jié)果

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

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

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

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

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

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

圖6 采用Euler方程的顫振動壓計算結(jié)果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方程計算結(jié)果

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

3.2 高超聲速壁板顫振

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

如圖9結(jié)果為邊界層無量綱厚度δ／h＝0.087，Ma＝8.0，Δz／h＝1，λ＝4 806時的3個典型相位的無量速度云圖，由圖可知邊界層厚度基本保持不變。中心點位移為最小值、零、最大值的3個典型相位的位移如圖10所示，其運動形態(tài)上和圖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 不同馬赫數(shù)下N－S方程和Euler方程顫振動壓結(jié)果對比Fig.11 Comparison of flutter dynamic pressure at different Mach number using N－S and Euler equations

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

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

4 結(jié) 論

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

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

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

［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］ 楊智春，夏巍.壁板顫振分析模型、數(shù)值求解方法和研究進展［J］.力學(xué)進展，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，F(xiàn)riedmann 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］張兵，韓景龍.多場耦合計算平臺與高超聲速熱防護結(jié)構(gòu)傳熱問題研究［J］.航空學(xué)報，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.