付為剛， 廖 喆， 熊煥杰， 馬駿馳

(中國民用航空飛行學(xué)院 航空工程學(xué)院， 四川 廣漢 618307)

0 引言

板殼結(jié)構(gòu)在航空航天[1,2]、武器裝備[3,4]以及包裝制品[5]等領(lǐng)域應(yīng)用廣泛.板殼結(jié)構(gòu)在受到的壓縮載荷達(dá)到臨界值時(shí)會(huì)發(fā)生屈曲失穩(wěn)破壞，從而造成災(zāi)難性的事故.作為影響結(jié)構(gòu)安全性的重要設(shè)計(jì)因素，板殼結(jié)構(gòu)的屈曲失穩(wěn)特性一直是人們關(guān)注的熱點(diǎn).

近年來，許多學(xué)者對(duì)含簡支與固支邊界條件下矩形板的屈曲失穩(wěn)特性進(jìn)行了深入研究[6,7].Wang等[8]采用辛疊加方法得出了四種簡支與固支組合邊界條件下的解析解，其結(jié)果可作為其它數(shù)值方法的基準(zhǔn).Salama[9]采用Raylei-Ritz法研究了端部和中部承受單向軸壓作用下加載邊簡支、非加載邊固支矩形板的屈曲失穩(wěn)特性.該方法依賴于事先假定滿足邊界條件的撓度函數(shù).葛威延等[10]采用高階有限條傳遞矩陣法，分析了兩對(duì)邊簡支、兩對(duì)邊固支邊界條件下四邊受壓矩形板的屈曲失穩(wěn)問題，研究結(jié)果表明，該方法在有限條傳遞矩陣法的基礎(chǔ)上進(jìn)一步提高了計(jì)算精度.然而，在復(fù)雜的工程應(yīng)用環(huán)境中，矩形板四邊也存在含自由邊界條件的現(xiàn)象，研究含有自由邊界條件矩形板的屈曲失穩(wěn)問題具有重要意義.鮑四元等[11]基于最小勢能原理，計(jì)算了單向軸壓作用下任意彈性邊界矩形板的屈曲載荷系數(shù)，其中包括三邊固支、一邊自由等包含自由邊界的情況，并給出了模擬不同邊界條件時(shí)彈簧剛度的合理取值.劉寅華等[12]采用能量法以及非線性有限元法對(duì)加載對(duì)邊簡支，非加載邊一邊固支、一邊自由的矩形薄板的屈曲強(qiáng)度進(jìn)行了分析，并對(duì)比歐拉失穩(wěn)和極值點(diǎn)失穩(wěn)的臨界應(yīng)力，發(fā)現(xiàn)兩者差異不大.該方法僅適用于特定邊界條件，并不具有普適性.于海軍等[13]利用辛彈性力學(xué)方法推導(dǎo)出了矩形帶材局部縱向屈曲區(qū)域承受加載對(duì)邊簡支、非加載邊一邊簡支、一邊自由約束條件時(shí)的臨界屈曲應(yīng)力和屈曲撓度函數(shù).而Li等[14,15]采用辛疊加方法，先后得出了四邊自由和具有兩個(gè)相鄰自由邊的雙軸壓縮矩形薄板屈曲失穩(wěn)問題的解析解.與通常采用的半逆解法相比，辛彈性力學(xué)方法無需對(duì)解的形式進(jìn)行假設(shè)，但其推導(dǎo)過程較為復(fù)雜.

有限差分法是一種數(shù)值解法，其有著簡單、通用性強(qiáng)及在計(jì)算機(jī)程序上實(shí)現(xiàn)簡便等優(yōu)勢.Ravari等[16]使用有限差分法研究了矩形納米板的屈曲失穩(wěn)特性，并與文獻(xiàn)解進(jìn)行比較，驗(yàn)證了有限差分法的有效性.袁堅(jiān)鋒等[17]基于Levy形式級(jí)數(shù)和有限差分法提出了面內(nèi)彎曲載荷作用下兩邊簡支、兩邊固支復(fù)合材料層合板的數(shù)值解法.該方法理論推導(dǎo)復(fù)雜，且只適用于特定邊界條件，不能適用于包含自由邊界條件矩形板的屈曲失穩(wěn)問題求解.

