超聲速氣流下黏彈性夾層壁板非線性顫振分析

呂海煒, 李映輝, 李中華, 李 亮

（1.西南交通大學(xué) 力學(xué)與工程學(xué)院，成都 610031；2.中國石油寶雞石油機(jī)械有限責(zé)任公司，寶雞 721002）

0 引言

物體由軸向運(yùn)動(dòng)誘發(fā)產(chǎn)生的橫向振動(dòng)[1-3]，對(duì)系統(tǒng)的振動(dòng)頻率和穩(wěn)定性有很大影響，而且可能導(dǎo)致運(yùn)動(dòng)失穩(wěn)。有研究表明，為了改善鴨式飛機(jī)（V-173）的橫向穩(wěn)定性，通常在其垂直尾翼下安裝阻尼夾層板[3]。這是一個(gè)典型高速飛行黏彈性阻尼夾層板的例子，其他飛行夾層結(jié)構(gòu)還有例如飛機(jī)蒙皮、航天器展開附件、高速飛行導(dǎo)彈外殼等。由軸向運(yùn)動(dòng)引起橫向振動(dòng)并導(dǎo)致失穩(wěn)的問題有廣泛背景，因此對(duì)超聲速氣流下黏彈性夾層壁板的顫振進(jìn)行研究具有重要的工程應(yīng)用價(jià)值。

壁板顫振問題研究始于20世紀(jì)50年代，學(xué)者們使用不同方法對(duì)該問題進(jìn)行了研究。Dowell[4]對(duì)1970年以前學(xué)者們?cè)诒诎孱澱穹治鲋兴龅墓ぷ鬟M(jìn)行了全面總結(jié)，主要包括使用的結(jié)構(gòu)理論和氣動(dòng)力理論，同時(shí)還詳細(xì)分析了上述各種理論的缺點(diǎn)并提出了相應(yīng)的改進(jìn)方法；Deman 和Dowell[5]基于二維不可壓縮渦格氣動(dòng)力理論從理論及實(shí)驗(yàn)上研究了兩端固支兩端自由的二維壁板的極限環(huán)振動(dòng)和線性顫振邊界；Jinsoo 和Younhyuck[6]使用非定常三維面法研究了機(jī)翼的超聲速顫振；Bogdan 等[7]基于von-Karmen 大變形理論和線性活塞理論，使用 有限差分法、Galerkin 法及正交分解法研究了非定常超聲速氣流下壁板的非線性振動(dòng)；Plaut 等[8]研究了各態(tài)歷經(jīng)靜態(tài)高斯隨機(jī)載荷的長(zhǎng)彈性板在超聲速氣流下的穩(wěn)定性；Bismarck-Nasr 和Bones[9]使用攝動(dòng)法得到固支平板的極限環(huán)振動(dòng)解，研究了結(jié)構(gòu)阻尼對(duì)非線性壁板顫振的遲滯影響；Beldica 等[10]作了大量關(guān)于黏彈性壁板顫振的研究，主要關(guān)注材料的性能和結(jié)構(gòu)的蠕變壽命時(shí)間；Kiiko[11]使用B-G 法和平均法研究了黏彈性矩形板的動(dòng)力穩(wěn)定性及黏性系數(shù)對(duì)臨界參數(shù)的影響；Gordnier 和Visbal[12]使用數(shù)值算法研究了非線性壁板顫振的三維黏彈性解；Kyo-Nam 和Woo-Seok[13]基于哈密頓原理及線性活塞理論得到了復(fù)合材料壁板的超聲速顫振方程，研究了遲滯和氣動(dòng)阻尼對(duì)顫振的影響；Bolotin 等[14]研究了初始條件對(duì)非線性彈性壁板的后臨界行為的影響，主要關(guān)注彈性壁板分叉邊界和混沌吸引子；張?jiān)品搴蛣⒄忌鶾15]研究了黏彈性材料壁板在超聲速來流作用下顫振時(shí)的分叉及混沌等復(fù)雜動(dòng)力學(xué)特性；肖艷平等[16]基于Kelvin-voigt 黏彈性本構(gòu)模型，根據(jù)von-Karmen 大應(yīng)變-位移關(guān)系和一階活塞氣動(dòng)力理論，建立了二維黏彈性壁板顫振方程，采用數(shù)值方法分析了黏彈性阻尼、面內(nèi)壓力及壁板幾何尺寸對(duì)黏彈性壁板顫振的影響，并研究了黏彈性壁板顫振時(shí)的分叉及混沌特性；李映輝 等[17]基于薄板小變形理論和一階氣動(dòng)力活塞理論建立了三維黏彈性夾層壁板的氣動(dòng)彈性顫振方程，并用Galerkin 截?cái)喾椒ㄑ芯苛讼到y(tǒng)的不同參數(shù)對(duì)顫振特性的影響。

