基于二次特征值法的矩形薄板的動力穩(wěn)定性研究*

鐘子林，劉愛榮，盧漢文，黃友欽

(廣州大學(xué)-淡江大學(xué)工程結(jié)構(gòu)災(zāi)害與控制聯(lián)合研究中心，廣東 廣州 510006)

基于二次特征值法的矩形薄板的動力穩(wěn)定性研究*

鐘子林，劉愛榮，盧漢文，黃友欽

(廣州大學(xué)-淡江大學(xué)工程結(jié)構(gòu)災(zāi)害與控制聯(lián)合研究中心，廣東 廣州 510006)

基于Von-Karman薄板大撓度理論，利用Galerkin法得到面內(nèi)周期荷載作用下四邊簡支矩形薄板的二階常微Mathieu-Hill型參數(shù)振動方程；運用二次特征值法分別求出矩形薄板線性參數(shù)振動方程周期為2T和T時的主要與次要動力不穩(wěn)定域，并用有限元數(shù)值分析方法驗證了二次特征值法的精確性，同時定性地分析了主要參數(shù)共振下非線性彈性對系統(tǒng)定態(tài)振幅的影響。分析結(jié)果表明：① 當(dāng)激發(fā)力頻率近薄板兩倍自振頻率時，薄板發(fā)生強烈的橫向參數(shù)共振；② 二次特征值法可精確計算矩形薄板發(fā)生動力不穩(wěn)定時對應(yīng)的頻率和激發(fā)系數(shù)；③ 隨著薄板振幅的增長，非線性的存在抑制了定態(tài)振動幅值的無限增長，牽引系統(tǒng)向大頻率方向振動，導(dǎo)致振幅穩(wěn)定增加或迅速增大的復(fù)雜振動狀態(tài)。

矩形薄板；二次特征值法；參數(shù)共振；動力不穩(wěn)定域；周期荷載

板殼結(jié)構(gòu)在航空、建筑、橋梁、機械等結(jié)構(gòu)中有著廣泛的應(yīng)用，因此對板殼結(jié)構(gòu)力學(xué)行為的研究受到學(xué)者的普遍關(guān)注。結(jié)構(gòu)在壓力作用下處于平衡狀態(tài)，當(dāng)壓力值達到臨界荷載時，微小的外部干擾將會使其發(fā)生失穩(wěn)現(xiàn)象[1]。

例如，2004年，溫州大學(xué)體育場看臺的板殼結(jié)構(gòu)在臺風(fēng)“云娜”的襲擊下發(fā)生整體破壞的事故；2005年，寧波北侖體藝中心板殼結(jié)構(gòu)屋頂在臺風(fēng)“麥莎”的作用下出現(xiàn)撕毀破壞事故[2]。所以揭示板結(jié)構(gòu)在地震、風(fēng)荷或車輛等動力荷載作用下的穩(wěn)定性能亦是當(dāng)前研究的重要課題。

關(guān)于板結(jié)構(gòu)動力穩(wěn)定性能的研究已有廣泛涉及[3-13]。例如，鮑洛金[8]在其著作中探討了面內(nèi)受周期荷載作用的板在一定條件下會發(fā)生劇烈的參數(shù)共振問題；Konishi等[9]基于薄板小撓度理論，研究了不同邊界條件下面內(nèi)受周期變化荷載作用的矩形薄板的動力穩(wěn)定性；Ganapathi等[10]利用Newmark積分和迭代法研究了面內(nèi)受周期荷載作用的各向同性和各向異性彈性板的非線性動力不穩(wěn)定特性；Dey等[11]研究了面內(nèi)周期荷載作用下復(fù)合材料層合板的動力穩(wěn)定特性，并利用Bolotin法獲得不同設(shè)計參數(shù)薄板的主要和次要不穩(wěn)定域；Amabili[12]分別對不同邊界條件下受到面內(nèi)周期荷載的矩形薄板的非線性振動進行理論與實驗的研究分析；傅衣銘等[13]運用諧波增量法研究了四邊簡支矩形板在縱橫諧載共同作用下非線性動力穩(wěn)定性問題；彭凡[14]采用Boltzman積分型本構(gòu)關(guān)系，運用諧波增量法獲得粘彈性板在不同設(shè)計參數(shù)下的主要不穩(wěn)定域。