目前，針對(duì)含對(duì)邊自由邊界矩形板屈曲失穩(wěn)特性方面的研究工作較少.本文基于有限差分法對(duì)含對(duì)邊自由邊界矩形板的屈曲失穩(wěn)特性展開研究.首先，針對(duì)不考慮剪切效應(yīng)的矩形板屈曲控制微分方程及邊界條件進(jìn)行離散化，將求解臨界屈曲載荷的問題轉(zhuǎn)換為求解撓度系數(shù)矩陣廣義特征值的問題；對(duì)比不同網(wǎng)格結(jié)點(diǎn)數(shù)下有限差分?jǐn)?shù)值解、有限元仿真解與文獻(xiàn)解之間的相對(duì)誤差，驗(yàn)證本文求解方法的精確性.在此基礎(chǔ)上，系統(tǒng)研究矩形板幾何參數(shù)中邊長比與厚寬比同有限差分法求解精度之間的耦合作用規(guī)律.針對(duì)厚寬比較大的矩形板，采用最小二乘法進(jìn)行曲線擬合給出有限差分法精確求解的幾何參數(shù)適用范圍計(jì)算公式.最后，經(jīng)過實(shí)例參數(shù)分析驗(yàn)證表明，有限差分法能在擬合公式限定的幾何參數(shù)范圍內(nèi)精確求解含對(duì)邊自由邊界矩形板的屈曲失穩(wěn)問題.

1 基本方程

受軸向壓縮載荷作用下含對(duì)邊自由邊界矩形板的示意圖如圖1所示.在圖1中：b為矩形板的寬度；l為矩形板的長度；F代表自由邊界條件；Nx為其x方向上的軸向載荷.

基于彈性板的小撓度彎曲理論，各向同性矩形板的屈曲控制微分方程[18]：

(1)

式(1)中：D=Et3/[12(1-v2)]為矩形板的彎曲剛度；E為板的彈性模量；v為泊松比；t為板厚；w為撓度；Nx、Ny以及Nxy分別為x方向上的軸向載荷、y方向的上的軸向載荷、xy平面內(nèi)的剪切載荷.

考慮矩形板僅承受x方向上的均勻分布軸向壓縮載荷，即Ny=Nxy=0，則式(1)簡化為

(2)

圖1 含對(duì)邊自由邊界矩形板受載示意圖

矩形板的典型邊界條件包括固支、簡支以及自由邊界條件，圖1中邊AC、BD為自由邊界，邊AB、CD受簡支或固支邊界條件約束.以邊x=0為例，以上三種邊界條件可表示為

簡支邊：w=0,Mx=0

2 控制方程的有限差分法求解

有限差分法是應(yīng)用離散的思想，用一個(gè)插值多項(xiàng)式及其微分來代替偏微分方程的解，首先需對(duì)求解區(qū)域進(jìn)行網(wǎng)格劃分，然后對(duì)方程中的導(dǎo)數(shù)進(jìn)行充分近似，最后得到網(wǎng)格結(jié)點(diǎn)上未知函數(shù)的線性代數(shù)方程組并對(duì)其求解.

采用有限差分法對(duì)簡化后的屈曲控制微分方程式(2)進(jìn)行離散化處理.分別沿x、y方向?qū)⒕匦伟鍏^(qū)域均勻劃分成長度為Δx、寬度為Δy的網(wǎng)格，如圖2所示.對(duì)于矩形板內(nèi)任一結(jié)點(diǎn)的撓度wi,j的各階差分公式可通過圖2中相應(yīng)的12個(gè)結(jié)點(diǎn)處的撓度按式(3)～(6)表示

一階：

(3)

二階：

(4)

三階：

(5)

四階：

(6)

圖2 結(jié)點(diǎn)劃分示意圖

將式(3)～(6)代入式(2)，得到板內(nèi)各個(gè)結(jié)點(diǎn)處的以撓度為未知量的差分方程式(7).同時(shí)，根據(jù)矩形板及其邊界條件的對(duì)稱性，可對(duì)結(jié)點(diǎn)進(jìn)行簡化，減少獨(dú)立結(jié)點(diǎn)數(shù)量.

i=1,2,…,n;j=1,2,…,m

(7)

式(7)中：n為矩形板y方向上結(jié)點(diǎn)數(shù)，m為矩形板x方向上結(jié)點(diǎn)數(shù).

下文以軸向壓縮載荷作用下非加載對(duì)邊自由、加載對(duì)邊簡支矩形板為例，研究含對(duì)邊自由邊界矩形板的屈曲失穩(wěn)特性.通過式(3)～(6)導(dǎo)出邊界條件的差分形式：

自由邊AC、BD：

