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

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

引言

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

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

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

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

1.1 矩形薄板的屈曲方程

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

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

面內(nèi)可移的四邊簡支可等效地表示為[2]

面內(nèi)可移的四邊固支可等效地表示為[2]

1.2 小波Galerkin法的離散格式

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

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

為了便于處理控制方程中耦合導(dǎo)數(shù)項(xiàng)，首先忽略無量綱化變量的上標(biāo)，且定義下列新函數(shù)：

那么方程（5）可寫為

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

其中系數(shù)矩陣和未知向量分別為

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

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

2.1 Jacobi 矩陣的計(jì)算

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

其中g(shù)表示非線性算子.由向量值函數(shù)的廣義Taylor 定理，可在狀態(tài)W0附近對上式進(jìn)行展開：

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

引入另一個應(yīng)力函數(shù)fp，且令面內(nèi)載荷表示為

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

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

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

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

和二階擾動方程

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

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

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

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

2.2 二次屈曲分析的方法

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

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

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

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

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

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

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

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

3 數(shù)值算例與討論

3.1 Euler 屈曲算例

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

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

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

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

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

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

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

圖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

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

圖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 結(jié)論

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

1）基于Coiflet 的小波Galerkin 法在求解屈曲臨界載荷時仍然具有較好的收斂性，且與穩(wěn)定性實(shí)驗(yàn)結(jié)果一致；

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

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

附錄

文獻(xiàn)[2,13]中建立的基于Coiflet 尺度函數(shù)的修正基函數(shù)表達(dá)式為

修正系數(shù)為

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