但目前關(guān)于板的動力穩(wěn)定性的研究主要是采用近似法來獲得動力系統(tǒng)的不穩(wěn)定域，當(dāng)動力系統(tǒng)的激發(fā)系數(shù)較大以及迭代次數(shù)較多時，此方法的計算結(jié)果誤差較大[14]；此外，許多文獻中并未考慮阻尼對矩形薄板動力穩(wěn)定性的影響，事實上，研究表明阻尼對結(jié)構(gòu)動力穩(wěn)定域的結(jié)果影響較大，尤其體現(xiàn)在次要不穩(wěn)定域的結(jié)果中[15]。本文針對前人研究的不足，采用二次特征值法精確求出矩形薄板線性參數(shù)振動方程周期為2T和T時的主要與次要不穩(wěn)定域，并用有限元數(shù)值法驗證了二次特征值法的精確性，同時定性地分析了主要參數(shù)共振下非線性彈性對系統(tǒng)定態(tài)振幅的影響。

1 動力系統(tǒng)模型

為探討薄板在周期荷載作用下的動力穩(wěn)定性能，本文選擇四邊簡支、長為a，寬為b，厚度為h的各向同性矩形薄板為研究對象(圖1)。并假設(shè)：矩形薄板的加載邊在面內(nèi)可以移動和轉(zhuǎn)動，非加載邊在面內(nèi)不可動；加載邊受周期簡諧荷載Px=Po+Ptcos(t)作用；板的四邊無剪切荷載，忽略板的面內(nèi)慣性力的作用。

基于Von-Karman薄板大撓度理論，考慮面外慣性作用以及阻尼的影響，得到板在外部激勵作用下系統(tǒng)的動力方程組為［16］:

(1)

式中，ω為板的橫向位移函數(shù)；φ為板應(yīng)力函數(shù)；ω為板的密度；ε為阻尼比，ε=ΩΔ/2π，Δ、Ω分別為承受縱向力不變分量的板的固有振動阻尼縮減率和固有振動頻率；D=Eh3/12(1-ν2)，E和ν為材料的彈性模量和泊松比。

圖1 動力系統(tǒng)模型Fig.1 Dynamic system model

板的非線性振動的邊界條件包括側(cè)向邊界條件和面內(nèi)邊界條件。

側(cè)向邊界條件(簡支)

(2)

面內(nèi)邊界條件

(3)

(4)

由薄板振動理論可知，板在振動時沿板長、寬方向形成屈曲半波，根據(jù)屈曲半波數(shù)選擇其橫向振動撓度函數(shù)［17］:

(5)

本文采用一階模態(tài)Galerkin逼近進行分析，選取橫向位移函數(shù)展開式中的第一項，即

(6)

式中，ω11和f11分別為基本模態(tài)的橫向位移和振動幅值，下標(biāo)表示板振動時沿板長、寬方向各形成一個屈曲半波，為簡單起見，在下文表述中省略其下標(biāo)。

將式(5)和式(6)代入方程組(1)的第二式，可得

(7)

根據(jù)應(yīng)力邊界條件(3)和(4)求得滿足式(6)的應(yīng)力函數(shù)

(8)

將應(yīng)力函數(shù)φ代入式(1)的第一式，并進行伽遼金積分，得

(9)

于是，方程(8)可表示成以下含非線性項的Mathieu-Hill方程

(10)

(11)

設(shè)f具有周期為2T的周期解展開的傅里葉級數(shù)形式為［8］:

(12)

同理，可設(shè)f具有周期為T的周期解展開的傅里葉級數(shù)形式為式(14)[8]。

為使方程組(13)、 (15)具有非零解，則其系數(shù)行列式應(yīng)等于零，從其系數(shù)行列式中可知周期為2T和周期為T時Mathieu-Hill 方程的激發(fā)系數(shù)μ和外荷載頻率與2倍結(jié)構(gòu)固有頻率比θ/2Ω之間存在一定的關(guān)系，通過進一步的分析可得到周期為2T和T時的動力不穩(wěn)定域[18]。