由上可見，研究超聲速氣流下彈性結(jié)構(gòu)和黏彈性結(jié)構(gòu)顫振方面的文獻(xiàn)相當(dāng)多，但對(duì)黏彈性阻尼夾層結(jié)構(gòu)顫振研究較少。鑒于此，本文將研究超聲速氣流下黏彈性阻尼夾層結(jié)構(gòu)的顫振特性。

1 黏彈性夾層板顫振方程

圖1為超聲速氣流下黏彈性夾層板的模型，其幾何尺寸和材料參數(shù)與實(shí)際的軸向運(yùn)動(dòng)黏彈性夾層板一致，不同的是夾層板上表面有沿x方向的超聲速氣流時(shí)，夾層板將會(huì)受氣動(dòng)力作用。

圖1 黏彈性夾層壁板模型Fig.1 Model of the viscoelastic sandwich panel

在文獻(xiàn)[17]中若考慮幾何非線性，即黏彈性夾層板應(yīng)變?chǔ)舩、εy、γxy與撓度w的關(guān)系為：

且超聲速氣流作用下產(chǎn)生的氣動(dòng)載荷p采用三階活塞理論[18]可表示為

其中：qa=ρ∞V2/2 表示動(dòng)壓，ρ∞為來流空氣密度；M∞為馬赫數(shù)；V為來流速度；Cit和Cix(i=1, 2, 3)取0或1，就可分別討論氣動(dòng)力相應(yīng)子項(xiàng)對(duì)壁板氣動(dòng)彈性穩(wěn)定性的影響。

則可得超聲速氣流下黏彈性夾層壁板顫振的非線性振動(dòng)方程為

其系數(shù)見文后附錄。

2 Galerkin 離散

對(duì)四邊簡(jiǎn)支壁板，邊界條件為

設(shè)方程(3)的解形式為

其中：φmn(x,y)是滿足邊界條件的試函數(shù)；M、N為 截?cái)嚯A數(shù)。在四邊簡(jiǎn)支條件下通常取

將式(5)代入方程(3)，得系統(tǒng)的殘差R(x,y,t)（詳見文獻(xiàn)[19]）。使用Galerkin 方法，得

2.1 一階非線性分析

2.1.1 一階截?cái)喾匠?/p>

由文獻(xiàn)[18]可知，由于壁板的兩邊被簡(jiǎn)支，又與來流垂直，高階模態(tài)不易被激發(fā)，在y方向的截?cái)嗳∫浑A模態(tài)就可以體現(xiàn)此方向上變形的影響，故在垂直來流y方向上，N值取為1。

令M=1，得到一階截?cái)喾匠獭Ａ顇1=q11，x2=，得到二維一階常微分方程組：

其系數(shù)見附錄。

2.1.2 系統(tǒng)的平衡點(diǎn)及其穩(wěn)定性 系統(tǒng)的平衡點(diǎn)為

求得方程(8)有3 個(gè)平衡點(diǎn)，分別為：

要使平衡點(diǎn)X2和X3有意義，須滿足即來流速度應(yīng)滿足

當(dāng)來流速度超過此速度時(shí)，系統(tǒng)發(fā)生靜態(tài)分叉，平衡點(diǎn)的個(gè)數(shù)將由1 個(gè)變?yōu)? 個(gè)。

