基于小波Galerkin法的矩形薄板二次屈曲分析*

張 磊，張文明，王 林，李世斌

（1.國防科技大學 空天科學學院，長沙 410073；2.湖北三江航天集團 紅陽機電有限公司，湖北 孝感 432000）

引言

自20世紀90年代小波分析理論建立以來，小波數值方法因其在非線性和多尺度計算方面的應用潛力，得到了持續不斷的關注和研究，尤其是在非線性力學問題求解方面的應用[1]，例如結構大撓度分析[2]、智能結構振動控制、應力集中和裂紋擴展等應力奇異性問題[3].為了繼續發掘小波方法在非線性力學問題求解方面的計算優勢，本文在矩形薄板大撓度問題小波解[2]的基礎之上，繼續采用小波方法對矩形薄板非線性穩定性問題求解與分析.

對于結構完善的彈性平板，當面內壓縮載荷超過Euler 屈曲臨界值，它們就由平直的初始狀態進入橫向變形的后屈曲狀態，并以屈曲波形的形式保持穩定[4].隨著載荷的增加，屈曲波形保持不變而其幅值呈非線性增大.大量理論研究和穩定性實驗還表明，矩形薄板在后屈曲變形增大的過程中，還可能存在另一個臨界狀態，初次的屈曲模態將由穩定變為不穩定，并伴隨著屈曲波形的跳躍現象[5]，這個臨界狀態被稱為二次屈曲.上述受壓結構也可能不發生二次屈曲，初始屈曲后直接進入塑性屈服而迅速發生破壞.對于二次屈曲之前的后屈曲階段，只要給定合適的初始值或攝動參數，就可以跨越Euler 屈曲臨界值而直接利用Newton 迭代法或者攝動法進行后屈曲分析[6-7].這類方法對于屈曲波形保持不變的初始后屈曲行為是有效的.然而對于二次屈曲和彈塑性屈曲等大范圍屈曲行為，則需借助擴展方程法和延拓法等非線性屈曲分析方法實現不同平衡狀態的連續追蹤和穩定性判斷[8-9].結構非線性屈曲分析方法中的空間離散一般使用經典的三角級數[10-11]或非線性有限元格式[12]，其中除Navier 解和Levy 解之外，三角級數在大范圍后屈曲階段的截斷誤差難以保證.而基于大型有限元程序，為了僅利用切線剛度矩陣而避免其導數計算帶來的單元改進和內存增加問題，又發展了基于解流形和切向量旋轉度等幾何近似方法來計算奇異的分岔點和分岔方向[12].

本文通過結合小波Galerkin 法和不同的非線性屈曲分析方法來求解矩形薄板二次屈曲問題.首先，給出了屈曲問題控制方程的小波離散格式；然后，提出了基于一個擾動方程的切線剛度矩陣計算方法，并討論了基于小波離散方程求解二次屈曲臨界載荷的特征方程和擴展方程，以及進行后屈曲追蹤的偽弧長格式；最后，較為詳細地討論了矩形薄板完整的二次屈曲平衡路徑和波形跳躍行為.

1 薄板屈曲方程的離散格式

1.1 矩形薄板的屈曲方程

設邊長為a和b，厚度為h的各向同性彈性矩形薄板，兩對邊上分別承受面內壓縮載荷px和py以及面內剪切載荷τxy，矩形薄板屈曲問題的控制方程如下[4]：

其中w(x,y)和f(x,y)為中面橫向位移和應力函數，D和E為矩形薄板的抗彎剛度和彈性模量，?4為雙調和算子，算子H(w,f)的表達式為

面內可移的四邊簡支可等效地表示為[2]

面內可移的四邊固支可等效地表示為[2]

1.2 小波Galerkin法的離散格式

本文采用基于Coiflet 小波的Galerkin 法，其中修正Coiflet 基函數作為Galerkin法的權函數與試函數.對于定義在有限區間上的函數f(x,y)∈L2[0,1]2，有逼近公式如下：

