復(fù)合材料層合板的雜交有限元方法

卿光輝，賈瑞升（中國民航大學(xué)航空工程學(xué)院，天津 300300）

復(fù)合材料層合板的雜交有限元方法

卿光輝，賈瑞升
（中國民航大學(xué)航空工程學(xué)院，天津 300300）

結(jié)合復(fù)合材料修正后的H-R混合變分原理，直接借助應(yīng)力-應(yīng)變關(guān)系，推導(dǎo)了新的應(yīng)力模式，建立了復(fù)合材料層合板的雜交等參有限元列式。利用Mathematica語言編程進(jìn)行數(shù)值實(shí)例分析，其計(jì)算結(jié)果與相關(guān)文獻(xiàn)的精確解以及Abaqus軟件建模分析結(jié)果對比，實(shí)例證明該方法所得到的各個(gè)靜力學(xué)量更接近精確解，并且可用較少的網(wǎng)格劃分得到較精確的解。

復(fù)合材料；雜交有限元；等參元；Mathematica

近年來復(fù)合材料由于其結(jié)構(gòu)本身的特性，已經(jīng)在航空領(lǐng)域得到了廣泛的應(yīng)用和發(fā)展。20世紀(jì)60年代初期假設(shè)位移法對于平板彎曲的問題遇到了困難，主要是假定的形函數(shù)不容易在邊界上滿足C1連續(xù)性條件，即可以保證位移的連續(xù)性，但不容易保證轉(zhuǎn)角的連續(xù)性。而雜交應(yīng)力有限元以假定單元的應(yīng)力試解為基本特色，很容易實(shí)現(xiàn)單元間位移協(xié)調(diào)性的要求，這對于克服板殼中要求撓度一階導(dǎo)數(shù)連續(xù)的困難有重要的意義。1980年以后，Chen和T.H.H.Pian[1]運(yùn)用等參坐標(biāo)來表示單元應(yīng)力，使用Hellinger-Reissner原理推導(dǎo)時(shí)，可以不必考慮平衡條件[2-4]。本文在前人工作的基礎(chǔ)上，直接借助應(yīng)力-應(yīng)變關(guān)系，推導(dǎo)了新的應(yīng)力模式，建立了復(fù)合材料層合板的雜交有限元和等參元及其數(shù)值分析模型，并利用Mathematica語言編程進(jìn)行了具體數(shù)值實(shí)例的分析，數(shù)值結(jié)果與Abaqus有限元軟件模型和相關(guān)文獻(xiàn)的精確解進(jìn)行了對比，發(fā)現(xiàn)本文的方法對網(wǎng)格的劃分不是很敏感，即網(wǎng)格劃分?jǐn)?shù)目較少就可以得到比較精確的解，并且隨著網(wǎng)格的加密計(jì)算所得的結(jié)果變化不大。

1 變分原理和狀態(tài)方程

1.1 修正的H-R變分原理

假設(shè)復(fù)合材料是各向同性或正交各向異性的，其本構(gòu)關(guān)系可由下式給出

其中：σx、σy、σz、σyz、σxz、σxy分別為應(yīng)力分量；C11、C12、C13、C22、C23、C33、C44、C55、C66是復(fù)合材料的剛度系數(shù)；Sx、Sy、Sz、Syz、Sxz、Sxy是應(yīng)變分量。

應(yīng)變-位移關(guān)系為

其中：u、v、w為沿x、y、z方向的位移分量。

在單元內(nèi)假設(shè)獨(dú)立的應(yīng)力場

其中：P={σ1σ2… σM}為應(yīng)力矩陣；σi為假設(shè)應(yīng)力模式；βi為對應(yīng)的應(yīng)力參數(shù)。

在單元內(nèi)部假設(shè)的位移場

為三維八節(jié)點(diǎn)等參單元的形函數(shù)，在這里自然坐標(biāo)系ξ，η，ζ∈[-1，1]，而ξi、ηi、ζi是i節(jié)點(diǎn)在自然坐標(biāo)系中相應(yīng)節(jié)點(diǎn)的坐標(biāo)值，i代表某一層中三維八節(jié)點(diǎn)單元的節(jié)點(diǎn)號(hào)，j代表層合板單元的某一層。

引入記號(hào)

是每層子單元的位移矢量。

當(dāng)位移已得到滿足，赫林格-賴斯納（Hellinger-Reissner）廣義變分原理的泛函為

其中：Ωe為單元區(qū)域；S為材料柔度矩陣；b為體力；t為面力；D為微分操作符；Γt為應(yīng)力邊界區(qū)域。

把式（3）、式（5）代入式（6）可得

其中：H為單元柔度矩陣；G為杠桿矩陣；fe為等效節(jié)點(diǎn)載荷，分別為

故式（9）可簡化為

其中：單元?jiǎng)偠染仃?/p>

對式（14）泛函對qe變分，即利用駐值條件0可以得到位移的矩陣狀態(tài)控制方程從而可以利用一般有限元的分析方法求得未知節(jié)點(diǎn)的位移qe，進(jìn)一步由式（4）求得單元位移場，并且由式（3）和式（13）求得單元的應(yīng)力場

1.2 假設(shè)應(yīng)力模式

因?yàn)樾纬呻s交應(yīng)力元?jiǎng)偠染仃嚂r(shí)必須對柔度矩陣H求逆，所以雜交應(yīng)力元存在的充分必要條件是柔度矩陣可逆，即柔度矩陣是非奇異的，也就是要求假設(shè)的應(yīng)力模式線性無關(guān)。