鮑洛金[7]應(yīng)用連分式的數(shù)值方法，對方程組(13)和(15)的系數(shù)行列式進行逐次迭代求得周期為2T和T時的動力不穩(wěn)定邊界。但是，這種方法在求解過程中需要較大的計算量，每次迭代只能得出一個不穩(wěn)定區(qū)域的邊界線，且在激發(fā)系數(shù)較大時，誤差較大。為較為準(zhǔn)確地得出振動系統(tǒng)的不穩(wěn)定域，本文采用二次特征值法求解周期為2T和T時振動系統(tǒng)的動力不穩(wěn)定邊界。

2 二次特征值法求解方程組

觀察周期為2T時振動系統(tǒng)的方程組，對其系數(shù)行列式進行分析。

(13)

(14)

(15)

(16)

(17)

(18)

方程組(14)可利用上面的矩陣表示為：

(19)

(20)

運用特征值法便可求得周期為2T時振動系統(tǒng)的不穩(wěn)定域邊界，同時觀察到周期為T時的系統(tǒng)方程不能簡單地按上述方法建立求解特征值的相應(yīng)方程。由于方程組(15)的系數(shù)矩陣右下部分B0的奇異性大，影響矩陣的特征值求解，可對該矩陣進行一次特征值計算，然后將求解得到的特征值重新加入到總矩陣中并運用上述求解周期為2T時動力不穩(wěn)定邊界的方法求解得振動系統(tǒng)周期為T時的動力不穩(wěn)定邊界。

對矩陣B0進行一次特征值求解，由于該矩陣對角元素含有實數(shù)1，求解時會出現(xiàn)無窮大的特征值，將無窮大的特征值舍去，得到含有特征值解的矩陣B：其中kn為特征值，將化簡后的矩陣B重新加入到總矩陣中，運用二次特征值法求解系數(shù)行列式(23)便可得到周期為T時的動力不穩(wěn)定邊界。

(21)

(22)

(23)

3 動力不穩(wěn)定性分析

基于以上理論分析，利用MATLAB編制程序，計算圖1所示薄板模型的動力不穩(wěn)定域。本文選取如下計算參數(shù)，即四邊簡支板薄板的長a=1 m，寬b=0.8 m，厚度h=5 mm;材質(zhì)為鋁材：彈性模量E=6.9×1010Pa，泊松比v=0.33，密度ρ=2 700 kg/m3；阻尼比ε=0.01。有限元模型見圖2。

圖3、4為矩形薄板一對邊在面內(nèi)受到周期荷載作用并考慮阻尼影響，周期解分別為2T(θ=2Ω/k，k=1，3，5，7，9)和T(θ=2Ω/k，k=2，4，6，8，10)時前五階動力不穩(wěn)定域(陰影部分)，其中橫坐標(biāo)表示外荷載作用頻率θ與2倍薄板固有頻率Ω比值的變化，縱坐標(biāo)表示外荷載pt與薄板屈曲荷載Pcr比值的變化，虛線是利用近似迭代法得到的主要不穩(wěn)定域解[7]。由圖3和圖4可知，面內(nèi)受周期荷載作用的矩形薄板除了會發(fā)生周期為T的強迫共振外，還會發(fā)生周期為2T的強烈橫向參數(shù)共振；此外，主要不穩(wěn)定域和次要不穩(wěn)定域的激發(fā)系數(shù)起始值大于零，這是由于考慮阻尼影響導(dǎo)致的，通常，若不考慮阻尼影響不穩(wěn)定域激發(fā)系數(shù)的起始值為零[15]；此外，從圖3和圖4中可以看出，在阻尼相同時，主要不穩(wěn)定域的起始值與激發(fā)系數(shù)原點的距離遠小于次要不穩(wěn)定與激發(fā)系數(shù)原點的距離，即阻尼對次要不穩(wěn)定域的影響較主要不穩(wěn)定域明顯；圖3和圖4亦表明隨著頻率比的降低，不穩(wěn)定域逐漸減小，當(dāng)阻尼一定時，較小的激發(fā)系數(shù)就可使薄板出現(xiàn)主要不穩(wěn)定域，而不會出現(xiàn)次要不穩(wěn)定域，即只有當(dāng)激發(fā)系數(shù)較大時才會出現(xiàn)次要不穩(wěn)定域；當(dāng)激發(fā)系數(shù)較大時，由二次特征值法求解得到的解與近似法求解結(jié)果差異較大，在求解過程中迭代法用先前計算得到的結(jié)果代入到求解下一個不穩(wěn)定區(qū)域的計算當(dāng)中，使得誤差逐級累加并且所需計算量較大，而二次特征值法的精度隨著所選取的諧波數(shù)的增加而增加。

