劉燦禮 袁麗蕓 王俊鵬 施國(guó)強(qiáng)

摘 要：針對(duì)平行四邊形薄板的自由振動(dòng)問(wèn)題，采用傳遞矩陣法，結(jié)合攝動(dòng)理論分析平行四邊形板，并通過(guò)與有限元法對(duì)比驗(yàn)證該方法的精確性.首先，利用坐標(biāo)變換公式將平行四邊形板的求解域轉(zhuǎn)換到矩形板域，并將控制方程及其對(duì)應(yīng)的邊界條件變換到矩形板域內(nèi)；隨后，結(jié)合Fourier級(jí)數(shù)展開(kāi)以及攝動(dòng)理論，將平形四邊行板的振動(dòng)控制方程寫(xiě)為一階常微分矩陣的形式；再運(yùn)用傳遞矩陣法對(duì)該一階常微分矩陣方程進(jìn)行求解，得到平形四邊行板的振動(dòng)固有頻率.算例結(jié)果表明，在平行四邊形板偏角較小的范圍內(nèi)，所提出的傳遞矩陣法可精確應(yīng)用于平行四邊形板的動(dòng)力學(xué)特性分析.

關(guān)鍵詞：平行四邊形板；傳遞矩陣法；攝動(dòng)法；一階常微分矩陣方程

中圖分類(lèi)號(hào)：U462.31 文獻(xiàn)標(biāo)志碼：A

0 引言

在彈性薄板的彎曲振動(dòng)問(wèn)題中，形狀規(guī)則的矩形板和圓形板較易獲得解析解.但對(duì)于平行四邊形板的振動(dòng)分析，因邊界條件不易滿足，其解析解極難獲得，因而常采用數(shù)值解法或近似函數(shù)解法進(jìn)行近似逼近.在數(shù)值解法中，有限元法邊界適用性最佳，能得到較精確的結(jié)果，但進(jìn)行參數(shù)化分析通常需要重復(fù)建模.在近似函數(shù)解法中，常用的有變分法.文獻(xiàn)[1]用康托洛維奇變分法對(duì)平行四邊形板進(jìn)行了分析，但由于在近似函數(shù)中采用了均布載荷下的梁函數(shù)，因而僅適用于均布載荷的情況，其他載荷情形需要重新求解，且求高階近似解的計(jì)算工作量也相當(dāng)大.楊柳等[2]采用GD方法對(duì)平行四邊形板的靜力學(xué)問(wèn)題進(jìn)行了研究，該方法數(shù)學(xué)原理嚴(yán)謹(jǐn)，精度較高，可拓展到動(dòng)力學(xué)問(wèn)題的求解，但其本質(zhì)仍然是用差分代替微分，存在一定的數(shù)值誤差.

上個(gè)世紀(jì)20年代以來(lái)，由于具有力學(xué)概念清晰、精度高、易于編程計(jì)算等優(yōu)點(diǎn)，傳遞矩陣法在結(jié)構(gòu)振動(dòng)中廣泛應(yīng)用[3] ，該方法的求解思路是將結(jié)構(gòu)控制微分方程寫(xiě)成一階常微分方程組的形式，則起點(diǎn)與終點(diǎn)的狀態(tài)向量間即可用傳遞矩陣建立簡(jiǎn)單關(guān)系，利用邊界條件對(duì)結(jié)構(gòu)振動(dòng)進(jìn)行求解[4].由于傳遞矩陣法中要計(jì)算矩陣運(yùn)動(dòng)微分方程中系數(shù)矩陣的指數(shù)函數(shù)的冪，導(dǎo)致其精度不高，為此，向宇[5]運(yùn)用微分方程與矩陣分析理論，結(jié)合精細(xì)積分技術(shù)，精確地計(jì)算了指數(shù)函數(shù)的冪，大大提高了傳遞矩陣法的計(jì)算精度.此后，傳遞矩陣法被廣泛應(yīng)用于復(fù)合多層梁[6]以及旋轉(zhuǎn)殼[7]等規(guī)則結(jié)構(gòu)的振動(dòng).該方法要求結(jié)構(gòu)振動(dòng)方程可降維寫(xiě)成一階常微分矩陣方程的形式，因此在形狀不規(guī)則的結(jié)構(gòu)中的應(yīng)用受到一定的限制.

