一種新型板彎單元及其在幾何非線性下的應(yīng)用

黃志強，周 云

(1.太原科技大學(xué)應(yīng)用科學(xué)學(xué)院數(shù)學(xué)系，太原 030024;2.西北工業(yè)大學(xué)理學(xué)院，西安 710129;3.西安應(yīng)用光學(xué)研究所，西安 710065)

無人機因其具有成本低、機動性好、能實現(xiàn)“零傷亡”等諸多優(yōu)良特性而被廣泛使用，世界各航空強國均投入大量人力、財力以加速先進無人機系統(tǒng)的研制。從機體平臺方面來看，結(jié)構(gòu)輕質(zhì)化、低可探測性、高機動性、長航時、遠航程等特性成為無人機研制的主要技術(shù)要求和未來發(fā)展方向[1-2]。從結(jié)構(gòu)輕質(zhì)化角度考慮，由于先進復(fù)合材料具有比強度高、比剛度高、抗疲勞能力強等優(yōu)良性質(zhì)，故世界上先進的無人機均大量使用復(fù)合材料。在大展弦比機翼結(jié)構(gòu)設(shè)計時，由于機翼展向尺寸大，故翼根將承受較大彎矩，同時也要求設(shè)計人員仔細考慮機翼蒙皮由于彎曲引起的變形和應(yīng)力。對金屬蒙皮而言，應(yīng)用經(jīng)典的板殼理論(如Kirchhoff薄板理論)可較為精確地分析其彎曲變形及應(yīng)力;而對復(fù)合材料蒙皮而言，使用經(jīng)典的彎曲變形理論對其進行分析將導(dǎo)致較大誤差，故有必要使用更加精細、恰當(dāng)?shù)膹澢Ｐ蛠砻枋鰪?fù)合材料板殼結(jié)構(gòu)的變形。

本文將介紹一種高階剪切變形理論及與之對應(yīng)的四邊形有限元，并將該單元用于復(fù)合材料層合板的應(yīng)力分析、變形模擬和結(jié)構(gòu)的幾何非線性分析。數(shù)值算例表明這種新型單元保證非線性有限元算法的收斂性，能夠比較精確計算復(fù)合材料板受彎時產(chǎn)生的撓度和應(yīng)力。

1 Reddy高階剪切變形理論

較為熟知的板彎曲理論有Kirchhoff薄板理論及一階剪切變形理論。前者很好地滿足力邊界條件，但對中厚板而言，橫向剪切變形的影響不應(yīng)忽略;而對于復(fù)合材料層合板，當(dāng)面內(nèi)彈性模量之比較大(如E11/E22＞25)時，即使其寬厚比滿足薄板條件，橫向剪切變形的影響也不容忽略。為將橫向剪切變形的影響考慮進來，Reissner和Mindlin等人提出了一階剪切變形理論。但該理論不能精確滿足力邊界條件，需要引入剪切因子作為修正。對于復(fù)合材料層合板，剪切因子和每一層的材料參數(shù)、鋪層方式、幾何形狀、邊界條件以及其它因素有關(guān)，確定起來相當(dāng)麻煩[3]。

為了不引進剪切因子并保證板上、下表面的力邊界條件得到滿足，種種高階剪切變形理論被提出，其中Reddy提出的高階剪切變形理論非常有代表性[4](以下稱“Reddy高階剪切變形理論”)。該理論假設(shè)位移函數(shù)為如下形式:

其中 ui(i=1，2，3)分別為某一點處 x、y、z方向的總位移;u、v代表參考面上x、y方向位移，一般由縱向載荷引起;θx和θy為原來垂直中面的纖維在彎后形成的曲線在中面處的轉(zhuǎn)角。不難驗證由式(1)所示的位移模式導(dǎo)出的橫向剪應(yīng)變在板上、下表面為零，即:

結(jié)合復(fù)合材料層合板的本構(gòu)方程可知在板上、下表面橫向剪應(yīng)力也為零，這就達到了無需引入剪切因子就能保證板上、下表面力邊界條件精確滿足的目的。

3 用于幾何非線性計算的HSQ8P單元

由于式(1)所示的位移模式含有撓度的一階偏導(dǎo)數(shù)，這要求在協(xié)調(diào)的位移型有限元設(shè)計時撓度是C1連續(xù)的，即撓度本身及其偏導(dǎo)數(shù)在計算區(qū)域內(nèi)連續(xù):

對四邊形單元而言，構(gòu)造C1連續(xù)的撓度插值非常困難。于是比式(2)更弱的連續(xù)性條件被用來構(gòu)造撓度，其中一種思路為加權(quán)積分意義下的撓度協(xié)調(diào)。

式(3)中 ρi(i=1，2，3)為權(quán)函數(shù)，可通過對應(yīng)變能泛函進行變分來確定。從式(3)出發(fā)，文獻[5]針對四邊形網(wǎng)格提出了一種新型撓度插值。以這種撓度插值為基礎(chǔ)，文獻[6]構(gòu)造了一種新單元HSQ8P.該單元很好地解決了Reddy高階剪切變形理論位移型有限元設(shè)計時遇到的撓度協(xié)調(diào)性問題。數(shù)值算例顯示，在小變形條件下，當(dāng)網(wǎng)格加密時計算結(jié)果收斂于精確解，在四邊形網(wǎng)格下HSQ8P元有較高的計算精度。