為了驗證二次特征值法的正確性，利用有限元軟件ANSYS對該簡支薄板模型進行瞬態(tài)仿真模擬分析。根據(jù)圖1所示薄板模型的四邊簡支邊界條件，ANSYS建模時可約束四邊面外的自由度Uz，并分別約束平行于X軸的對稱軸面內(nèi)平動自由度Uy和平行于Y軸的對稱軸面內(nèi)平動自由度Ux。結(jié)合ANSYS在求解過程中的計算量與精確度，將薄板單元網(wǎng)格的大小設(shè)為25mm。

考慮到在次要不穩(wěn)定域中的薄板振動幅值較小，不能凸顯薄板在動力荷載作用下的不穩(wěn)定性，故在圖3、4中的主要不穩(wěn)定域中選取點A(0.65,0.6)、B(1,26,0.6)、C(0.36,0.6)、D(0.52,0.6)來驗證本文方法的準(zhǔn)確性。通過計算得到該矩形薄板的自振頻率為3.11Hz，以上所

選取點對應(yīng)的外荷載作用頻率分別為θA=4.04 Hz，θB=7.84 Hz，θC=2.24 Hz，θD=3.23 Hz，根據(jù)各點對應(yīng)的激發(fā)系數(shù)取周期外荷載幅值PtA=PtB=PtC=PtD=34 N，令外荷載恒值P0等于零，分別編寫相應(yīng)的ANSYS命令流，對受到簡諧荷載作用下的矩形薄板進行定頻模擬分析，得到點A、B、C、D的時域圖5-圖8。

圖5 點B的瞬態(tài)響應(yīng)Fig.5 Transient response of point B

圖6 點A的瞬態(tài)響應(yīng)Fig.6 Transient response of point A

圖7 點D的瞬態(tài)響應(yīng)Fig.7 Transient response of point D

圖8 點C的瞬態(tài)響應(yīng)Fig.8 Transient response of point C

由圖3和圖4可知：點B和點D位于本文方法計算得到的不穩(wěn)定域之內(nèi)而位于近似法求解得到的不穩(wěn)定域之外，結(jié)合此兩點的響應(yīng)圖可知，在與點B對應(yīng)的激發(fā)系數(shù)和頻率比相同的動力荷載作用下，薄板處于不穩(wěn)定的跳動狀態(tài)，當(dāng)該動力荷載持續(xù)作用于薄板時將會發(fā)生動力失穩(wěn)[8]，同樣，D點也是如此；同時觀察到點A和點C位于本文方法計算得到的不穩(wěn)定域之外而位于近似法求解得到的不穩(wěn)定域之內(nèi)，結(jié)合此兩點的響應(yīng)圖可知，在與點A對應(yīng)的激發(fā)系數(shù)和頻率比相同的動力作用下，薄板發(fā)生相對于板厚為微小振幅的穩(wěn)定狀態(tài)，沒有發(fā)生失穩(wěn)現(xiàn)象，C點亦然。綜上可知，二次特征值法比近似法更能準(zhǔn)確地反映動力系統(tǒng)發(fā)生失穩(wěn)時激發(fā)系數(shù)與外荷載作用頻率的關(guān)系。

(24)

式中，A=a2+b2，表示定態(tài)振動的振幅；Δ、Ω同式(1)，γ同式(9)。通過計算可得：λ=6.22×108；Ω=2π×3.11=19.54;μ=pt/2Pcr=0.6，取阻尼比為0.01，將相關(guān)系數(shù)代入式(24)中，可得含定態(tài)振幅與頻率比的非線性變化趨勢。

