張 兵,韓景龍,錢 凱
(南京航空航天大學(xué)航空宇航學(xué)院,江蘇 南京210016)
壁板顫振是壁板結(jié)構(gòu)在高速氣流中產(chǎn)生的一種自激振蕩現(xiàn)象,屬于典型的流固耦合問題。幾乎所有的高速飛行器的蒙皮都屬于壁板結(jié)構(gòu),發(fā)生壁板顫振會引發(fā)結(jié)構(gòu)疲勞破壞,給飛行器安全帶來嚴(yán)重威脅。在過去的幾十年中,國內(nèi)外學(xué)者已經(jīng)就壁板顫振問題開展了大量研究[1~3],但到目前為止,和試驗(yàn)結(jié)果完全一致的分析結(jié)果還很少見[4]。其原因在于影響壁板顫振的因素較復(fù)雜,使得試驗(yàn)條件難以準(zhǔn)確模擬。
早在1963年,F(xiàn)ung就指出在壁板顫振的諸多影響因素中,湍流邊界層厚度效應(yīng)是計(jì)算結(jié)果和試驗(yàn)不一致的主要原因[5]。后來Muhlstein等專門進(jìn)行了低超聲速(馬赫數(shù)1.0~1.4)條件下不同湍流邊界層厚度壁板顫振的試驗(yàn)研究,結(jié)果表明隨著馬赫數(shù)增加,湍流邊界層厚度對顫振邊界影響逐漸增大,并在馬赫數(shù)為1.2時(shí)達(dá)到最大,而后隨馬赫數(shù)增加而迅速減小[6]。Dowell通過求解線性擾動(dòng)方程來計(jì)算表面壓力場,并通過假定法向速度分布滿足1/7次方指數(shù)率來模擬邊界層,其計(jì)算結(jié)果和試驗(yàn)結(jié)果也存在較大誤差[7]。近年來,隨著CFD技術(shù)的發(fā)展,CFD/CSD耦合分析方法逐漸成為解決流固耦合問題的主要手段。Friedmann和Zhong Xaolin等曾分別采用三階活塞理論、Euler以及N-S方程對給定運(yùn)動(dòng)的二維壁板進(jìn)行對比計(jì)算,計(jì)算馬赫數(shù)為10.05[8]。結(jié)果發(fā)現(xiàn)Euler方程和三階活塞理論的氣動(dòng)力十分吻合,但采用N-S的結(jié)果與二者有著本質(zhì)的不同,并給出了高超條件下氣動(dòng)力模型需要進(jìn)一步考核的結(jié)論。隨后,Bendisken,Gordnier,Selvam先后采用了CFD方法對跨聲速、超聲速壁板顫振問題進(jìn)行研究,對比了無粘和粘性結(jié)算結(jié)果,但缺少試驗(yàn)驗(yàn)證,且沒有針對邊界層效應(yīng)進(jìn)行討論[9~12]。最近,Atsushi等采用 CFD方法以及基于von Kármán板理論的有限元方法對Muhlstein等的試驗(yàn)?zāi)P瓦M(jìn)行了計(jì)算[13],得到了和試驗(yàn)一致的結(jié)果,證明了采用 RANS(Reynolds-Averaged Navier-Stokes equations)方程可以獲得較準(zhǔn)確的氣動(dòng)力,但該文計(jì)算最大馬赫數(shù)為2.4,并未考慮高超聲速情形。因此,對于高超聲速壁板顫振問題,有必要采用能反映流場粘性特征的N-S方程來開展研究。
本文通過采用求解RANS方程實(shí)現(xiàn)非定常氣動(dòng)力計(jì)算,并和結(jié)構(gòu)有限元分析軟件耦合實(shí)現(xiàn)壁板顫振時(shí)域計(jì)算。建立和文獻(xiàn)[6]相同的壁板顫振模型,首先計(jì)算低超聲速下不同湍流邊界層厚度的顫振邊界,并和試驗(yàn)結(jié)果進(jìn)行比較驗(yàn)證;然后將計(jì)算方法推廣到高超聲速,分析壁板顫振中的湍流邊界層效應(yīng),探討其作用機(jī)理。
基于CFD/CSD的流固耦合計(jì)算分主要包括3部分:結(jié)構(gòu)分析、氣動(dòng)力計(jì)算、載荷及位移插值。下面介紹本文采用的方法。
以前研究較多采用von Kármán板理論對結(jié)構(gòu)進(jìn)行建模,考慮到目前結(jié)構(gòu)動(dòng)力學(xué)分析的有限元方法已經(jīng)較為完善,以ANSYS,NASTRAN為代表的通用有限元分析軟件可以完成包含材料非線性、幾何大變形等非線性因素在內(nèi)的結(jié)構(gòu)分析。本文采用ANSYS軟件來實(shí)現(xiàn)結(jié)構(gòu)計(jì)算。
采用有限元方法對壁板結(jié)構(gòu)進(jìn)行離散,其動(dòng)力學(xué)方程為
式中x為節(jié)點(diǎn)自由度向量,M為質(zhì)量矩陣,C為結(jié)構(gòu)阻尼矩陣,K為結(jié)構(gòu)剛度矩陣,F(xiàn)為載荷向量。當(dāng)存在幾何非線性效應(yīng)時(shí),剛度矩陣K不再是常數(shù)矩陣,而是自由度x的函數(shù),需要采用迭代法求解。在ANSYS中采用全瞬態(tài)分析,開啟大變形非線性效應(yīng),計(jì)算過程采用Newmark方法進(jìn)行時(shí)間推進(jìn)。
采用迎風(fēng)格式的有限體積方法是目前CFD計(jì)算的主流方法,本文采用自主開發(fā)的基于上述方法的三維N-S方程求解程序來實(shí)現(xiàn)氣動(dòng)力計(jì)算。采用有限體積方法離散后的N-S方程為
式中Ω為控制體體積,W為控制體守恒變量,F(xiàn)c和Fv分別為控制體第i個(gè)面上的對流通量和粘性通量,Q為控制體源項(xiàng),Si為控制體第i個(gè)面的面積向量,nf為控制體面的個(gè)數(shù)。
為了精確捕捉激波,方程(2)中的對流項(xiàng)采用Roe格式,并通過MUSCL插值使空間離散精度達(dá)到二階。粘性通量采用二階中心差分格式計(jì)算。
時(shí)間離散方法采用二階Euler向后差分格式,如下式
式中 上標(biāo)表示時(shí)間步序號,Δt為時(shí)間步長,R為右端殘值項(xiàng),形式如下
采用Newton-Raphon方法求解非線性方程(3),通過將Wn+1在Wn處線性化,并引入偽時(shí)間步得到隱式時(shí)間格式
式中 上標(biāo)*表示偽時(shí)間步。每一個(gè)偽時(shí)間步采用LUSGS方法求解。
由于流體方程是采用Euler描述法,而結(jié)構(gòu)則采用Lagrangian描述法,因此對于包含網(wǎng)格運(yùn)動(dòng)的N-S方程其對流通量需要采用Arbitary Lagrangian Eulerian(ALE)公式計(jì)算
式中Vt為控制體面的運(yùn)動(dòng)速度。
此外為了減小時(shí)間離散誤差,控制體體積、面積矢量、面運(yùn)動(dòng)速度還需滿足幾何守恒率[14]。
為準(zhǔn)確計(jì)算湍流邊界層,湍流模型選取Menter’s SST兩方程剪切應(yīng)力模型[15],考慮到高超聲速高雷諾數(shù)的特點(diǎn),能量方程中需要包含湍流動(dòng)能部分,此外湍流方程和質(zhì)量、動(dòng)量及能量方程全耦合求解。為避免駐點(diǎn)處湍流動(dòng)能的累計(jì)誤差,湍流動(dòng)能方程求解過程中的需要對產(chǎn)生項(xiàng)進(jìn)行限制,本文取小于10倍的耗散項(xiàng)。
數(shù)據(jù)傳遞、載荷插值及迭代計(jì)算方法是耦合計(jì)算的3個(gè)關(guān)鍵問題。為了避免傳統(tǒng)I/O數(shù)據(jù)傳遞方式的效率瓶頸,本文在自主開發(fā)的基于共享內(nèi)存的多場耦合計(jì)算平臺上進(jìn)行多個(gè)軟件之間的數(shù)據(jù)交互[16]。
載荷插值精度直接決定了計(jì)算結(jié)果的精度,常用的載荷插值方法有樣條插值、形函數(shù)插值、守恒插值等??紤]到本文計(jì)算模型較簡單,結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)和流體網(wǎng)格節(jié)點(diǎn)采用一對一形式,這樣可以避免由邊界載荷及位移差值帶來的計(jì)算誤差。
耦合計(jì)算流程依據(jù)耦合特點(diǎn)可分為松耦合、緊耦合兩種,松耦合方法各個(gè)物理場邊界條件采用另一側(cè)上個(gè)時(shí)間步的值,Gordnier曾指出這種兩個(gè)物理場間的不同步,在長物理時(shí)間計(jì)算中會引起較大累積誤差,建議采用緊耦合方式計(jì)算[10]。緊耦合方式需要在一個(gè)時(shí)間步對結(jié)構(gòu)和流場進(jìn)行內(nèi)迭代求解來滿足時(shí)間步上的一致,其缺點(diǎn)是計(jì)算量較大,本文采用緊耦合方法計(jì)算。
計(jì)算結(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,計(jì)算結(jié)果的第4,5階模態(tài)分別和試驗(yàn)的第5,4階對應(yīng),這是由于試驗(yàn)邊界并非完全固支的原因[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)壓取流場初始計(jì)算結(jié)果作用力的平均值,為了激發(fā)壁板振動(dòng),在結(jié)構(gòu)所有節(jié)點(diǎn)施加初始z方向3m/s的速度擾動(dòng)。
流場計(jì)算區(qū)域如圖3所示,流體網(wǎng)格一共分為12個(gè)結(jié)構(gòu)網(wǎng)格塊,流場沿x,y,z三個(gè)方向的網(wǎng)格節(jié)點(diǎn)數(shù)為111×81×81,壁面第一層網(wǎng)格高度為1.0×10-6m,法向網(wǎng)格增長率為1.1。板所在平面(z=0)內(nèi)各個(gè)面均為絕熱壁面,為了防止入口處奇異,在距離入口0.1a距離的面設(shè)置為對稱面。其余各個(gè)邊界設(shè)置為超聲速遠(yuǎn)場。
圖3 流場計(jì)算區(qū)域示意圖Fig.3 Computational domain of flow-field
邊界層厚度通過改變壁板前方到對稱面距離l來調(diào)整,其值取板中心位置處的厚度值,計(jì)算公式采用速度邊界層厚度公式
式中U為流體速度,U∞為來流速度。
計(jì)算來流條件采用試驗(yàn)值[6],馬赫數(shù)范圍為1.05~1.4,時(shí)間步長取0.000 1s,最大子迭代步取6,子迭代步收斂殘差為0.01。結(jié)果中的顫振無量綱動(dòng)壓由下式計(jì)算
式中υ為材料泊松比,Es為彈性模量,ρ∞為來流密度,U∞為來流速度。
3.1.1 Euler方程計(jì)算結(jié)果
圖4 Ma=1.2時(shí)板中心點(diǎn)垂直位移響應(yīng)時(shí)間歷程Fig.4 Vertical displacement response history of the center point at Ma=1.2
如圖4為Ma=1.2,λ=105時(shí)的板中心點(diǎn)垂直方向無量綱位移(Δz/h)的時(shí)間歷程,壁板顫振形態(tài)為極限環(huán)振蕩,中心點(diǎn)的無量綱幅值為1.0,振動(dòng)頻率為116.3Hz。如圖5為3個(gè)典型相位(中心點(diǎn)位值、零、最大值的3個(gè)相位的邊界層無量綱速度云圖,其中x值坐標(biāo)范圍0~0.228 6為壁板,0.98的等值線為邊界層厚度。3個(gè)狀態(tài)下邊界層厚度基本不變化,但要超過來流厚度的10%,中心點(diǎn)的無量綱顫振幅值Δz/h在1左右,而邊界層厚度為24倍的板厚,因此振幅和邊界層厚度相比很小。
考慮板的振動(dòng)位移,板中心點(diǎn)的邊界層厚度在一個(gè)周期內(nèi)的變化規(guī)律相當(dāng)于幅值為Δz的正弦運(yùn)動(dòng),此時(shí),邊界層起到了振動(dòng)抑制的作用。
圖8所示為馬赫數(shù)1.2時(shí)無量綱顫振動(dòng)壓隨邊界層厚度的變化曲線,為了便于比較,將Euler方程移分別為最小值、零和最大值)壁板表面的無量綱位移分布圖,由位移圖及頻率值可知,顫振主模態(tài)發(fā)生在一階,這與文獻(xiàn)[7]的結(jié)論吻合較好。
圖6為顫振邊界隨馬赫數(shù)的變化曲線,本文計(jì)算結(jié)果和文獻(xiàn)[7]結(jié)果一致,和試驗(yàn)值吻合較好[6]。馬赫數(shù)為1.4時(shí)的顫振邊界和試驗(yàn)值差別較大,是由于試驗(yàn)結(jié)構(gòu)邊界并不是完全固支,可以從結(jié)構(gòu)模態(tài)計(jì)算結(jié)果看出。
圖5 Ma=1.2,λ=105時(shí)3個(gè)相位的壁板表面位移分布Fig.5 Panel displacement contours of three typical phases at Ma=1.2,λ=105
圖6 采用Euler方程的顫振動(dòng)壓計(jì)算結(jié)果Fig.6 Flutter dynamic pressure results using Euler euqations
圖7 Ma=1.2,λ=180時(shí)3個(gè)典型相位的無量綱速度云圖Fig.7 Dimensionless velocity contours of three typical phases at Ma=1.2,λ=180
圖8 Ma=1.2時(shí)顫振動(dòng)壓隨邊界層厚度變化曲線Fig.8 Flutter dynamic pressure varying with boundary layer thickness at Ma=1.2
3.1.2 N-S方程計(jì)算結(jié)果
通過求解RANS實(shí)現(xiàn)氣動(dòng)力計(jì)算,計(jì)算了馬赫數(shù)為1.2時(shí)不同邊界層厚度下的壁板顫振動(dòng)壓。圖7所示為λ=180,δ/h=0.11時(shí)中心點(diǎn)位移為最小的計(jì)算結(jié)果(δ=0)也標(biāo)記在圖中。其中空心圓點(diǎn)數(shù)據(jù)點(diǎn)為按照試驗(yàn)結(jié)果外插得到的無粘結(jié)果。結(jié)果表明本文結(jié)果和試驗(yàn)吻合較好,馬赫數(shù)為1.2時(shí),邊界層厚度對顫振動(dòng)壓影響非常大。當(dāng)邊界層的無量綱厚度為0.1時(shí),考慮邊界層效應(yīng)的計(jì)算顫振動(dòng)壓超過無粘結(jié)果的160%。
試驗(yàn)中的來流馬赫數(shù)最大為1.4[6],并未考慮高馬赫數(shù)下的情況,文獻(xiàn)[13]也僅計(jì)算到馬赫數(shù)為2.4。為考察高馬赫數(shù)下的邊界層效應(yīng),計(jì)算馬赫數(shù)選取了2.0,4.0,6.0,8.0四個(gè)工況,不考慮真實(shí)氣體效應(yīng),壁面按照絕熱壁面考慮,來流溫度為288 K。
如圖9結(jié)果為邊界層無量綱厚度δ/h=0.087,Ma=8.0,Δz/h=1,λ=4 806時(shí)的3個(gè)典型相位的無量速度云圖,由圖可知邊界層厚度基本保持不變。中心點(diǎn)位移為最小值、零、最大值的3個(gè)典型相位的位移如圖10所示,其運(yùn)動(dòng)形態(tài)上和圖5所示的低超聲速有較大不同。
圖9 Ma=8,λ=4 806時(shí)3個(gè)典型相位的無量綱速度云圖Fig.9 Dimensionless velocity contours of three typical phases at Ma=8,λ=4 806
圖10 Ma=8,λ=4 806時(shí)3個(gè)典型相位的無量綱位移云圖Fig.10 Dimensionless displacement contours of three typical phases at Ma=8,λ=4 806
圖11 不同馬赫數(shù)下N-S方程和Euler方程顫振動(dòng)壓結(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時(shí),Euler方程和N-S方程的顫振動(dòng)壓結(jié)果在絕對值上差別逐漸增大。圖12為二者結(jié)果的相對差別,馬赫數(shù)超過2后,相對差別基本保持在20%左右,且這個(gè)差別基本不隨馬赫數(shù)變化。高馬赫數(shù)下,采用N-S方程的計(jì)算顫振動(dòng)壓高于Euler結(jié)果,此時(shí)邊界層主要起到了振動(dòng)抑制的作用,在一定程度上提高了顫振穩(wěn)定性。
(1)針對文獻(xiàn)[6]壁板顫振試驗(yàn)?zāi)P?,通過求解RANS方程實(shí)現(xiàn)氣動(dòng)計(jì)算,和ANSYS軟件耦合實(shí)現(xiàn)壁板顫振的時(shí)域計(jì)算,計(jì)算結(jié)果和試驗(yàn)值吻合較好,表明本文計(jì)算模型正確,計(jì)算方法精度較高。
(2)計(jì)算結(jié)果表明低超聲速階段湍流邊界層厚度對壁板顫振影響顯著,邊界層無量綱厚度為0.1時(shí),其顫振動(dòng)壓和Euler結(jié)果相比,差別最高可達(dá)160%,和試驗(yàn)結(jié)果吻合較好。
(3)高超聲速情況下,邊界層厚度對顫振動(dòng)壓仍有較大影響,其差別維持在20%左右,基本不隨馬赫數(shù)變化。湍流邊界層總體上起到了提高系統(tǒng)穩(wěn)定性作用。因此準(zhǔ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ìn)展[J].力學(xué)進(jìn)展,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]張兵,韓景龍.多場耦合計(jì)算平臺與高超聲速熱防護(hù)結(jié)構(gòu)傳熱問題研究[J].航空學(xué)報(bào),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.