基此，針對(duì)四邊簡(jiǎn)支平行四邊形板的自由振動(dòng)問(wèn)題，本文擬采用坐標(biāo)變換將平行四邊形板的求解域變換為矩形域，此后將借助攝動(dòng)理論，忽略方程中的小項(xiàng)，采用傅立葉級(jí)數(shù)展開(kāi)和無(wú)量綱化處理，得到平行四邊形板自由振動(dòng)的一階常微分矩陣方程，并擬借助高精度的傳遞矩陣法來(lái)求解.此方法為該類(lèi)結(jié)構(gòu)的參數(shù)化研究提供一個(gè)有力的手段.

1 平行四邊形板的動(dòng)力學(xué)問(wèn)題描述

考慮一厚度為h的各向同性均勻薄板，坐標(biāo)系位于中面（x，y）上.板中面x，y，z三個(gè)方向的位移分量分別記為us，vs，ws；u，v，w分別為板內(nèi)任意一點(diǎn)的3個(gè)位移分量.由薄板理論，可得其與中面位移分量的關(guān)系如下：

u=us-z■， v=vs-z■， w=ws （1）

與文獻(xiàn)[8]不同的是，式（1）中薄板中面內(nèi)的位移不為0，即薄板在彎曲振動(dòng)動(dòng)中，中面內(nèi)的各點(diǎn)位移有平行于中面的位移，且變形后的中面不再是平面.

結(jié)合薄板的幾何方程，以及廣義胡克定律，將其內(nèi)力向中面進(jìn)行簡(jiǎn)化可得內(nèi)力-位移關(guān)系：

（2）

其中，K=■，D=■ ；E，μ分別為板的彈性模量和泊松比，Nx，Ny分別為x，y方向的中面單位長(zhǎng)度的薄膜內(nèi)力，Nxy為中面單位長(zhǎng)度的薄膜剪力；Mx，My分別為x，y方向的中面單位長(zhǎng)度彎矩，Mxy為中面單位長(zhǎng)度扭矩.

由板殼振動(dòng)理論[9]，對(duì)薄板進(jìn)行受力分析可得平衡方程：

（3）

式（3）中，Qx，Qy分別為沿x，y方向的中面單位長(zhǎng)度橫向剪力；pi（i=x，y，z）為作用在單元體中的單位面積體力.

引入等效剪力：

Sx=Nxy，Vx=Qx+■ （4）

諧激勵(lì)下，式（1）～式（4）中的各物理量可展開(kāi)成幅值與諧激勵(lì)因子ejwt乘積的形式，即：

f（x，y，t）=■（x，y，w）ejwt （5）

則可得諧激勵(lì)下板的平衡方程如下：

（6）

同理可得， 內(nèi)力-位移關(guān)系（式（2））和等效剪力（式（4））在諧激勵(lì)下的形式，僅需將各物理量用其幅值替代即可.由于篇幅有限，不再一一列出.

為了便于使用邊界條件，需將平形四邊行板的求解域進(jìn)行坐標(biāo)變換成矩形板域.記α為平形四邊行板的偏角，采用如圖1所示坐標(biāo)變換，斜坐標(biāo)（x，y）與直角坐標(biāo)系（ξ，η）之間的轉(zhuǎn)換關(guān)系為：

x=ξ+ηsinα， y=ηcosα （7）

由鏈?zhǔn)角髮?dǎo)法則，有：

■=■ ■， ■=■ ■+■ ■ （8）

將式（8）代入諧激勵(lì)作用下的式（2）、式（4）、式（6）可得（ξ，η）坐標(biāo)系下平行四邊形板的內(nèi)力-位移關(guān)系如下：

（9）

運(yùn)動(dòng)控制方程如下：

（10）

等效剪力為：

■ξ η=■ξη ，■ξ=■ξ -tanα■+secα■ （11）

由于式（9）～式（11）均為偏微分方程，無(wú)法直接轉(zhuǎn)換為常微分方程，且在α較小時(shí)，tanα→0，可將上述方程中所含tanα項(xiàng)消除；因此，采用攝動(dòng)理論，取tanα為攝動(dòng)參數(shù)，將式（9）～式（11）中各物理量展開(kāi)為W=W0+W1tanα+W2tan2α+…的形式，并只保留第一項(xiàng)W0作為主項(xiàng)，誤差應(yīng)該不大.

采用矩形板類(lèi)似的求解方法，令ξ=Lc，（L=a），然后將各物理量沿η方向進(jìn)行傅里葉級(jí)數(shù)展開(kāi)并進(jìn)行無(wú)量綱化可得：

