仲繼澤， 徐自力

(西安交通大學(xué) 機械結(jié)構(gòu)強度與振動國家重點實驗室，西安 710049)

基于動網(wǎng)格降階算法的機翼顫振邊界預(yù)測

仲繼澤， 徐自力

(西安交通大學(xué) 機械結(jié)構(gòu)強度與振動國家重點實驗室，西安 710049)

基于彈性體動網(wǎng)格技術(shù)，發(fā)展了一種用于機翼流場網(wǎng)格變形的降階算法。將流場網(wǎng)格所包圍的空間區(qū)域視為虛擬彈性體。以虛擬彈性體變形的靜力平衡方程為基礎(chǔ)，結(jié)合機翼的振動控制方程，推導(dǎo)了機翼與虛擬彈性體的整體的振動控制方程。通過模態(tài)疊加方法計算機翼和流場網(wǎng)格節(jié)點的位移，進而得到變形后的流場網(wǎng)格。考慮到機翼顫振多為1階彎曲和扭轉(zhuǎn)振動，所以在流場網(wǎng)格節(jié)點位移的計算中只需考慮1階彎曲和扭轉(zhuǎn)振型。為了保證計算精度，在計算中同時考慮了2階彎曲和扭轉(zhuǎn)振型。RANS方程為流體控制方程，采用Spalart-Allamras湍流模型，結(jié)合動網(wǎng)格降階算法，對AGARD Wing 445.6顫振邊界進行了流固耦合計算。計算結(jié)果相對于實驗值的偏差小于2%，且與已有的彈性體動網(wǎng)格方法比，計算時間減少了54.8%。

機翼顫振；流固耦合；動網(wǎng)格；虛擬彈性體；模態(tài)疊加

機翼顫振是一種典型的流固耦合問題，它會使機翼產(chǎn)生損壞。美國空軍F-117戰(zhàn)機機翼因顫振發(fā)生破壞，導(dǎo)致該戰(zhàn)機墜毀[1]。美國洛克希德-馬丁公司的驗證機在飛行實驗中，機翼因發(fā)生顫振而折斷[2]。因此，機翼與空氣之間的流固耦合問題一直是航空工程領(lǐng)域研究的熱點。

傳統(tǒng)的機翼顫振預(yù)測的數(shù)值方法主要有V-g和p-k方法[3]。在傳統(tǒng)方法中，假設(shè)機翼做簡諧振動，并采用線性氣動模型分析機翼周圍流場。對于亞音速流動，采用偶極子格網(wǎng)法[4]分析機翼周圍的流動。對于超音速流動，采用piston理論[5]分析機翼周圍的流動。這些傳統(tǒng)方法的計算效率高，而且對亞音速及超音速顫振的預(yù)測精度較高。但是采用傳統(tǒng)方法預(yù)測機翼跨音速顫振時，需要通過縮比模型的風(fēng)洞試驗對所預(yù)測的顫振流速進行修正。而縮比模型的設(shè)計及其實驗數(shù)據(jù)的分析，通常需要1年時間[6]。這會嚴重拖慢飛行器的研發(fā)進度。鑒于上述原因，航空工程領(lǐng)域的研究人員引入了基于CFD(Computational Fluid Dynamics)/CSD(Computational Structural Dynamics)的流固耦合方法，以預(yù)測機翼顫振邊界[7-8]。流固耦合方法可以分為弱耦合和強耦合兩種算法[9]。弱耦合算法中，流場分析在時間域上先于振動分析一個時間步，容易累積誤差，引起發(fā)散。強耦合算法中，流場與機翼振動通過每一迭代步的氣動力及位移傳遞實現(xiàn)時間同步求解，具有更好的魯棒性。