式(8)的Jacobi 矩陣為

X1=(0,0)處的特征方程為

其中λ1為X1的特征根，求得

λ1=令C22- 4C1= 0，得系統(tǒng)的臨界來流速度為

由C1＞0，C2＞0 可知該平衡點(diǎn)的兩個(gè)特征根均具有負(fù)實(shí)部，所以是穩(wěn)定的。

與討論平衡點(diǎn)X1的計(jì)算方法一樣，可計(jì)算平衡點(diǎn)X2、X3處Jacobi 矩陣的特征值，從而判斷其穩(wěn)定性。

2.1.3 數(shù)值模擬

超聲速氣流下黏彈性夾層壁板的幾何參數(shù)見表1，材料參數(shù)見表2。氣動(dòng)力參數(shù)為：馬赫數(shù)M∞=5，比熱比γ=1.4，C1x=C1t=C2x=C2t=C3x=C3t=1。

表1 黏彈性夾層板幾何參數(shù)Tabel 1 Geometrical size of viscoelastic sandwich panel

表2 黏彈性夾層板材料參數(shù)Tabel 2 Materials parameters of viscoelastic sandwich panel

圖2給出了平衡點(diǎn)隨來流速度變化的分叉圖。從圖2中可以看到，當(dāng)速度小于7571 m/s 時(shí)，系統(tǒng)只有一個(gè)平衡點(diǎn)X1；當(dāng)速度大于7571 m/s 時(shí)，系統(tǒng)的平衡點(diǎn)個(gè)數(shù)由1 個(gè)變?yōu)? 個(gè)。

圖2 平衡點(diǎn)隨來流速度變化Fig.2 The equilibrium points against the flow speed

圖3為平衡點(diǎn)X1的特征值λ1隨來流速度變化情況。從圖中看到，當(dāng)來流速度小于27 036 m/s 時(shí)，隨著來流速度的增加，特征值λ1的虛部絕對(duì)值逐漸減小，實(shí)部（始終為負(fù)）絕對(duì)值逐漸增大，故平衡點(diǎn)X1為穩(wěn)定的焦點(diǎn)；當(dāng)來流速度大于27 036 m/s時(shí)，特征值虛部為0，實(shí)部則分為兩支，但數(shù)值均小于0，說明此時(shí)平衡點(diǎn)仍然是穩(wěn)定的，只是平衡點(diǎn)類型由穩(wěn)定的焦點(diǎn)變?yōu)榉€(wěn)定的結(jié)點(diǎn)。

圖4為平衡點(diǎn)X2的特征值λ2隨來流速度變化情況。從圖中看到，特征值λ2虛部始終為0，且為一正一負(fù)，說明平衡點(diǎn)X2為不穩(wěn)定的鞍點(diǎn)。

圖4 特征值λ2 隨來流速度變化Fig.4 The eigenvalue (λ2) against the flow speed

圖5為特征值λ3隨來流速度變化情況。因特征值λ3實(shí)部與虛部絕對(duì)值相差較大，在圖5中不能完全反映λ3實(shí)部的真實(shí)走向，故單獨(dú)畫出其實(shí)部隨來流速度變化圖（圖6）。從圖5和圖6看到，特征值λ3虛部絕對(duì)值隨來流速度增大而增大，當(dāng)來流速度在區(qū)間(7571, 18 161) m/s 中，實(shí)部為負(fù)數(shù)，平衡點(diǎn)X3為穩(wěn)定的焦點(diǎn)；當(dāng)來流速度等于18 161 m/s 時(shí)，實(shí)部等于0，特征值為純虛數(shù)，此時(shí)平衡點(diǎn)X3為中心；當(dāng)來流速度超過18 161m/s，特征值實(shí)部變?yōu)檎龜?shù)，并隨來流速度的增大而增大，此時(shí)平衡點(diǎn)X3的類型為不穩(wěn)定焦點(diǎn)。