現(xiàn)考慮將HSQ8P元用于幾何非線性條件下的變形、應(yīng)力分析。對板彎曲問題而言，幾何非線性效應(yīng)主要來自結(jié)構(gòu)較大的剛體平動和剛體轉(zhuǎn)動，而應(yīng)變往往較小。這種情況下為便于有限元計算，可采用 Green應(yīng)變和第二類 Piola-Kirchhoff應(yīng)力(PK2)來描述應(yīng)變和應(yīng)力，它們之間的本構(gòu)關(guān)系和小位移條件下本構(gòu)關(guān)系完全一致，相應(yīng)的計算格式為完全 Lagrange格式(T.L.格式)。類似于文獻[7]針對實體單元的分析，現(xiàn)推導(dǎo)HSQ8P元用于幾何非線性計算的T.L.格式有限元方程。

對應(yīng)于T.L.格式的虛功原理為:外力對虛增量位移做的功等于虛應(yīng)變能增量。假設(shè)對第k載荷步的計算已經(jīng)完成，則第k+1步的應(yīng)變?yōu)?

其中Ek為第k步的應(yīng)變，可由第k步位移表出;虛應(yīng)變增量可記為:

結(jié)合式(1)不難給出Bk的顯式表達。和應(yīng)變一樣，應(yīng)力也可以寫成增量形式:Sk+1=Sk+△Sk.故在單元上虛功方程為:

應(yīng)力、應(yīng)變之間滿足關(guān)系:

其中D矩陣和幾何線性情況下所使用的本構(gòu)關(guān)系矩陣一致。將式(6)按單元組裝即可得到有限元方程。

值得指出的是，由于有限元方程中的系數(shù)矩陣涉及到應(yīng)力的計算，故位移解的結(jié)果也間接反應(yīng)出應(yīng)力的計算精度。

3 數(shù)值算例

算例1 復(fù)合材料層合板應(yīng)力計算[8]

正方形層合板幾何尺寸、網(wǎng)格剖分如圖1所示:h1=h3=h/4;三層的鋪層角為分別為0°、90°、0°;材料常數(shù):E11/E22=25，G12=G13=0.5E22，G23=0.2E22，ν12=0.25;邊界條件為四邊簡支;載荷為雙正弦分布q=q0sin(πx/a)sin(πy/a).某些關(guān)鍵點處的應(yīng)力和撓度數(shù)值解與精確解之比見表1.從表1不難看出，基于Reddy高階剪切變形理論的單元HSQ8P能精確地求解層合板的撓度和應(yīng)力。

圖1 網(wǎng)格剖分及層合板形狀Fig.1 Mesh and shape of the laminated plate

表1 應(yīng)力和撓度數(shù)值解與精確解之比/%Tab.1 Numerical solutions of stress and deflection with the exact solution ratio/%

算例2 各向同性方板的非線性彎曲

計算正方形板在縱向拉伸載荷、橫向彎曲載荷耦合作用下的非線性變形。板邊長為1，厚度0.2，剖分成16個單元;E=1000，ν=0.3;一端固支，另一端的三個節(jié)點施加拉伸-彎曲耦合載荷:f拉=f彎=0.3(如圖2 所示)。

圖2 正方形板的網(wǎng)格剖分、邊界條件、載荷示意圖Fig.2 Square board，mesh，boundary conditions and load diagram

表2 拉-彎載荷下方板自由端中點軸向位移(×10－3)Tab.2 Pull-bending loads below the plate free end of the mid-point axial displacement(×10－3)

表3 拉-彎載荷下方板自由端中點橫向位移(×10－2)Tab.3 Pull-under the bending load plate free end of the midpoint of the horizontal displacement(×10－2)

表2-表3給出了HSQ8P元以及商用FEA軟件Nastran進行非線性計算的結(jié)果。數(shù)據(jù)顯示HSQ8P有很高的計算精度。

4 結(jié)論

將HSQ8P單元的應(yīng)用范圍從小變形條件下位移求解擴展到應(yīng)力計算、大變形分析，數(shù)值算例顯示:HSQ8P單元具有較高的位移、應(yīng)力數(shù)值精度，當(dāng)板結(jié)構(gòu)發(fā)生大變形時也能保證結(jié)果的收斂性和計算精度。

[1]尹晶.無人機及其復(fù)合材料結(jié)構(gòu)優(yōu)化設(shè)計[D].西安:西北工業(yè)大學(xué)，2009.

[2]馬永忠，崔小朝，陶海瑞，等.復(fù)合材料錨桿托板的有限元分析與實驗研究[J].太原科技大學(xué)學(xué)報，2007，28(4):327-330.

[3]沈惠申.板殼后屈曲行為[M].上海:上?？茖W(xué)技術(shù)出版社，2002.

[4]REDDY J N.A refined nonlinear theory of plates with transverse shear deformation[J].International Journal of Solids and Structures，1984，20:881-896.

[5]周云，孫秦.一種四邊形廣義協(xié)調(diào)薄板彎曲單元[J].計算機仿真，2010，27(3):322-325.

[6]周云，孫秦.一種基于高階剪切變形板模型的有限元算法[J].航空工程進展.2010，1(1):71-75.

[7]凌道盛，徐興.非線性有限元及程序[M].杭州:浙江大學(xué)出版社，2004.

[8]PAGANO N J，HATFIELD S J.Elastic behavior of multilayered bidirectional composites[J].AIAA Journal，1972，10:931-933.