機翼振動會改變流場的邊界，在流固耦合計算中為了考慮邊界變動對流場的影響，必須對流場網(wǎng)格節(jié)點坐標(biāo)進行更新。目前所發(fā)展的網(wǎng)格變形方法主要有彈簧法[10]、彈性體方法[11]、溫度體方法[12]及徑向基函數(shù)方法[13]。與其他動網(wǎng)格方法相比，彈性體方法是一種基于物理模型的網(wǎng)格變形方法，把流場網(wǎng)格變形作為一個彈性體變形問題進行求解，更能貼近結(jié)構(gòu)彈性變形的物理實際。另外，彈性體方法可以通過設(shè)置不均勻和各向異性彈性剛度等具體措施，保證變形過程中流場邊界層網(wǎng)格的質(zhì)量。因此，與其他的動網(wǎng)格方法相比，彈性體方法具有較強的網(wǎng)格變形能力，且變形后網(wǎng)格質(zhì)量較高[14]。但是彈性體方法比較復(fù)雜，計算量大。尤其是以強耦合的方式對機翼顫振進行計算時，每一時間步都需要對流場網(wǎng)格變形進行多次迭代求解。所以在流固強耦合計算中，采用彈性體方法更新流場網(wǎng)格會耗費大量計算資源。鑒于此，本文在現(xiàn)有彈性體動網(wǎng)格技術(shù)[15-17]的基礎(chǔ)上，發(fā)展了一種能夠快速計算機翼流場網(wǎng)格節(jié)點位移的動網(wǎng)格降階算法。采用頻率及振型描述機翼及流場網(wǎng)格變形。通過求解模態(tài)坐標(biāo)下解耦的動力學(xué)方程得到模態(tài)位移，然后采用模態(tài)疊加方法計算機翼流場網(wǎng)格節(jié)點的位移。考慮Wing 445.6的前4階模態(tài)，采用流固強耦合方法計算了其無因次顫振流速，計算結(jié)果與實驗值吻合很好。

1 控制方程

機翼附近的氣流一般為黏性可壓縮流動，可以采用RANS(Reynolds Averaged Navier Stokes)方程作為控制方程。

(1)

式中：Ω為有限控制體；?Ω為有限控制體的表面；n為控制體表面外法向；W為守恒性流體向量；Fc為無黏通量；Fv為黏性通量。

(2)

式中：ρ為流體密度；u為流速；ur為流體相對流場網(wǎng)格的流速；E為流體總能；p為流體壓力；τ為流體剪切應(yīng)力張量；q為熱流向量。

ur=u-ug

(3)

式中,ug為流場網(wǎng)格速度。

RANS方程不是封閉的，需要結(jié)合湍流模型來封閉方程。本文采用Spalart-Allamras模型[18]作為湍流模型。

機翼表面的氣流壓力，使機翼受到氣動力的作用。在氣動力作用下，機翼會產(chǎn)生振動。機翼振動的有限元控制方程為

(4)

式中：Mw為機翼總體質(zhì)量矩陣；Cw為機翼總體阻尼矩陣；Kw為機翼總體剛度矩陣；Xw為機翼結(jié)構(gòu)網(wǎng)格節(jié)點位移向量；Ffsi為機翼流固耦合面受到的氣動力。

在流固耦合面處，即機翼表面，機翼振動會改變流場邊界。因此，需要變更流場網(wǎng)格節(jié)點的位置。HUO等將流場網(wǎng)格所包圍的空間區(qū)域視為虛擬彈性體，通過求解虛擬彈性體變形的靜力平衡方程計算流場網(wǎng)格節(jié)點的位移。虛擬彈性體變形的靜力平衡方程為

KeXe=0

(5)

式中：Ke為流場網(wǎng)格域虛擬彈性體的總體剛度矩陣；Xe為流場網(wǎng)格節(jié)點位移向量。

2 流場動網(wǎng)格的降階算法

為了實現(xiàn)流場網(wǎng)格計算的降階，本文將流場網(wǎng)格域視為虛擬彈性體，并將機翼和虛擬彈性體作為一個整體系統(tǒng)。本小節(jié)將推導(dǎo)機翼與虛擬彈性體的整體振動控制方程，給出動網(wǎng)格降階算法的理論基礎(chǔ)。

機翼網(wǎng)格節(jié)點可以分為流固耦合面節(jié)點和非流固耦合面節(jié)點兩個部分，據(jù)此將機翼振動控制方程式(4)寫成分塊矩陣的形式。

(6)

式中：Xw_fsi為機翼網(wǎng)格流固耦合面節(jié)點位移向量；Xw_n_fsi為機翼網(wǎng)格非流固耦合節(jié)點位移向量。

流場網(wǎng)格節(jié)點也可以分為流固耦合面節(jié)點和非流固耦合面節(jié)點兩個部分，據(jù)此將流場網(wǎng)格域虛擬彈性體變形的靜力平衡方程式(5)寫成分塊矩陣的形式。

(7)

式中：Xe_fsi為流場網(wǎng)格流固耦合面節(jié)點位移向量；Xe_n_fsi為流場網(wǎng)格非流固耦合節(jié)點位移向量。