(8(a))

(8(b))

對(duì)于簡支邊AB、CD，Mx=0可以簡化為?2w/?x2=0:

wi,j=0

(9(a))

(9(b))

考慮到矩形板邊界處結(jié)點(diǎn)的差分方程中包含板外虛結(jié)點(diǎn)撓度，應(yīng)用邊界條件，將板外虛結(jié)點(diǎn)用板內(nèi)結(jié)點(diǎn)表示.

如圖3所示，S、F分別代表簡支、自由邊界條件.考慮矩形板的兩種邊界條件，對(duì)自由邊外由內(nèi)到外第一行虛結(jié)點(diǎn)上任意結(jié)點(diǎn)4，第二行虛結(jié)點(diǎn)上任意結(jié)點(diǎn)12，以及簡支邊外第一列虛結(jié)點(diǎn)上任意結(jié)點(diǎn)20處的撓度可如式(10)表示

(10)

將式(10)帶入差分方程式(7)，得到一個(gè)線性方程組，其矩陣形式為

(A+ηB)w=0

(11)

式(11)中：η=Nx/(DΔx2)；其中矩陣A、B由邊界條件、網(wǎng)格劃分、彎曲剛度等決定.求得式(11)撓度系數(shù)矩陣廣義特征值的最小值ηcr，進(jìn)而得到矩形板的臨界屈曲載荷：

Ncr=ηcrDΔx2

(12)

圖3 邊界外虛結(jié)點(diǎn)與邊界內(nèi)結(jié)點(diǎn)示意圖

3 收斂性分析

基于MATLAB平臺(tái)編程實(shí)現(xiàn)有限差分法求解矩形板的臨界屈曲載荷，程序運(yùn)行流程如圖4所示.

圖4 MATLAB有限差分法程序流程圖

計(jì)算中取x、y方向上的結(jié)點(diǎn)數(shù)相等(m=n).單向軸壓作用下非加載對(duì)邊自由、加載對(duì)邊簡支矩形板寬度b=1 m，厚寬比t/b=0.01.本文材料彈性模量E=70 GPa，泊松比v=0.3.分別采用有限差分法、有限元軟件ABAQUS的buckle模塊求解其臨界屈曲載荷.ABAQUS有限元模型(r=1)如圖5所示.ABAQUS有限元模型采用4結(jié)點(diǎn)減縮積分殼單元(S4R)，每個(gè)單元包含4個(gè)結(jié)點(diǎn)，在每個(gè)結(jié)點(diǎn)處有6個(gè)自由度(3個(gè)平動(dòng)自由度，3個(gè)轉(zhuǎn)動(dòng)自由度).對(duì)于矩形板的左右兩端采取簡支約束并在右端施加軸向壓縮載荷，上下兩端不設(shè)置約束，使其作為自由邊.

圖5 有限元模型

當(dāng)邊長比r=1時(shí)，基于有限元法在不同結(jié)點(diǎn)數(shù)下求解單向軸壓作用下非加載對(duì)邊自由、加載對(duì)邊簡支矩形板的臨界屈曲載荷與文獻(xiàn)解之間的相對(duì)誤差，結(jié)果如圖6所示.隨著結(jié)點(diǎn)數(shù)的增加，有限元仿真解逐漸收斂，結(jié)點(diǎn)數(shù)10 201、12 544分別對(duì)應(yīng)單元尺寸的取值為0.01 m、0.009 m，當(dāng)結(jié)點(diǎn)數(shù)取上述兩值時(shí)，由圖6可見，基于文獻(xiàn)[19]求解結(jié)果的相對(duì)誤差均為0.44%，而基于文獻(xiàn)[20]求解結(jié)果的相對(duì)誤差僅由0.36%變到0.35%，此時(shí)有限元仿真解具有較高精確性及穩(wěn)定性，由此認(rèn)為單元尺寸取值為0.01 m(結(jié)點(diǎn)數(shù)為10 201)時(shí)有限元仿真解收斂，后續(xù)有限元仿真相關(guān)計(jì)算均采用此單元尺寸.

圖6 有限元仿真解相對(duì)誤差隨結(jié)點(diǎn)數(shù)變化曲線