其中二維平方可積函數空間的基函數由一維空間基函數的張量積生成，一維基函數的表達式[2,13]具體見文后附錄.為了適用于四階板殼方程，選擇小波函數消失矩數目為6 和尺度函數一階矩為7 的Coiflet 小波修正基函數.

為了便于處理控制方程中耦合導數項，首先忽略無量綱化變量的上標，且定義下列新函數：

那么方程（5）可寫為

式中Φm,i(x)表示滿足邊界條件的基函數，代入新函數的表達式且應用Galerkin 法，可得

其中系數矩陣和未知向量分別為

代入方程(8)，且應用Galerkin 法可得

2 二次屈曲分析的小波方法

2.1 Jacobi 矩陣的計算

Von Kármán 方程是個單向耦合問題，利用逆算子建立應力函數f和撓度函數w的映射關系，并代入到離散代數方程組（11）中，可化為如下非線性方程的一般形式：

其中g表示非線性算子.由向量值函數的廣義Taylor 定理，可在狀態W0附近對上式進行展開：

其中DwG(W0)和DwwG(W0)即為W0處的Jacobi 矩陣和Hesse 矩陣.

引入另一個應力函數fp，且令面內載荷表示為

那么無量綱化方程(5)可表示為

再將平衡路徑上某個平衡狀態記為(w0,f0)，且假設該狀態的擾動量為(wh,fh)，那么控制方程（15）平衡解的擾動形式則為[11]

將式（16）代入方程(15)，展開上式且考慮狀態(w0,f0)處的平衡關系，則有

化簡后可得原方程的一階擾動方程

和二階擾動方程

離散一階擾動方程（18），可得

由式(20)中第二式可得展開式（13）中的一階擾動項，即方程（12）的Jacobi 矩陣：

離散二階擾動方程(19)，可得

由式(22)中第二式可得展開式（13）中的二階擾動項，即方程（12）的Hesse 矩陣：

2.2 二次屈曲分析的方法

由隱函數定理可知，當Jacobi 矩陣DwG(W0,λx,λy,λxy)非奇異時，方程存在唯一解.而依據屈曲問題的多解性質可知其Jacobi 矩陣奇異，即有

方程（24）為非線性特征方程，一般為超越方程，直接求解非常困難.考慮僅有λx的單向屈曲問題，假設后屈曲平衡狀態w0為基本狀態，可得到由位移增量表示的線性特征方程：

因此二次屈曲的臨界條件也可表述為后屈曲路徑上某處載荷值恰好為該點處Jacobi 矩陣的特征值[10].

擴展方程法是一種適用于高維分岔問題的直接求解方法，而且能同時獲得分岔點及其特征.對于單向壓縮屈曲問題，可采用擴展方程為[9]

偽弧長法是一種具有大范圍收斂性的間接延拓方法，在分岔點附近對初始值的選取要求不高.式（12）則可進一步化為常微分方程的Cauchy 問題[12]：

式（28）為Euler 預報格式，式（29）為Newton 修正格式，迭代序列收斂于解曲線y(s)上的某一點，DG(yk)則為預報點yk處的Fréchet 導數，DG的表達式為

值得指出的是，擴展方程中的特征向量是求解分岔方向的前提，而分岔方向是進行解支轉換的前提，聯合使用擴展方程法和偽弧長法有利于開展大范圍的后屈曲行為追蹤來獲取較為真實的失穩過程.

對于多參數屈曲問題，例如px和py雙向壓縮矩形薄板，可采用比例加載的方式求解臨界參數，給定λx，λy=kλx，則臨界參數為λx,cr，λy,cr=kλx,cr，通過改變比例系數k可得雙參數平面上的相關曲線.平衡狀態穩定性由Lagrange 定理判定，即Jacobi 矩陣為負定矩陣.