虛擬彈性體隨機翼振動發(fā)生變形，即在流固耦合面處，虛擬彈性體與機翼的位移具有連續(xù)性。可以采用有限單元位移插值算法[19]，對機翼結(jié)構(gòu)網(wǎng)格流固耦合面的節(jié)點位移Xw_fsi進行插值得到機翼流場流固耦合面節(jié)點的位移Xe_fsi。

Xe_fsi=NXw_fsi

(8)

式中,N插值系數(shù)矩陣。

因此，虛擬彈性體變形的靜力平衡方程可以改寫為

(9)

機翼和虛擬彈性體的整體的振動控制方程可以通過疊加機翼振動控制方程式(6)和虛擬彈性體變形的靜力平衡方程式(9)得到。

(10)

采用模態(tài)分析得到機翼與虛擬彈性體同步振動系統(tǒng)的模態(tài)頻率ω及質(zhì)量歸一化的正則振型Ψ，并通過坐標(biāo)變換得到正則坐標(biāo)下的解耦的動力學(xué)方程。

(11)

式中：I為單位矩陣；ξ為模態(tài)位移；ζ為模態(tài)阻尼比；Q為模態(tài)氣動力。

(12)

(13)

(14)

動力學(xué)方程式(11)的分量形式為

(15)

可以通過Wilson-θ方法求解上述的模態(tài)方程計算出模態(tài)位移ξi，然后采用模態(tài)疊加方法計算機翼及流場網(wǎng)格節(jié)點的位移。

(16)

式中,n為參與機翼流固耦合振動的模態(tài)數(shù)。

本文將機翼和虛擬彈性體作為一個整體系統(tǒng)，通過模態(tài)分析得到機翼和虛擬彈性體的整體的模態(tài)，然后采用模態(tài)疊加方法計算機翼和流場網(wǎng)格節(jié)點位移。考慮到機翼顫振多為一階彎曲或彎扭耦合振動[20]，所以在流場網(wǎng)格節(jié)點位移的計算中只需考慮一階彎曲或扭轉(zhuǎn)等低階模態(tài)，從而達到降階的目的。

3 流固耦合算法

本文采用強耦合算法對機翼顫振進行流固耦合分析。為加快收斂，采用穩(wěn)態(tài)流場的計算結(jié)果作為瞬態(tài)流場的初值。流固耦合分析的流程為

(1)采用SIMPLE算法[21]求解RANS方程，計算出流場的流速，壓力等變量。

(2)根據(jù)機翼流固耦合面上的流場壓力分布計算機翼的氣動力Ffsi，然后采用式(14)計算機翼的模態(tài)氣動力Qi。

(3)通過求解模態(tài)動力學(xué)方程式(15)計算模態(tài)位移ξi，然后采用模態(tài)疊加法計算機翼流場網(wǎng)格節(jié)點位移。

(4)根據(jù)流場網(wǎng)格節(jié)點位移更新流場網(wǎng)格節(jié)點坐標(biāo)，在新網(wǎng)格的基礎(chǔ)上再次求解RANS方程，計算流場流速，壓力等變量。

(5)觀察計算的收斂情況，如果計算不收斂，則重復(fù)(2)～(4)，直到計算收斂，如果計算收斂，則進行下一時間步的計算。

流場壓力、速度和總能的殘差小于某一特定值時，我們認為計算收斂。這一特定值的選取是與具體問題相關(guān)的。在本文的AGARD Wing 445.6算例中，當(dāng)壓力殘差<5.2×10-5，速度殘差<2.4×10-5，總能殘差<2.7×10-7時，我們認為計算達到收斂。

4 機翼和流場網(wǎng)格及網(wǎng)格變形

AGARD Wing 445.6[22]是由美國國家航天局蘭利研究中心(NASA Langley Research Center)設(shè)計的顫振實驗?zāi)Ｐ?，具有完整的風(fēng)洞實驗數(shù)據(jù)，是目前國際上用于機翼流固耦合算法驗證的標(biāo)準(zhǔn)算例。AGARD Wing 445.6采用NACA 65A004對稱翼型，展長0.762 m，根弦長0.558 7 m，1/4弦線后掠角為45°，展弦比1.65，根梢比0.66，見圖1。該機翼為木質(zhì)材料，其展向彈性模量為3.15 GPa，展向垂直方向的彈性模量為0.416 2 GPa，剪切模量為0.439 2 GPa，泊松比為0.31，密度為393.5 kg/m3。考慮機翼的材料阻尼，選取的材料阻尼比為0.03。機翼網(wǎng)格見圖2，共有8節(jié)點六面體單元5 974個，節(jié)點7 500個。