當(dāng)邊長比r=1時(shí)，基于有限差分法在不同結(jié)點(diǎn)數(shù)下求解單向軸壓作用下非加載對(duì)邊自由、加載對(duì)邊簡支矩形板的臨界屈曲載荷與文獻(xiàn)解之間的相對(duì)誤差，結(jié)果如圖7所示.隨著結(jié)點(diǎn)數(shù)的增加，有限差分?jǐn)?shù)值解呈收斂趨勢，基于文獻(xiàn)解的相對(duì)誤差的絕對(duì)值隨著結(jié)點(diǎn)數(shù)增加而逐漸減小，當(dāng)結(jié)點(diǎn)數(shù)取20×20(m×n)時(shí)，基于文獻(xiàn)[19]、[20]求解結(jié)果的相對(duì)誤差分別為-0.15%、-0.23%.隨著結(jié)點(diǎn)數(shù)的增加有限差分?jǐn)?shù)值解精度逐漸提高，同時(shí)曲線斜率逐漸減小，當(dāng)結(jié)點(diǎn)數(shù)由20×20變化到22×22時(shí)，各曲線中相對(duì)誤差絕對(duì)值的變化值最大僅為0.04%，而后者所需計(jì)算時(shí)間卻為前者的133.67%，為了在保證計(jì)算精度的前提下節(jié)約計(jì)算成本，后續(xù)計(jì)算中有限差分法離散結(jié)點(diǎn)數(shù)均取20×20.有限元仿真求解中，單元尺寸取值為0.01 m時(shí)包含10 201(r=1)個(gè)結(jié)點(diǎn)，且在每個(gè)結(jié)點(diǎn)處有6個(gè)自由度，其收斂速度以及計(jì)算效率遠(yuǎn)不如本文有限差分法求解程序.

圖7 有限差分?jǐn)?shù)值解相對(duì)誤差隨結(jié)點(diǎn)數(shù)變化曲線

為進(jìn)一步驗(yàn)證結(jié)點(diǎn)數(shù)取20×20時(shí)，有限差分法的適用性，在不同幾何參數(shù)下，將有限差分?jǐn)?shù)值解與文獻(xiàn)解進(jìn)行對(duì)比.如表1所示，有限差分?jǐn)?shù)值解與文獻(xiàn)[19]、文獻(xiàn)[20]求解結(jié)果非常接近.

表1 不同幾何參數(shù)下的臨界屈曲載荷(N/m)

有限元仿真解易于獲得，且能保證較高的精度，后文所述相對(duì)誤差均基于有限元仿真解.當(dāng)邊長比0.1≤r≤5時(shí)，分別采用有限差分法與有限元法求得的臨界屈曲載荷如圖8所示.兩種方法所得結(jié)果非常接近，有限差分?jǐn)?shù)值解與有限元仿真解相對(duì)誤差絕對(duì)值的最大值僅為0.84%，說明有限差分法具有較高的精確性.當(dāng)邊長比r=1時(shí)，有限元法得出臨界屈曲載荷為60 467 N/m，有限差分法得出臨界屈曲載荷為60 109 N/m，有限差分?jǐn)?shù)值解相對(duì)誤差絕對(duì)值僅為0.59%.邊長比r=1時(shí)，有限差分法求解結(jié)果與有限元模型的一階屈曲失穩(wěn)模態(tài)對(duì)比如圖9所示.由圖9可見，有限差分法求解的屈曲失穩(wěn)模態(tài)與有限元模型一致.

圖8 臨界屈曲載荷隨邊長比變化曲線

圖9 t/b=0.01,r=1時(shí)矩形板的一階屈曲模態(tài)圖

4 誤差分析

矩形板寬度b=1，厚寬比t/b分別取0.01～0.1，當(dāng)邊長比0.1≤r≤5時(shí)，以有限差分法求解單向軸壓作用下非加載對(duì)邊自由、加載對(duì)邊簡支矩形板為例研究有限差分?jǐn)?shù)值解相對(duì)誤差的變化規(guī)律，并考慮有限差分法精確求解矩形板屈曲臨界載荷的幾何參數(shù)適用范圍.

4.1 幾何參數(shù)變化對(duì)有限差分法求解精度的影響