3 數值算例與討論

3.1 Euler 屈曲算例

圖1 四邊簡支板屈曲載荷的收斂階數Fig.1 Convergence orders of buckling loads for 4-edge simply supported plates

對于方形板，四邊簡支情形的后屈曲平衡路徑由二次攝動法獲得[7]，四邊固支情形由非線性有限單元法和基于Bior 3.1 的小波Galerkin 法獲得其后屈曲平衡路徑[6]，如圖2所示.本文利用基于Coiflet 的小波Galerkin 法求得的后屈曲路徑與上述結果是一致的，四邊簡支情形（m=n=4，5）兩者的誤差在最大撓度2 倍范圍內誤差不超過5%，四邊固支情形（m=4）與已有結果的差別更小.基于Bior 3.1 的小波Galerkin 法涉及尺度函數的三重積分以及采用混合法處理邊界，而基于Coiflet 的小波Galerkin 法僅僅需要計算尺度函數的二重積分.需要指出的是，Euler 后屈曲分析中通常假設模態的波形不變，因此其分析范圍不宜取得過大.

圖2 方形薄板單向壓縮的后屈曲路徑：(a)四邊簡支；(b)四邊固支Fig.2 Post-buckling paths of square plates under uniaxial compression:(a)4-edge simply supported;(b)4-edge fixed

3.2 二次屈曲算例

Euler 屈曲系數按大小順序記為ki,i=1,2,···.最小特征值k1附近，平凡解由穩定變為不穩定，高階特征值ki(i>1)附近，平凡解均為不穩定.從ki處出發的非平凡解支上仍然有可能發生穩定性的變化，記該解支上的二次屈曲系數為ki2.本文采用拓展方程（26）直接求解二次屈曲臨界載荷，分解尺度分別取(m，n)=(4,4)，(5,4)，(5,5)和(6,5)時的方形薄板屈曲系數如圖3所示，表明小波Galerkin 法在求解二次屈曲臨界載荷時仍然具有良好的一致性.

圖3 四邊固支方形薄板的二次屈曲系數Fig.3 Secondary buckling coefficients of 4-edge fixed square plates

不同長寬比下單向受壓四邊簡支矩形薄板的二次屈曲系數見表1，分解尺度取m=6，n=4.表中結果顯示：對于β=2.0 的簡支板，屈曲系數k1=4.0，比值k12/k1和k22/k1分別為2.78 和1.34；對于β=2.0 的固支板，屈曲系數k1=7.467，比值k12/k1為6.15.β=2.0 的鋁板受壓屈曲實驗則顯示，簡支板在加載壓力為初始臨界壓力1.3 倍左右時屈曲波形由兩個半波迅速跳躍成三個半波，固支板從一開始就出現三個半波，直至出現塑性破壞也沒有出現波形跳躍現象[2].由半解析半經驗公式可得，β=2.0 簡支板和固支板的平均破壞應力與Euler 屈曲應力的比值分別為5.3 和3.5[14].可見，若k12對應的二次臨界載荷小于屈曲破壞載荷，則實驗能觀察到波形跳躍現象，不過觀察到的是不穩定非平凡解支上的失穩點，即k22.這是由于波形跳躍是一個變化迅速的動態過程.Nakamura 和Uetani 獲得的長寬比為2 簡支板的二次屈曲系數為12.45，更接近k12，這是由于其判定條件為穩定的非平凡解支失穩[8].由非線性有限元獲得固支方板的k12和k22的值分別為17.46 和12.53，與本文結果也較為一致[12].

表1 單向受壓彈性矩形薄板的二次屈曲系數Table 1 Secondary buckling coefficients of rectangular plates under uniaxial compression