圖5 特征值λ3 隨來流速度變化Fig.5 The eigenvalue (λ3) against the flow speed

圖6 特征值λ3 實(shí)部隨來流速度變化Fig.6 The real part of eigenvalue (λ3) against the flow speed

表3給出了關(guān)鍵來流速度下平衡點(diǎn)的特征值及其平衡點(diǎn)類型。

表3 關(guān)鍵來流速度下平衡點(diǎn)的特征值及其類型Table 3 The eigenvalue at the equilibrium point and its type under typical flow speeds

2.2 二階非線性分析

2.2.1 二階截?cái)喾匠?/p>

令M=2，得二階截?cái)喾匠探M。令則有

其中相關(guān)系數(shù)C1，C2，．．．，C39見文后附錄。將式(15)寫為矢量形式為=f(x)，x∈R4。

2.2.2 系統(tǒng)發(fā)生Hopf 分叉的臨界來流速度

由f(x)=0，即可求得系統(tǒng)的平衡點(diǎn)。顯然X0= [0，0，0，0]是系統(tǒng)的平衡點(diǎn)。在X0處的Jacobi矩陣為

由分叉理論，考慮式(16)的特征值可得發(fā)生分叉的邊界曲線方程，則其特征方程為

其中：

由于馬赫數(shù)M∞、來流速度、黏彈性系數(shù)大于0，可知a1＞0，a2＞0，a3＞0，a4＞0。式(18)的Hurwitz行列式為

其中Δ1＞0，Δ2＞0，令Δ3=0，可解出臨界來流速度Vcr。

2.2.3 數(shù)值模擬

2.2.3.1 平衡點(diǎn)的穩(wěn)定性

非零平衡點(diǎn)的穩(wěn)定性也可通過上述方法判斷，因f(x)為高階非線性方程組，無法求其解析解，本文通過數(shù)值方法求其數(shù)值解。所使用的幾何參數(shù)見表1，材料參數(shù)見表2。

圖7給出了在平衡點(diǎn)X0處Δ3隨來流速度變化圖。數(shù)值分析結(jié)果表明，隨著來流速度的變化，系統(tǒng)始終只有平衡點(diǎn)X0。可以看到當(dāng)來流速度為1700 m/s時(shí)，Δ3=0。當(dāng)來流速度未達(dá)到臨界來流速度時(shí)，Δ3＞0，根據(jù)Hurwitz 定理判斷在此區(qū)間的平衡點(diǎn)X0是穩(wěn)定的焦點(diǎn)；當(dāng)來流速度為1700 m/s 時(shí)，Δ3=0，此時(shí)平衡點(diǎn)為穩(wěn)定極限環(huán)；超過臨界速度后Δ3＜0，平衡點(diǎn)變?yōu)椴环€(wěn)定焦點(diǎn)。表4給出了關(guān)鍵來流速度 下平衡點(diǎn)X0處Jacobi 矩陣的特征值。

圖7 Δ3 隨來流速度變化Fig.7 Δ3 against the flow speed

表4 關(guān)鍵來流速度下平衡點(diǎn)X0 特征值及類型Table 4 The eigenvalue at X0 and its stability under different flow speeds

2.2.3.2 Hopf 分叉圖

圖8給出了二階截?cái)嘞孪到y(tǒng)隨來流速度變化的Hopf 分叉圖，其中圖8(b)為來流速度在[1600, 2000]m/s 區(qū)域的分叉圖。

圖8 二階截?cái)嘞孪到y(tǒng)隨來流速度分叉圖Fig.8 Bifurcation diagram under the second order truncation

2.2.3.3 不同來流速度下系統(tǒng)的時(shí)間歷程圖和 相圖

當(dāng)V＜1700 m/s 時(shí)，系統(tǒng)僅有穩(wěn)定的平衡點(diǎn)X0，不論初值大小，最后都趨于此平衡點(diǎn)。圖9、圖10給出了二階截?cái)嘞聛砹魉俣萔＝1500 m/s、1690 m/s，初值取(0.000 01, 0.01, 0.000 01, 0.01)時(shí)系統(tǒng)的時(shí)間歷程圖和相圖?？梢钥吹较到y(tǒng)作衰減運(yùn)動(dòng)，說明系統(tǒng)是穩(wěn)定的。