圖1 AGARD Wing 445.6 幾何尺寸Fig.1 Geometric configuration of AGARD Wing 445.6

圖2 AGARD Wing 445.6網(wǎng)格Fig.2 Mesh for AGARD Wing 445.6

邊界層網(wǎng)格的大小直接影響流場的在機翼表面的微觀流動特性，因此選取合適的邊界層網(wǎng)格至關(guān)重要。本文根據(jù)機翼的特征長度(弦長)、流場的雷諾數(shù)和機翼壁面的Y+值，采用NASA所提供的經(jīng)驗公式確定邊界層第一層網(wǎng)格的厚度。特征長度取翼尖弦長(0.368 2 m)。最大馬赫數(shù)(1.141 Ma)時的流場雷諾數(shù)為6.5×105。為本文采用Spalart-Allmaras湍流模型，該模型要求Y+<2，本文取Y+=1。據(jù)此本文計算得到邊界層第一層網(wǎng)格的厚度為1.0×10-5m。設(shè)定邊界層每一層的網(wǎng)格厚度以1.2的比例逐層增加。

圖3給出了采用不同數(shù)量的流場網(wǎng)格時計算得到的翼尖截面處上側(cè)的壓力分布?？梢钥闯觯S著流場網(wǎng)格的增多，計算結(jié)果逐漸趨于一致。流場網(wǎng)格數(shù)量為563 472時的計算結(jié)果與網(wǎng)格數(shù)量為649 872時的計算結(jié)果的相對偏差為0.02%。我們認為，當(dāng)流場網(wǎng)格數(shù)量為563 472時，增加網(wǎng)格數(shù)量基本不會提升計算精度，反而會增加計算規(guī)模。因此，本文選取網(wǎng)格數(shù)量為563 472的流場網(wǎng)格為本文的計算網(wǎng)格。機翼流場網(wǎng)格見圖4，共有節(jié)點588 080個。

圖3 Wing 445.6翼尖截面處上側(cè)壓力分布Fig.3 Pressure on the upside of the tip section of Wing 445.6

圖4 AGARD Wing 445.6 流場網(wǎng)格Fig.4 Flow mesh for AGARD Wing 445.6

當(dāng)設(shè)定虛擬彈性的為單一彈性模量時，流場網(wǎng)格變形過程中容易引發(fā)網(wǎng)格畸變(見圖5(a))。為了避免流場邊界層網(wǎng)格在變形過程中發(fā)生畸變，保證變形過程中流場網(wǎng)格的質(zhì)量，本文設(shè)定流場網(wǎng)格域虛擬彈性體的彈性模量與流場網(wǎng)格的大小相關(guān)，見式(17)。

(17)

在流固耦合面處，流場邊界層網(wǎng)格體積較小，其內(nèi)部的虛擬彈性體的彈性模量較大。離邊界層越遠網(wǎng)格體積就越大，網(wǎng)格內(nèi)部的彈性體彈性模量就越小。虛擬彈性體的彈性模量隨網(wǎng)格體積的增加以指數(shù)規(guī)律減小。當(dāng)流場網(wǎng)格隨機翼振動發(fā)生變形時，邊界層網(wǎng)格基本不發(fā)生彈性變形，而是隨著機翼振動作整體的剛體運動，從而避免流場網(wǎng)格變形過程中邊界層網(wǎng)格的畸變問題(見圖5(b))。

圖5 畸變和非畸變流場網(wǎng)格變形Fig.5 Flow mesh deformation with and without distortion

流場網(wǎng)格內(nèi)部的虛擬彈性體的彈性模量的最大值為50 Pa，遠小于機翼的彈性模量0.416 2 GPa。因此，與機翼的整體剛度矩陣相比，虛擬彈性體的總體剛度矩陣是小量。可以認為，雖然機翼在流固耦合面處與流場網(wǎng)格域的虛擬彈性體關(guān)聯(lián)，但是虛擬彈性體不會改變機翼的固有振動頻率及振型。采用本文的模態(tài)方法計算出的機翼的固有頻率(見表1)，本文的計算結(jié)果與實驗值吻合較好。機翼前4階模態(tài)振型對應(yīng)的流場網(wǎng)格變形(見圖6)，可以看出流場網(wǎng)格變形整體比較平順，沒有出現(xiàn)網(wǎng)格畸變。