在周期荷載作用下，由線性系統(tǒng)所求的振幅是無限增長的，但這與實際的情況不相符合。由點B的響應(yīng)圖可知，動力系統(tǒng)發(fā)生參數(shù)共振時薄板的最大的位移為1.225mm，根據(jù)點B的定態(tài)振幅與頻率比的關(guān)系曲線(圖9)，可知點B位于穩(wěn)定解的上方；同理可知處于平穩(wěn)狀態(tài)的點A則位于穩(wěn)定解的下方，表明當(dāng)薄板發(fā)生參數(shù)共振時動力系統(tǒng)表現(xiàn)出了強烈的非線性特征，有效抑制了定態(tài)振動幅值的無限增長，牽引系統(tǒng)向大頻率方向振動，導(dǎo)致系統(tǒng)處于振幅穩(wěn)定增加或迅速增大的復(fù)雜振動狀態(tài)。

圖9 非線性系統(tǒng)下定態(tài)振幅與頻率比的關(guān)系Fig.9 The relation between steady amplitude and frequency ratio in nonlinear system

4 結(jié) 論

本文基于Von-Karman薄板大撓度理論，利用Galerkin法得到面內(nèi)周期荷載作用下四邊簡支矩形薄板的二階常微分Mathieu-Hill型參數(shù)振動方程，運用二次特征值法和ANSYS有限元數(shù)值方法分析了矩形薄板的動力不穩(wěn)定性，可以得出以下結(jié)論：

1)二次特征值法可精確計算矩形薄板發(fā)生動力失穩(wěn)時的主要不穩(wěn)定域和次要不穩(wěn)定域，較其他計算方法精度高。

2)當(dāng)外荷載作用的頻率與薄板的固有頻率的比值在兩倍附近時，面內(nèi)受到周期荷載作用的矩形薄板除了會發(fā)生強迫共振外，還可能引發(fā)強烈的橫向參數(shù)共振現(xiàn)象。

3)隨著頻率比的降低，矩形薄板的不穩(wěn)定域逐漸減?。蛔枘釋Υ我环€(wěn)定域的影響比較明顯，當(dāng)阻尼一定時，較小的激發(fā)系數(shù)就可使薄板出現(xiàn)主要不穩(wěn)定域，而不會出現(xiàn)次要不穩(wěn)定域，即只有當(dāng)激發(fā)系數(shù)較大時才會出現(xiàn)次要不穩(wěn)定域。

4)在薄板振動過程中，非線性因素隨著振動系統(tǒng)振幅的增長開始發(fā)生作用，其主要的作用是抑制定態(tài)振動幅值的無限增長，牽引系統(tǒng)向大頻率方向振動，逐漸遠離發(fā)生參數(shù)共振的頻率段，導(dǎo)致系統(tǒng)處于振幅逐漸減小或增大的復(fù)雜振動狀態(tài)。

[1] 禹奇才，劉愛榮，肖才濤，等. 基于尖點突變理論的平行組拼雙肋拱側(cè)傾失穩(wěn)臨界荷載計算新方法[J]. 中山大學(xué)學(xué)報 (自然科學(xué)版)，2014，53 (2)：6-11.

[2] 沈世釗, 武岳. 大跨空間結(jié)構(gòu)抗風(fēng)研究新進展與展望[C]∥ 第十三屆全國結(jié)構(gòu)風(fēng)工程學(xué)術(shù)會議.大連，2007.

[3]BIRMANV.Dynamicstabilityofunsymmetricallylaminatedrectangularplates[J].MechanicsResearchCommunications, 1985, 12(2):81-86.

[4]LEEHP,NGTY.Dynamicstabilityofamovingrectangularplatesubjecttoin-planeaccelerationandforceperturbations[J].AppliedAcoustics, 1995, 45(94):47-59.

[5]SINHASK.Commentson``dynamicstabilityofarectangularplatesubjectedtodistributedin-planedynamicforce''[J].JournalofSound&Vibration, 1989, 128(3):505-507.

[6]CHIACY.Largeamplitudevibrationsoflaminatedrectangularplates[J].FibreScience&Technology, 1982, 17(2):123-131.

[7]RAOGV,RAJUIS,RAJUKK,etal.Afiniteelementformulationforlargeamplitudeflexuralvibrationsofthinrectangularplates[J].Computers&Structures, 1976, 6(6):163-167.