下面根據(jù)等函數(shù)方法給出三維八節(jié)點(diǎn)正交異性復(fù)合材料單元的一種假設(shè)應(yīng)力模式的推導(dǎo)[5-6]。

對于一般的三維八節(jié)點(diǎn)單元，可假設(shè)位移場如下

在自然坐標(biāo)系中，依據(jù)位移-應(yīng)變關(guān)系，可得

將式（19）代入式（1）可以求得

將式（19）代入式（20），調(diào)整系數(shù)后可以得到

σξ=β1+β2ξ+β3η+β4ηζ+β5ξζ+β6ξη

即

因此該假設(shè)應(yīng)力模式為

2 實(shí)例分析

例題1 有一懸臂梁，梁長150 mm，寬2.5 mm，高5 mm；一端固定；自由端中心承受5 N的載荷；材料的楊氏模量E為70 GPa，泊松比為0[7]。其計(jì)算結(jié)果如表1所示。

從表1的計(jì)算結(jié)果可以看出本文的解與文獻(xiàn)[7]的解析解的誤差很小，明顯好于使用Abaqus有限元軟件計(jì)算的結(jié)果，并且該方法使用較少的網(wǎng)格劃分（2×12）就可以得到很精確的解。

表1 懸臂梁自由端的最大撓度Tab.1 The maximum deflection of the cantilever

例題2 考慮一三層板，且a=b，在x=0，a處固支，在y=0，b兩邊簡支，板的總厚度為h，h1=h3= 0.01h，h2=0.08h，第一層、第三層的材料參數(shù)相同，第二層材料參數(shù)為，板的上表面受均布?jí)簭?qiáng)q的作用。材料參數(shù)如下

計(jì)算結(jié)果如表2所示，其中I1、I2和I3分別為第一、第二和第三層所分的薄層數(shù)，從表1可以看出本文的解與文獻(xiàn)[8]的精確解比較吻合，并且顯然好于Abaqus采用的C3D8R單元類型有限元計(jì)算的結(jié)果。

表2 三層層合板的撓度和應(yīng)力Tab.2 Three laminated plates of deflection and stress

3 結(jié)語

通過修正的H-R混合變分原理，直接借助應(yīng)力-應(yīng)變關(guān)系，推導(dǎo)了新的應(yīng)力模式。在此基礎(chǔ)上，建立了復(fù)合材料雜交等參元的有限元列式，計(jì)算步驟簡單，數(shù)據(jù)累計(jì)誤差小，提高了計(jì)算精度，同時(shí)，我們發(fā)現(xiàn)該方法對網(wǎng)格的劃分不是很敏感，即可用較少的網(wǎng)格劃分就能得到比較精確的解，隨著網(wǎng)格的加密所計(jì)算的結(jié)果變化不大，因此該方法不僅可節(jié)約大量的計(jì)算時(shí)間，而且還能降低硬件的配置，從而節(jié)省了資源。

[1]卞學(xué)璜.雜交應(yīng)力有限元法的研究進(jìn)展[J].力學(xué)進(jìn)展，2001，31（3）：344-349.

[2] SZE K Y，GHALL A.Hybrid hexahedral element for solids，plates，shells and beams by selective scaling[J].International Journal of Numerical Methods in Engineering，1993，36：1519-1540.

[3]CHEUNG Y C，DI SHENGLIN.Analysis of laminated composite plate by hybrid stress isoperimetric element[J].International Journal of Solids Structures，1993，30（20）：2843-2857.

[4] SPILKER R L.Hybrid stress eight-node elements for thin and thick multilayer laminated plates[J].International Journal for Numerical Methods in Engineering，1982，18：801-828.

[5] FENG W，HOA S V，Huang Q.Classification of stress modes in assumed stress fields of hybrid finite elements[J].International Journal for Numerical Methods in Engineering，1997，40：4313-4339.

[6]張東鋒.壓電層合板的雜交元方法[D].成都：西南交通大學(xué)，2002.

[7]莊 茁，張 帆，岑 松，等.ABAQUS非線性有限元分析與實(shí)例[M].北京：科學(xué)出版社，2005.

[8]范家讓.強(qiáng)厚度疊層板殼的精確理論[M].北京：科學(xué)出版社，1996：130-351.

（責(zé)任編輯：楊媛媛）

Hybrid finite element method for laminated composite plate

QING Guang-hui，JIA Rui-sheng
（College of Aeronautical Engineering，CAUC，Tianjin 300300，China）

In this paper based on modified H-R mixed variational principle for composite materials，with the stressstrain relations directly，derivation of a new mode of stress，the hybrid and isoparametric finite element formulation for the laminated composite plate was established.Then the Mathematica was applied for the programming and calculating of a numerical example.Compared with the modeling analysis result using Abaqus software and the exact solution provided in relevant literatures concerning some mechanical quantities，the result obtained in this way is proved to be closer to their exact solutions and satisfactory precision can be obtained with less mesh.

composite materials；hybrid finite element；isoparametric element；Mathematica

O242.21

A 文章編號(hào)：1674－5590（2013）01－0082－03

2012-06-08；

2012-07-11

中國民航大學(xué)科研基金項(xiàng)目（2012kye07）

卿光輝（1968—），男，湖南新化人，教授，博士，研究方向?yàn)榻Y(jié)構(gòu)力學(xué).