表1 AGARD Wing 445.6頻率

圖6 AGARD Wing 445.6 不同模態(tài)下流場網(wǎng)格變形Fig.6 Flow mesh deformation against vibration mode for AGARD Wing 445.6

5 機翼顫振邊界預(yù)測

通常采用無因次流速表征機翼顫振邊界，表達式為

(18)

式中：vf為來流流速；b為翼根半弦長；ωα為機翼一階扭轉(zhuǎn)振動角頻率；μ為質(zhì)量比。

(19)

式中：m為機翼質(zhì)量；ρ為空氣密度；v(=0.13 m3)為參考體積。

在初始時刻，設(shè)置機翼第1～4階的模態(tài)初始位移為0，初始速度為0.023，0.033，0.074和0.054。在來流馬赫數(shù)為0.499的條件下，當(dāng)無因次流速為0.424時，機翼各階模態(tài)位移時間曲線(見圖7)。各階模態(tài)響應(yīng)的模態(tài)位移幅值隨著時間的推進逐漸減小，呈收斂趨勢。即，機翼不會發(fā)生顫振。當(dāng)無因次流速為0.468時，機翼各階模態(tài)位移時間曲線(見圖8)。各階模態(tài)響應(yīng)的模態(tài)位移幅值隨著時間的推進逐漸增加，呈發(fā)散趨勢。即機翼發(fā)生顫振。

圖7 無因次流速為0.424時的模態(tài)位移時間曲線Fig.7 Modal displacement vs time at reduced velocity 0.424

圖8 無因次流速為0.468時的模態(tài)位移時間曲線Fig.8 Modal displacement vs time at reduced velocity 0.468

機翼振動的穩(wěn)定性，可以采用對數(shù)衰減率(Logarithmic Decrement Ratio, LDR)進行評估。

(20)

式中：ln為自然對數(shù)符號；Ai+1為第個i+1周期的振動幅值；Ai為第i個周期的振動幅值。

在馬赫數(shù)為0.499的條件下，機翼振動的對數(shù)衰減率與無因次流速的關(guān)系曲線(見圖9)。在虛線的上方，LDR>0表示振動幅值隨時間增大，機翼發(fā)生顫振。在虛線的下方，LDR<0表示振動幅值隨時間減小，機翼不會發(fā)生顫振。曲線與虛線的交點處，LDR=0，表示振動幅值不隨時間變化，可以認為機翼處于顫振臨界點。在該顫振臨界點處，無因次流速為0.452，機翼各階模態(tài)位移時間曲線見圖10。本文設(shè)置的時間步長為4×10-4s，采用4核心2.9 GHz CPU 計算出圖10所示響應(yīng)的一個周期，耗時14.1 h。而采用HUO等所提出的彈性體方法，耗時為31.2 h。本文的方法可以減少計算時間54.8%。

圖9 馬赫數(shù)為0.499時對數(shù)衰減率-無因次流速曲線Fig.9 LDR vs. reduced velocity at the Mach number 0.499

圖10 無因次流速為0.452時的模態(tài)位移時間曲線Fig.10 Modal displacement vs. time at reduced velocity 0.452

本文計算出的Wing 445.6顫振邊界，見圖11?？梢钥闯?，同時考慮第1階和第2階的彎曲及扭轉(zhuǎn)模態(tài)的計算結(jié)果優(yōu)于只考慮第1階模態(tài)的計算結(jié)果。將考慮第1和第2階的計算結(jié)果與文獻的結(jié)果對比如下：在亞音速區(qū)域(Ma=0.499,0.678)，本文的結(jié)果與CHEN等研究的結(jié)果基本一致，且均與實驗值吻合較好。在跨音速區(qū)域(Ma=0.901，0.960)，本文的結(jié)果比CHEN等研究的結(jié)果更接近實驗值。本文采用Spalart-Allamras湍流模型，CHEN等采用Baldwin-Lomax湍流模型。由此可見，Spalart-Allamras模型比Baldwin-Lomax模型更好的模擬機翼的跨音速流動。這與文獻[24]的結(jié)論是一致的。本文和CHEN等采用RANS方程作為流體的控制方程，計算出的機翼的超音速(Ma=1.141)顫振邊界，與實驗值均有較大偏離。IM等[25]對這一現(xiàn)象進行了研究，發(fā)現(xiàn)RANS方程不能準(zhǔn)確捕捉超音速流動中激波邊界層干擾及其所引發(fā)的流動分離現(xiàn)象，導(dǎo)致計算出的超音速顫振邊界偏差較大。整體上，本文計算的顫振邊界相對于實驗值的偏差<2%。