（■s，■s ）=L■（■s n，■s n）sin■，（■s ）=L■（■s n）cos■，（■ξ， ■η， ■ξ， ■ξ ）=K■ （■ξ n，■η n，■ξ n，■ξ n）sin■，

（■ξ η，■η ξ，■η ，■ξ ）=K■（■ξ ηn， ■η ξn， ■η n， ■ξ n）cos■， （■ξ ， ■η ）=（LK）■（■ξ n， ■η n）sin■，

（■ξ η， ■η ξ）=（LK）■（■ξ η n，■η ξn）cos■， （■ξ ， ■z ）=■（■ξ n， ■z n）sin■， （■η ）=■■η ncos■，

（■ξ ）=■■x ncos■ （12）

把式（12）代入只保留主項(xiàng)后的式（9）～式（11）可得平形四邊行板的內(nèi)力-位移關(guān)系：

（13）

平行四邊形板的平衡方程：

（14）

等效剪力：

■c n0=■c ηn0，■c n0=■c n0-■■cηn0 （15）

并引入無(wú)量綱轉(zhuǎn)角狀態(tài)變量的Fourier級(jí)數(shù)展開(kāi)分量：

■c n0=-■ （16）

從式（13）～式（16）中消去其他中間變量，經(jīng)推導(dǎo)整理后可得到平行四邊形板的一階常微分矩陣方程：

■=AZ0+F （17）

其中，Z0={■sn0，■sn0，■s n0，■c n0，■c n0，■c n0，■c n0，■c n0}T為板中面狀態(tài)變量的Fourier級(jí)數(shù)展開(kāi)分量，

A=

為系數(shù)矩陣，其非零元素如下：

A12=■，A15=1 ，A21=-■ ，A26=■ ，A34=-1 ，A43=-■ ，

A48=■ ，A51=-■ ，A56=■ ，A62=-■+■ ，

A65=-■，A73=-■+■ ，A78=■ ，A84=■ ，A87=1

F8×1=■[0，0，0，0，■c n0，■η n0，■z n0，0]T為所受外激勵(lì).

2 傳遞矩陣法的動(dòng)力學(xué)問(wèn)題求解

本文研究平行四邊形板的自由振動(dòng)特性不考慮外力，此時(shí)，式（17）中的F為零向量，由傳遞矩陣法，邊界兩端變量之間存在的傳遞關(guān)系：

Z1=TZ0 （19）

其中，Z0為起始狀態(tài)變量，Z1為終端狀態(tài)變量；T=eA為傳遞矩陣；引入文獻(xiàn)[10]描述的精細(xì)積分法，當(dāng)系數(shù)矩陣A已知時(shí)，可精確計(jì)算T的值.

通過(guò)坐標(biāo)變化后，四邊簡(jiǎn)支的平行四邊形薄板的求解域可轉(zhuǎn)換為一矩形域，其在矩形域兩邊的邊界條件為：

c=0，1+■sinα： us=0， Ncn0≠0， vsn0≠0， Scn0=0， ws=0， Vcn0≠0， Mcn0=0， θcn0≠0 （20）

結(jié)合式（19）和式（20），可得平形四邊行板自由振動(dòng)時(shí)的特征方程：

（21）

式（21）中，矩陣TT中的元素Tij表示傳遞矩陣T的第i行第j列元素.

要使特征方程（式（21））有非零解，則特征方程系數(shù)矩陣的行列式│TT│=det（TT）=0，這樣就可以得到一個(gè)關(guān)于板的固有圓頻率ω的方程，其求解過(guò)程可以利用MATALAB編程實(shí)現(xiàn).此后，可得平行四邊形薄板的振動(dòng)固有頻率如下：

f=■ （22）

3 算例分析

為了驗(yàn)證本文方法的正確性，下面將取一個(gè)偏角較小的平形四邊行薄板為例，進(jìn)行動(dòng)力學(xué)分析.

算例1 考慮尺寸為2.00 m×1.00 m×0.02 m，偏角α=π/36 的四邊簡(jiǎn)支平行四邊形薄板，其物理參數(shù)如表1所示.

通過(guò)本文方法和有限元仿真求得的（1，1）階振動(dòng)固有頻率的解如圖2和圖3所示.