為了確定波形跳躍的具體途徑，還應了解不同屈曲模態之間的耦合，因此需要開展更大范圍內的后屈曲平衡路徑分析.對于長寬比為1.0 和2.0 的四邊簡支和四邊固支矩形薄板，聯合使用擴展方程和偽弧長法獲得了大范圍的后屈曲平衡路徑追蹤并進行了穩定性判定，分別如圖4 和5所示.結果顯示，k12處既可能發生超臨界分岔也可能發生亞臨界分岔，k22處則一般為亞臨界分岔.四邊固支方形薄板的大范圍后屈曲路徑最為典型，k12處和k22處均發生亞臨界分岔且k12和k22之間存在連續的不穩定解，這個不穩定的解為前兩階模態的耦合，且投影到位移-轉角組成的平面上為一條閉合曲線.長寬比為2.0 的四邊固支板和簡支板，k12處分別發生亞臨界分岔和超臨界分岔.

圖4 四邊簡支薄板二次屈曲平衡路徑：(a)β=1.0；(b)β=2.0Fig.4 Secondary buckling paths of 4-edge simply supported plates:(a)β=1.0;(b)β=2.0

圖5 四邊固支薄板二次屈曲平衡路徑：(a)β=1.0；(b)β=2.0Fig.5 Secondary buckling paths of 4-edge fixed plates:(a)β=1.0;(b)β=2.0

大量的數值計算表明，從初始屈曲到二次屈曲的平衡路徑可分為兩種情形：一類是隨著載荷幅值增大，屈曲波形不發生改變，只發生波形幅值增長，如圖6 給出的方形薄板平衡曲面就是這種情形；一類是隨著載荷增大，屈曲波形和幅值均發生變化，如圖7 給出的β=2.0 矩形薄板平衡曲面則是這種情形.表1 中二次屈曲系數k12的值超過30 的矩形薄板，其二次屈曲附近的屈曲波形均會發生變化，即出現相鄰穩定屈曲波形的耦合波形.最后考慮了雙向壓縮k=0.2 下四邊固支方板的二次屈曲路徑，如圖8所示，可見載荷變化也使得屈曲模態的耦合通道不再閉合.

圖6 方形薄板（β=1.0）后屈曲平衡曲面：(a)λ=110；(b)λ=159.991Fig.6 Equilibrium surfaces of square platesβ=1.0）on postbuckling paths:(a)λ=110;(b)λ=159.991

圖7 矩形薄板（β=2.0）后屈曲平衡曲面：(a)λ=350；(b)λ=1811.765Fig.7 Equilibrium surfaces of rectangular plates (β=2.0)on postbuckling paths:(a)λ=350;(b)λ=1811.765

圖8 加載比例k=0.2 四邊固支方板雙向壓縮的二次屈曲路徑Fig.8 Secondary buckling paths of a 4-edge fixed square plate with load ratio k=0.2

4 結論

本文開展了基于Coiflet 的小波Galerkin 法在矩形薄板二次屈曲問題中的應用研究，并給出了一個非線性屈曲方程Jacobi 矩陣的簡便計算方法.數值計算結果表明：

1）基于Coiflet 的小波Galerkin 法在求解屈曲臨界載荷時仍然具有較好的收斂性，且與穩定性實驗結果一致；

2）在小范圍后屈曲平衡路徑追蹤方面，小波Galerkin 法與攝動法和有限單元法等現有結果也是一致的；

3）完善矩形薄板二次屈曲的多參數大范圍分析表明，在長寬比、邊界條件和雙向壓縮載荷的影響下其前兩階屈曲模態的閉合耦合通道較為少見，矩形薄板后屈曲狀態之間的跳躍可能會涉及更高階的穩定模態.

附錄

文獻[2,13]中建立的基于Coiflet 尺度函數的修正基函數表達式為

修正系數為

引入 β0和 β1是為說明邊界條件的施加過程，為了滿足x=0 處的齊次邊界，僅僅需要令β0,i=0，而其他元素保持不變即可.根據四點Malkov 數值微分公式推導了 ζ0和 ζ1：