圖9 二階截?cái)嘞聛砹魉俣萔=1500 m/s 時(shí)系統(tǒng)時(shí)間 歷程圖和相圖 （其中紅色*號(hào)表示初始點(diǎn)，綠色+號(hào)表示終結(jié)點(diǎn)，下同）Fig.9 The time history and phase portrait of the system when V=1500 m/s

圖10 二階截?cái)嘞聛砹魉俣萔=1690 m/s 時(shí)系統(tǒng)時(shí)間 歷程圖和相圖Fig.10 The time history and phase portrait of the system when V=1690 m/s

當(dāng)V≥1700 m/s時(shí)，平衡點(diǎn)X0變得不穩(wěn)定；但由于發(fā)生Hopf分叉產(chǎn)生穩(wěn)定的極限環(huán)，軌線均趨于穩(wěn)定的極限環(huán)。圖11、圖12分別給出了來流速度V=1700 m/s、2000 m/s，初值為(0.000 01, 0.01, 0.000 01, 0.01)時(shí)系統(tǒng)的時(shí)間歷程圖，可以看出，當(dāng)速度大于等于 1700 m/s時(shí)，系統(tǒng)作周期運(yùn)動(dòng)。

圖11 二階截?cái)嘞聛砹魉俣萔=1700 m/s 時(shí)系統(tǒng)時(shí)間歷程圖Fig.11 The time history of the system when V=1700 m/s

圖12 二階截?cái)嘞聛砹魉俣萔=2000 m/s 系統(tǒng)時(shí)間歷程圖Fig.12 The time history of the system when V=2000 m/s

3 結(jié)束語

本文基于von-Karmen 薄板大撓度理論、Kelvin-Voigt 黏彈性本構(gòu)方程和非線性活塞理論，建立了超聲速黏彈性夾層板顫振非線性方程，并使用Galerkin 方法研究該方程。對(duì)非線性系統(tǒng)主要研究了隨著來流速度的變化，系統(tǒng)的平衡點(diǎn)及穩(wěn)定性變化情況，得到如下結(jié)論：

1）對(duì)于一階截?cái)?，系統(tǒng)在來流速度小于臨界來流速度7571 m/s 時(shí)只有一個(gè)平衡點(diǎn)，且平衡點(diǎn)類型為穩(wěn)定焦點(diǎn)；當(dāng)來流速度超過臨界來流速度后，系統(tǒng)的平衡點(diǎn)個(gè)數(shù)，由1 個(gè)變?yōu)? 個(gè)。平衡點(diǎn)X2始終為不穩(wěn)定的鞍點(diǎn)；平衡點(diǎn)X3在來流速度 18 161 m/s 處穩(wěn)定性發(fā)生改變，由穩(wěn)定的焦點(diǎn)變?yōu)椴环€(wěn)定焦點(diǎn)。

2）對(duì)于二階截?cái)?，得到系統(tǒng)發(fā)生Hopf 分叉時(shí)的臨界速度為1700 m/s。使用數(shù)值方法研究了黏彈性夾層壁板隨來流速度變化的動(dòng)力學(xué)行為：來流速度在臨界速度前，系統(tǒng)作穩(wěn)定的衰減運(yùn)動(dòng)；臨界速度后，系統(tǒng)產(chǎn)生周期運(yùn)動(dòng)。

研究表明，對(duì)同一系統(tǒng)采用不同的截?cái)?，得到的臨界速度差別相當(dāng)大。應(yīng)該說，所取的截?cái)嚯A數(shù)越高，得到的結(jié)果越符合實(shí)際。

（References）

[1]Mote C D Jr.A study of band saw vibration[J].Journal of the Franklin Institute, 1965, 279(6)∶430-444

[2]Mote C D Jr.Dynamic stability of axially moving band[J].Journal of the Franklin Institue, 1968, 285(5)∶329-346