圖11 AGARD Wing 445.6 顫振邊界Fig.11 Flutter boundary of the AGARD Wing 445.6

6 結(jié) 論

本文發(fā)展了一種動網(wǎng)格降階算法，可以實現(xiàn)流固耦合計算中機翼流場網(wǎng)格節(jié)點位移的快速計算。采用該方法對AGARD Wing 445.6顫振邊界進行了預(yù)測，得出如下結(jié)論:

(1)采用本文的動網(wǎng)格方法計算出的機翼顫振邊界與實驗測量得到的顫振邊界吻合。

(2)通過與CHEN等計算出的顫振邊界比較，發(fā)現(xiàn)對于跨音速顫振的預(yù)測，Spalart-Allamras湍流模型比Baldwin-Lomax湍流模型的準(zhǔn)確性好。

(3)采用RANS方程作為流體控制方程，對機翼超音速顫振進行預(yù)測會產(chǎn)生較大的偏差。

(4)與已有的彈性體方法相比，本文的方法可以減少計算時間54.8%。

[ 1 ] CRICKMORE P. Nighthawk F-117 stealth fighter[M]. Zenith Imprint, 2003.

[ 2 ] BURNETT E, ATKINSON C,BERANEK J, et al. NDOF Simulation model for flight control development with flight test correlation[C]//AIAA Modeling and Simulation Technologies Conference. Toronto：American Institute of Aeronautics and Astronautics, 2010: 7780-7794.

[ 3 ] SCHUSTER D M, LIU D D, HUTTSELL L J. Computational aeroelasticity: success, progress, challenge[J]. Journal of Aircraft, 2003, 40(5): 843-856.

[ 4 ] ALBANO E,RODDEN W P. A doublet-lattice method for calculating lift distributions on oscillating surfaces in subsonic flows[J]. AIAA Journal, 1969, 7(2): 279-285.

[ 5 ] ASHLEY H. Piston theory-a new aerodynamic tool for the aeroelastician[J]. Journal of the Aeronautical Sciences (Institute of the Aeronautical Sciences), 1956, 23(12): 1109-1118.

[ 6 ] YURKOVICH R N, LIU D D, CHEN P C. The state-of-the-art of unsteady aerodynamics for high performance aircraft[C]//39th Aerospace Sciences Meeting and Exhibit. Reno：AIAA, 2001.

[ 7 ] 謝亮, 徐敏, 李杰, 等. 基于 CFD/CSD 耦合的顫振與動載荷分析方法[J]. 振動與沖擊, 2012, 31(3): 106-110. XIE Liang, XU Min, LI Jie, et al. Flutter and dynamic analysis based on CFD/CSD coupling method [J]. Journal of Vibration and Shock, 2012, 31(3): 106-110.

[ 8 ] KEYE S. Fluid-structure coupled analysis of a transport aircraft and flight-test validation[J]. Journal of Aircraft, 2011, 48(2): 381-390.

[ 9 ] CHEN X,ZHA G C, YANG M T. Numerical simulation of 3D wing flutter with fully coupled fluid-structural interaction[J]. Computers & Fluids, 2007, 36(5): 856-867.

[10] 史愛明, 楊青, 楊永年. 非結(jié)構(gòu)運動網(wǎng)格下的三維機翼顫振數(shù)值分析[J]. 振動與沖擊, 2006, 24(6): 27-28. SHI Aiming, YANG Qing, YANG Yongnian. Numerical flutter analysis of a 3D wing using unstructured dynamic mesh Euler method [J]. Journal of Vibration and Shock, 2006, 24(6): 27-28.

[11] TEZDUYAR T E. Stabilized finite element formulations for incompressible flow computations[J]. Advances in Applied Mechanics, 1991, 28: 1-44.

[12] 陳炎, 曹樹良, 梁開洪, 等. 基于溫度體模型的動網(wǎng)格生成方法及在流固耦合振動中的應(yīng)用[J]. 振動與沖擊, 2010, 29(4): 1-5. CHEN Yan, CAO Shuliang, LIANG Kaihong, et al. A new dynamic grids based on temperature analogy and its application in vibration engineering with fluid-solid interaction [J]. Journal of Vibration and Shock, 2010, 29(4): 1-5.