不同厚寬比下有限差分?jǐn)?shù)值解相對(duì)誤差隨邊長比r變化曲線如圖10所示.結(jié)果顯示，隨著邊長比r的增大，有限差分法求解精度逐漸提高，當(dāng)邊長比r取值較小時(shí)，相對(duì)誤差隨邊長比r增大而減小趨勢比較劇烈；當(dāng)邊長比r取值較大時(shí)，相對(duì)誤差隨邊長比r變化趨勢較為平緩；當(dāng)邊長比r取值相同時(shí)，相對(duì)誤差隨著厚寬比的增大而增大.這是由于求解時(shí)忽略了板厚方向上的剪切效應(yīng)，而隨著厚寬比的增加，其影響也逐漸增大直至不可忽略，導(dǎo)致了相對(duì)誤差的增大.從圖10(a)可以看出，厚寬比取值較小時(shí)，相對(duì)誤差隨邊長比r的增大先單調(diào)減小再單調(diào)增大，且相對(duì)誤差會(huì)小于0，這是由于有限差分?jǐn)?shù)值解偏小，如圖7所示，結(jié)點(diǎn)數(shù)的增加可以減弱其影響，而厚寬比取值較大時(shí)忽略剪切效應(yīng)的影響占主導(dǎo)，所以未呈現(xiàn)此現(xiàn)象，如圖10(b)所示.

圖10 邊長比與厚寬比同相對(duì)誤差之間的 耦合作用規(guī)律

4.2 幾何參數(shù)精確求解適用范圍計(jì)算公式的擬合

采用曲線擬合的最小二乘法得出各厚寬比下有限差分?jǐn)?shù)值解相對(duì)誤差與邊長比的關(guān)系式，基于以上關(guān)系式得到一系列使相對(duì)誤差略小于5%的矩形板幾何參數(shù).基于這一系列幾何參數(shù)再次進(jìn)行擬合，得到有限差分法精確求解矩形板屈曲失穩(wěn)特性對(duì)應(yīng)的幾何參數(shù)適用范圍計(jì)算公式：

(t/b)=0.126 6r+0.004 98

(13)

在已知厚寬比或邊長比情況下，通過式(13)得出的矩形板幾何參數(shù)可保證相對(duì)誤差小于5%.為驗(yàn)證式(13)的準(zhǔn)確性，取厚寬比分別為0.01，0.05，0.1，基于式(13)分別求得相應(yīng)邊長比；取邊長比分別為0.5，1，2，5，基于式(13)求得相應(yīng)厚寬比.記錄以上式(13)計(jì)算結(jié)果及其相對(duì)誤差.

不同厚寬比下有限差分?jǐn)?shù)值解相對(duì)誤差隨邊長比變化曲線如圖11所示，式(13)計(jì)算得出各點(diǎn)相對(duì)誤差均略小于5%，當(dāng)邊長比大于式(13)計(jì)算結(jié)果時(shí)，相對(duì)誤差均在5%以下.圖12顯示了不同邊長比下的有限差分?jǐn)?shù)值解相對(duì)誤差隨厚寬比變化曲線，同樣式(13)計(jì)算得出各點(diǎn)相對(duì)誤差均略小于5%，當(dāng)厚寬比小于式(13)計(jì)算結(jié)果時(shí)，相對(duì)誤差不超過5%.綜上，式(13)可精確計(jì)算有限差分法精確求解矩形板屈曲失穩(wěn)問題的幾何參數(shù)適用范圍.

圖11 當(dāng)厚寬比t/b=0.01,0.05,0.1時(shí),相對(duì)誤差 隨邊長比變化曲線

圖12 當(dāng)邊長比r=0.5,1,2,5時(shí),相對(duì)誤差 隨厚寬比變化曲線

5 結(jié)論

本文采用有限差分法，以單向軸壓作用下非加載對(duì)邊自由、加載對(duì)邊簡支矩形板為例，對(duì)含對(duì)邊自由邊界矩形板的屈曲失穩(wěn)特性進(jìn)行了研究，得到以下主要結(jié)論：

(1)有限差分?jǐn)?shù)值解隨結(jié)點(diǎn)數(shù)的增加逐漸收斂，且具有較高計(jì)算效率.通過與文獻(xiàn)解、精確的有限元仿真解之間的對(duì)比分析，驗(yàn)證了有限差分法的精確性.

(2)當(dāng)厚寬比一定時(shí)，有限差分法求解的相對(duì)誤差會(huì)隨著邊長比的增大而減??；當(dāng)邊長比一定時(shí)，由于剪切效應(yīng)的影響有限差分?jǐn)?shù)值解相對(duì)誤差會(huì)隨著矩形板厚寬比增大而增大.

(3)有限差分法能夠在擬合公式給出的幾何參數(shù)適用范圍內(nèi)，精確求解含對(duì)邊自由邊界矩形板的屈曲失穩(wěn)問題，使其相對(duì)誤差保持在給定范圍內(nèi).