[3]Eun H B, Hwang M S, Park W J, et al.Lateral stability verification of a canard air plane with and without a vertical panel[J].Transactions of Nanjing University of Aeronautics and Astronautics, 2003, 20(2)∶125-133

[4]Dowell E H.Panel flutter∶a review of the aeroelastic stability of plates and shells[J].AIAA Journal, 1970, 8(3)∶385-399

[5]Deman T, Dowell E H.Limit cycle oscillations of two-dimensional panels in low subsonic flow[J].International Journal of Non-Linear Mechanics, 2002, 37(7)∶1199-1209

[6]Jinsoo C, Younhyuck C.Supersonic flutter analysis of wings using an unsteady 3D panel method[J].Computers and Fluids, 2001, 30(2)∶237-256

[7]Bogdan I E, Liaosha S T, Michael P P.Coherent structures and their influence on the dynamics of aeroelastic panel[J].International Journal of Non-Linear Mechanics, 2004, 39(6)∶977-991

[8]Plaut R H, Infante E F.On the stability of continuous systems subjected to random excitations[J].Journal of Applied Mechanics, 1970, 37(3)∶623-628

[9]Bismarck-Nasr M N, Bones C A.Damping effects in nonlinear panel flutter[J].AIAA Journal, 2000, 38(4)∶711-713

[10]Beldica C E, Hilton H H, Kubair D.Viscoelastic panel flutter-stability, probabilities of failure and survival times[C]//Proceedings of the 42ndAIAA/ASME/ASCE/ AHS/ASC Structures, Structural Dynamics and Materials Conference.Seattle, WA, USA, 2001∶3423-3433

[11]Kiiko A.Flutter of a viscoelastic plate[J].Journal of Application Mathematics and Mechanics, 1996, 60(1)∶167-170

[12]Gordnier R E, Visbal M R.Development of three- dimensional viscous aeroelastic solver nonlinear panel flutter[J].Journal of Fluids and Structures, 2002, 16(4)∶497-527

[13]Kyo-Nam K, Woo-Seok H．Effects of hysteretic and aerodynamic damping on supersonic panel flutter of composite plates[J]．Journal of Sound and Vibration, 2004, 273(3)∶569-583

[14]Bolotin V V, Grishko A A, Kounadis A N, et al.Influence of initial conditions on the post-critical behavior of a nonlinear aeroelastic system[J].Nonlinear Dynamics, 1998, 15(1)∶63-81

[15]張?jiān)品? 劉占生.粘彈壁板顫振的非線性動(dòng)力特性[J].推進(jìn)技術(shù), 2007, 28(1)∶103-107 Zhang Yunfeng, Liu Zhansheng.Nonlinear dynamic analysis of viscoelastic panel flutter[J].Journal of Propulsion Technology, 2007, 28(1)∶103-107

[16]肖艷平, 楊翊仁, 葉獻(xiàn)輝.三維粘彈壁板顫振分析[J].振動(dòng)與沖擊, 2011, 30(1)∶82-86 Xiao Yanping, Yang Yiren,Ye Xianhui.Flutter analysis for a three-dimensional viscoelastic panel[J].Journal of Vibration and Shock, 2011, 30(1)∶82-86

[17]李映輝, 李中華.超聲速下粘彈性夾層壁板顫振分析[J].力學(xué)季刊, 2012, 33(3)∶449-455 Li Yinghui, Li Zhonghua.Flutteranalysis of viscoelastic sandwich panel under supersonic speed[J].Chinese Quarterly of Mechanics, 2012, 33(3)∶449-455

[18]葉獻(xiàn)輝.壁板非線性氣動(dòng)彈性顫振及穩(wěn)定性研究[D].成都∶西南交通大學(xué)博士學(xué)位論文, 2008

[19]李中華.軸向運(yùn)動(dòng)黏彈性夾層板振動(dòng)分析[D].成都∶西南交通大學(xué)碩士學(xué)位論文, 2012

附 錄

且