其它階次固有頻率的本文解和有限元解如表2所示.從表2可以看出，隨著階數(shù)的提高，兩種方法所得固有頻率的差別逐漸增大，但本文提出的傳遞矩陣法求得的解和有限元解的誤差均低于2%，表明傳遞矩陣法具有較高的求解精度.誤差存在的主要原因可能在于忽略了偏角帶來(lái)的其他小項(xiàng)；因此，需要進(jìn)一步討論偏角的影響.

算例2 考慮平形四邊行板的 a， b，h 不變，但偏角發(fā)生變化時(shí)的情形（偏角分別取為α=π/36，α=π/18，α=π/9，α=π/6）；因此可得不同角度下本文解和有限元解，如表3所示.

算例2表明，隨著角度的增大，傳遞矩陣法的解和有限元解的誤差總體上是逐漸增大的，并且階數(shù)越高，誤差越大.當(dāng)角度小于π/9 （20°）時(shí)，兩者的誤差可控制在2%左右，為了節(jié)省篇幅，其余角度進(jìn)行的算例驗(yàn)證不再一一列表給出.通過(guò)不同角度下本文解和有限元解的對(duì)比表明：角度一旦超過(guò)π/9 （20°），誤差較大.

4 結(jié)語(yǔ)

采用傳遞矩陣法和攝動(dòng)理論，本文提出了一種分析平行四邊形板的動(dòng)力學(xué)特性的新方法.算例表明，在偏角較?。?～20°） 范圍內(nèi)，該方法和有限元法求得的固有頻率基本吻合.說(shuō)明該方法可應(yīng)用于平行四邊形薄板動(dòng)力特性問(wèn)題的求解，為其參數(shù)化研究提供了一個(gè)便利的手段.

參考文獻(xiàn)

[1] 王磊.平行四邊形彎曲板的康托洛維奇法[J]. 固體力學(xué)學(xué)報(bào)，1983，4（3）：420-426.

[2] 楊柳，彭建設(shè).解平行四邊形板彎曲問(wèn)題的GD法[J]. 成都大學(xué)學(xué)報(bào)（自然科學(xué)版），2014，33（3）：230-233.

[3] 朱鏡清.結(jié)構(gòu)振動(dòng)計(jì)算中傳遞矩陣法的應(yīng)用問(wèn)題[J]. 地震工程與工程振動(dòng)，1984，4（1）：36-44.

[4] 萬(wàn)浩川，李以農(nóng)，鄭玲.改進(jìn)的結(jié)構(gòu)振動(dòng)傳遞矩陣法[J]. 振動(dòng)與沖擊，2013，32（9）：173-177，82.

[5] 向 宇.分析結(jié)構(gòu)自由振動(dòng)的傳遞矩陣精確形式[J]. 振動(dòng)與沖擊，1999，18（2）：69-74.

[6] 吳寒劍，向宇，袁麗蕓.PCLD梁的振動(dòng)和阻尼特性分析[J]. 廣西科技大學(xué)學(xué)報(bào)，2014，25（4）：40-45，77.

[7] 張又軍，陸靜.環(huán)向加肋旋轉(zhuǎn)殼自由振動(dòng)特性分析[J]. 廣西科技大學(xué)學(xué)報(bào)，2014，25（1）：54-58.

[8] 徐芝綸.彈性力學(xué)[M]. 北京：高等教育出版社，2006.

[9] 曹志遠(yuǎn).板殼振動(dòng)理論[M]. 北京：中國(guó)鐵道出版社，1989.

[10] 向宇，黃玉盈，曾革委.精細(xì)時(shí)程積分法的誤差分析與精度設(shè)計(jì)[J].計(jì)算力學(xué)學(xué)報(bào)，2002，8（3）：276-280，319.

Abstract：The transfer matrix method is used to analyze the free vibration of parallel quadrilateral plate， and the accuracy of the method is verified by comparing with the finite element method. Firstly， using coordinate transformation formula， the solution domain with the parallel quadrilateral shape is converted into the rectangular one， and the governing equations and the corresponding boundary conditions in the transformed rectangular domain are obtained； then， combined with the Fourier series expansion and perturbation theory， the governing equation for the plate is written in the first order ordinary differential matrix form； next， the transfer matrix method is used to solve the matrix equations to get the natural frequency of the parallel quadrilateral plate. Example results show that when the parallelogram plate's angle is in a small range， the transfer matrix method proposed in this paper can be applied to the analysis of dynamic characteristics of a parallel quadrilateral plate.

Key words：parallel quadrilateral plate； transfer matrix method； perturbation method； the first order ordinary differential matrix equation

（學(xué)科編輯：張玉鳳）