[13] 謝亮, 徐敏, 張斌, 等. 基于徑向基函數(shù)的高效網(wǎng)格變形算法研究[J]. 振動與沖擊, 2013, 32(10): 141-145. XIE Liang, XU Min, ZHANG Bin, et al. Space points reduction in grid deforming method based on radial basis functions [J]. Journal of Vibration and Shock, 2013, 32(10): 141-145.

[14] BOTTASSO C L, DETOMI D, SERRA R. The ball-vertex method: a new simple spring analogy method for unstructured dynamic meshes [J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39): 4244-4264.

[15] BAR-YOSEPH P Z, MEREU S, CHIPPADA S, et al. Automatic monitoring of element shape quality in 2D and 3D computational mesh dynamics [J]. Computational Mechanics, 2001, 27(5): 378-395.

[16] STEIN K, TEZDUYAR T, BENNEY R. Mesh moving techniques for fluid-structure interactions with large displacements[J]. Journal of Applied Mechanics, 2003, 70(1): 58-63.

[17] HUO S H, WANG F S, YAN W Z, et al. Layered elastic solid method for the generation of unstructured dynamic mesh[J]. Finite Elements in Analysis and Design, 2010, 46(10): 949-955.

[18] SPALART P R, ALLMARAS S R. A one-equation turbulence model for aerodynamic flows[C]//30th Aerospace Sciences Meeting and Exhibit.Reno：AIAA,1992.

[19] ZIENKIEWICZ O C, TAYLOR R L, FOX D D. The finite element method for solid and structural mechanics[M]. Amsterdam:Elsevier, 2005：563-588.

[20] LIVNE E, WEISSHAAR T A. Aeroelasticity of nonconventional airplane configurations-past and future[J]. Journal of Aircraft, 2003, 40(6): 1047-1065.

[22] YATES J R,CARSON E. AGARD standard aeroelastic configurations for dynamic response I-wing 445.6[R]. France:Advisory Group for Aerospace Research and Development Neuilly-Sur-Seine,1988.

[23] BEAUBIEN R J, NITZSCHE F, FESZTY D. Time and frequency domain flutter solutions for the AGARD 445.6 wing[C]//International Forum on Aeroelasticity and Structural Dynamics.Munich：CAI, 2005.

[24] RUMSEY C L, SANETRIK M D, BIEDRON R T, et al. Efficiency and accuracy of time-accurate turbulent Navier-Stokes computations[J]. Computers & Fluids, 1996, 25(2): 217-236.

[25] IM H S, CHEN X Y, ZHA G C. Prediction of a supersonic wing flutter boundary using a high fidelity detached eddy simulation[C]//56th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference.Florida：AIAA, 2012: 9-12.

Wing flutter prediction using a reduced dynamic mesh method

ZHONGJize,XUZili

(State Key Laboratory for Strength and Vibration of Mechanical Structures,Xi’an Jiaotong University, Xi’an 710049, China)

In this paper, a reduced method for flow mesh deformation around a wing was developed based on the elastic solid method. The flow mesh domain was assumed to be a pseudo elastic solid. The total vibration equation for the wing with the pseudo elastic solid together was derived using the static equilibrium equation of the pseudo elastic solid and the vibration equation of the wing. The nodal displacements for the wing and flow mesh were computed through modal superposition and the deformed flow mesh was obtained. Considering that wing flutter often appeared as the 1stbending and torsion flutter, the nodal displacements for the flow mesh could be calculated by modal superposition of the 1stbending and torsion mode. To ensure the computational accuracy, the 2ndbending and torsion mode were also considered. The flutter boundary of the AGARD Wing 445.6 was predicted using the present dynamic mesh method coupled with the RANS equations and the Spalart-Allamras turbulent model. The relative error of the calculated results to the experimental data was less than 2%. The computing time was reduced by 54.8% compared with the pre-existing elastic solid method.

wing flutter; fluid structure interaction; dynamic mesh; pseudo elastic solid; modal superposition

國家自然基金(51275385);國家973(2011CB706505)

2015-10-19 修改稿收到日期：2016-01-16

仲繼澤 男，博士生，1988年7月生

徐自力 男，博士，教授，1967年3月生

V211.3; O327

A

10.13465/j.cnki.jvs.2017.04.029