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

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

0 引言

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

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

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

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

1 基本方程

受軸向壓縮載荷作用下含對邊自由邊界矩形板的示意圖如圖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 含對邊自由邊界矩形板受載示意圖

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

簡支邊：w=0,Mx=0

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

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

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

一階：

(3)

二階：

(4)

三階：

(5)

四階：

(6)

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

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

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

(7)

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

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

自由邊AC、BD：

(8(a))

(8(b))

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

wi,j=0

(9(a))

(9(b))

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

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

(10)

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

(A+ηB)w=0

(11)

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

Ncr=ηcrDΔx2

(12)

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

3 收斂性分析

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

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

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

圖5 有限元模型

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

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

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

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

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

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

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

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

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

4 誤差分析

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

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

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

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

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

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

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

(13)

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

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

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

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

5 結(jié)論

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

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

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

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