[8] 鮑洛金. 彈性體系的動力穩(wěn)定性[M]. 北京：高等教育出版社, 1960.

[9]TAKAHASHIK,KONISHIY.Dynamicstabilityofarectangularplatesubjectedtodistributedin-planedynamicforce[J].JournalofSound&Vibration, 1988, 123(1):115-127.

[10]GANAPATHIM,PATELBP,BOISSEP,etal.Non-lineardynamicstabilitycharacteristicsofelasticplatessubjectedtoperiodicin-planeload[J].InternationalJournalofNon-LinearMechanics, 2000, 35(3):467-480.

[11]DEYP,SINGHAMK.Dynamicstabilityanalysisofcompositeskewplatessubjectedtoperiodicin-planeload[J].Thin-WalledStructures, 2006, 44(9):937-942.

[12]AMABILIM.Nonlinearvibrationsofrectangularplateswithdifferentboundaryconditions:theoryandexperiments[J].Computers&Structures, 2004, 82(s31/32):2587-2605.

[13] 傅衣銘, 張思進. 中厚矩形板的非線性動力穩(wěn)定性分析[J]. 工程力學(xué), 1995(3):97-106.

[14] 彭凡, 傅衣銘. 粘彈性板的非線性動力穩(wěn)定特性分析[J]. 固體力學(xué)學(xué)報, 2004, 25(1):115-118.

[15]HUANGYQ,LUHW,FUJY,etal.Dynamicstabilityofeulerbeamsunderaxialunsteadywindforce[J].MathematicalProblemsinEngineering, 2014, 2014(2):1-12.

[16] 袁尚平, 張建武, 王慶宇. 薄板3倍超諧振動的分析與試驗[J]. 上海交通大學(xué)學(xué)報, 2001, 35(7):955-961.

[17]AMABILIM.Theoryandexperimentsforlarge-amplitudevibrationsofrectangularplateswithgeometricimperfections[J].JournalofSound&Vibration, 2006, 291(3):539-565.

[18] 李康杰, 劉愛榮, 盧漢文,等. 圓弧淺拱面外動力穩(wěn)定性實驗研究[J]. 中山大學(xué)學(xué)報(自然科學(xué)版), 2015,54(3):30-35.

Dynamic stability research on rectangular thin plate based on the quadratic eigenvalue method

ZHONGZilin,LIUAirong,LUHanwen,HUANGYouqin

(Guangzhou University-Tamkang University Engineering Structure Disaster and Control Joint Research Center，Guangzhou University，Guangzhou 510006，China)

Based on the Von-Karman large deflection theory of thin plates, the second-order ordinary differential Mathieu-hill Parametric vibration equation was established for of simply supported rectangular thin plate under in-plane periodic loading using the Galerkin method. The quadratic eigenvalue method was adopted to obtain the main and secondary instability domains of the rectangular thin plate with the periodic solutions of 2T and T, then the accuracy of the method was verified by finite element method, meanwhile, the influence of nonlinear elasticity on stationary amplitude caused by the main Parametric resonance was qualitatively analyzed. The analytic results showed that strong transverse parametric resonance will occur when the frequency of excitation force is about twice as large as the natural vibration frequency of the thin plate. The quadratic eigenvalue method can be used to calculate accurately those parameters, such as the frequency and the excitation coefficient of the dynamic instability region of the rectangular thin plate. As the increasing of the amplitude, the infinite increase of vibration amplitude was controlled owing to the nonlinear factor which tows system to the direction of large vibration frequency. Furthermore, the nonlinear factor can bring the system into a complex vibration state in which the vibration amplitude increased stably or rapidly.

rectangular thin plate; quadratic eigenvalue method; parametric resonance; dynamic instability domain；periodic load

10.13471/j.cnki.acta.snus.2016.03.012

2015-12-30

國家自然科學(xué)基金資助項目(51578166) ;廣州市羊城學(xué)者資助項目(1201541551)；廣東省教育廳高層次人才資助項目

鐘子林( 1990年生)，男；研究方向：新型拱橋的靜、動力穩(wěn)定性；通訊作者： 劉愛榮;E-mail:liu-a-r@163.com

TU

A

0529-6579(2016)03-0068-09