復(fù)合材料層合板應(yīng)力分析的小波有限元法*

彭惠芬 孟廣偉 范森 周立明

(1.吉林大學(xué)機(jī)械科學(xué)與工程學(xué)院，吉林長(zhǎng)春130022;2.東北石油大學(xué)機(jī)械科學(xué)與工程學(xué)院，黑龍江大慶163318;3.東北石油大學(xué)石油工程學(xué)院，黑龍江大慶163318)

近十年來，隨著復(fù)合材料層合板在航空、航天、船舶、建筑等領(lǐng)域廣泛應(yīng)用，復(fù)合材料層合板有限元分析越來越受到工程界、學(xué)術(shù)界的廣泛重視.由于復(fù)合材料層合板具有拉－剪、彎－扭和拉剪－彎扭耦合效應(yīng)，并且在板邊緣附近應(yīng)力分布復(fù)雜、變化梯度大，給傳統(tǒng)有限元網(wǎng)格的劃分和高精度單元的建立帶來很大困難［1-4］.小波有限元是近年來發(fā)展起來的一種新型數(shù)值分析方法，已廣泛應(yīng)用于故障診斷、振動(dòng)分析、熱傳導(dǎo)、電磁場(chǎng)等各個(gè)工程領(lǐng)域.以小波函數(shù)或尺度函數(shù)替代傳統(tǒng)多項(xiàng)式插值，是優(yōu)于單元網(wǎng)格加密和階次升高的自適應(yīng)算法［5-6］.目前，小波有限元研究的熱點(diǎn)是如何構(gòu)造精度高、穩(wěn)定性好的小波單元以滿足各類工程需要，為此，許多學(xué)者開展了相關(guān)方面的研究，周又和等［7］利用Daubechies小波，構(gòu)造了小波梁?jiǎn)卧桶鍐卧?向家偉等［8］構(gòu)造了區(qū)間B樣條小波單元;Phoon等［9］構(gòu)造了一類用于求解大梯度問題的小波單元;然而，這些小波單元絕大部分適用于各向同性材料的靜動(dòng)力學(xué)分析.

文中基于層合板理論，利用樣條函數(shù)待定系數(shù)少、連續(xù)性強(qiáng)、逼近精度和計(jì)算效率高的優(yōu)點(diǎn)［10］，采用同尺度不同階數(shù)區(qū)間B樣條小波(BSWI)尺度函數(shù)張量積對(duì)層合板位移和撓度插值，構(gòu)造了BSWI二維C0型和C1型單元轉(zhuǎn)換矩陣，將小波系數(shù)轉(zhuǎn)化為節(jié)點(diǎn)物理空間自由度，并從虛功原理出發(fā)，推導(dǎo)了BSWI層合板單元的單元?jiǎng)偠确匠?

1 BSWI層合板單元轉(zhuǎn)換矩陣

基于經(jīng)典層合板理論，在單元的交界面上需同時(shí)滿足位移u、v，撓度w及其導(dǎo)數(shù)的兼容性和連續(xù)性.文中采用同尺度不同階數(shù)的BSWI尺度函數(shù)張量積插值，構(gòu)造層合板單元轉(zhuǎn)換矩陣，單元節(jié)點(diǎn)及自由度排列如圖1所示.

單元長(zhǎng)度分別為 lex、ley，單元節(jié)點(diǎn)數(shù)為 n×n(n=2j+1，j為BSWI尺度函數(shù)的尺度)，單元總自由度數(shù)為3n2+4n+4，單元節(jié)點(diǎn)自由度具體布置如下:編號(hào)為1、n、n2－n+1、n2的4 個(gè)角節(jié)點(diǎn)自由度為 ui、vi、wi、?wi/?x、?wi/?y、?2wi/?x?y;左、右邊上內(nèi)部節(jié)點(diǎn)自由度為:ui、vi、wi、?wi/?y;上、下邊上內(nèi)部節(jié)點(diǎn)自由度為:ui、vi、wi、?wi/?x;其他為內(nèi)部節(jié)點(diǎn)，自由度為:ui、vi、wi.其中，ui、vi、wi分別為第 i個(gè)節(jié)點(diǎn)沿坐標(biāo)軸x、y、z方向位移.

圖1 小波單元節(jié)點(diǎn)布置Fig.1 Node distribution of wavelet finite element

采用BSWI尺度函數(shù)二維張量積插值，在單元邊界節(jié)點(diǎn)上對(duì)未知場(chǎng)函數(shù) u0(ξ，η)、v0(ξ，η)進(jìn)行插值，構(gòu)造二維C0型單元轉(zhuǎn)換矩陣;同時(shí)對(duì)未知場(chǎng)函數(shù)w(ξ，η)及其導(dǎo)數(shù)進(jìn)行插值，構(gòu)造二維C1型單元轉(zhuǎn)換矩陣.未知場(chǎng)函數(shù)可表示為

式中:u0(ξ，η)，v0(ξ，η)表示層合板中面沿 x、y 方向位移;w(ξ，η)表示層合板撓度 ;ξ、η為單元局部坐標(biāo)，表示為:ξ=(x－x1)/lex、η =(y－y1)/ley，其中x1、y1為整體坐標(biāo)下單元起始坐標(biāo);Φ1、Φ2和 Φ3、Φ4分別是同尺度不同階數(shù)的一維區(qū)間B樣條尺度函數(shù)，即

αe、βe、γe為單元上待求的小波插值系數(shù)列向量，分別表示為

定義單元節(jié)點(diǎn)位移列陣:

將式(3)代入式(1)得到由節(jié)點(diǎn)位移表示的未知場(chǎng)函數(shù):

式中，C0型單元轉(zhuǎn)換矩陣為C1型單元轉(zhuǎn)換矩陣表示為其中:

式中，Φs，m(i)表示BSWI尺度函數(shù)Φs對(duì)m偏導(dǎo)數(shù)在相應(yīng)點(diǎn)i處的取值.其中

2 BSWI復(fù)合材料層合板單元構(gòu)造

根據(jù)層合板理論［11-16］，取其中面為參考面，忽略厚度方向的變形(εz=0)，根據(jù)小變形假設(shè)，層合板的位移－應(yīng)變關(guān)系可表示為

地產(chǎn)是指土地、建筑物及固著在土地、建筑物上不可分離的部分及其附帶的各種權(quán)益。地產(chǎn)可以分為一線、二線、三四線等；

式中:u0、v0和w分別為中面沿x、y和z軸方向的位移分量;z為距中面位移;εx、εy和 γxy分別為沿 x、y和z軸方向的應(yīng)變分量.

由于層合板各層剛度及材料主方向的不同，第k層應(yīng)力－應(yīng)變關(guān)系可表示為

式中，

其中:θ表示從 x軸轉(zhuǎn)向主軸的角度;C11=E1和E2分別為材料在1、2主方向上的彈性模量，υ12和 υ21為泊松比.

當(dāng)層合板在外力作用下處于平衡時(shí)，且產(chǎn)生符合約束的微小虛位移，則第k層的虛功方程為

式中:虛應(yīng)變列陣 δε =［δεxδεyδγxy］T，應(yīng)力列陣 σk=［σxσyxy］T.

由于層合板應(yīng)力的不連續(xù)分布，對(duì)所鋪層求和后，可得層合板總虛功為

定義單元等效節(jié)點(diǎn)力列陣:

式中:

其中，面力px和py分別為層合板x、y方向軸向拉力，q為垂直于層合板面的壓力，pj為集中載荷.

外力在虛位移上所作的虛功為

將式(7)－(10)聯(lián)立，并將局部坐標(biāo)下單元求解域 Ωs={ξ，η|ξ，η∈［0，1］}映射為整體坐標(biāo)下標(biāo)準(zhǔn)求解域 Ωe，由虛功原理可得層合板 BSWI剛度方程:

3 算例分析

為驗(yàn)證文中構(gòu)造的BSWI層合板單元的正確性和有效性，圖2給出了16層等厚層合板結(jié)構(gòu)示意圖.層合板由上、下8層材料主方向與坐標(biāo)軸夾角分別為－45°和 45°層合板組成，簡(jiǎn)記為 －45°8/45°8，單層板板厚t=0.125mm，長(zhǎng)度a=1m，寬度b=0.5m，E1=181.00GPa，E2=10.30 GPa，v21=0.28，G12=7.17GPa，在其兩側(cè)承受軸向載荷px=5×103N/m.求:

1)對(duì)角線OA上各點(diǎn)撓度隨距離的變化;

2)層合板應(yīng)力、應(yīng)變沿板厚的分布規(guī)律.

圖2 層合板結(jié)構(gòu)示意圖Fig.2 Structure diagram of laminated composite plates

采用階數(shù)m=2，尺度 j=3和階數(shù) m=4，尺度j=3的BSWI尺度函數(shù)(分別簡(jiǎn)記為 BSWI23和BSWI43)的張量積對(duì)板中面位移及撓度插值，構(gòu)造BSWI層合板單元，單元數(shù)量為1個(gè)，總自由度數(shù)為283個(gè).

表1中給出了層合板承受單向軸向拉伸時(shí)，分別采用BSWI小波有限元法和ANSYS數(shù)值分析法求解層合板1層、8(+45°)層和16層應(yīng)力值，并將BSWI法與解析解進(jìn)行誤差比較.由表中結(jié)果可見:BSWI法計(jì)算精度明顯高于ANSYS數(shù)值分析法，且BSWI法計(jì)算層合板最大應(yīng)力的相對(duì)誤差不超過5.03%，說明文中所構(gòu)造的BSWI層合板單元是正確可行的，同時(shí)也說明本文所構(gòu)造的BSWI層合板單元繼承了樣條函數(shù)逼近精度高，連續(xù)性強(qiáng)的優(yōu)點(diǎn)，可用較少單元和自由度數(shù)獲得較高的計(jì)算精度.

表1 BSWI法求解各層應(yīng)力值與解析解比較Table 1 Comparison of analytical results of stress with those obtained by BSWI method

圖3為對(duì)角線OA上各點(diǎn)撓度隨距O點(diǎn)距離的變化.由圖可見，由于具有拉剪耦合效應(yīng)，在均勻軸向拉伸載荷作用下，非對(duì)稱層合板發(fā)生了彎曲，撓度隨板面上距O點(diǎn)距離非線性變化.其中，O點(diǎn)撓度為0，4個(gè)角點(diǎn)撓度最大為5.42 mm.圖4為層合板各層應(yīng)力、應(yīng)變沿板厚的分布圖，由圖可見，層合板應(yīng)變沿厚度連續(xù)分布，而各層應(yīng)力不連續(xù)分布，這主要由于層合板各層剛度不同造成的.

圖3 對(duì)角線OA上點(diǎn)的撓度與距離的變化關(guān)系Fig.3 Relationship between deflection of nodes on diagonal OA and the distance to the center of plate

圖4 層合板應(yīng)變、應(yīng)力沿板厚分布圖Fig.4 Distributions of stress and strain of laminated composite plates in thickness direction

4 結(jié)論

文中基于經(jīng)典層合板理論，采用同尺度不同階數(shù)區(qū)間B樣條小波尺度函數(shù)的二維張量積插值，構(gòu)造二維C0型和C1型層合板單元轉(zhuǎn)換矩陣，從而將無明確物理意義的小波系數(shù)轉(zhuǎn)換為節(jié)點(diǎn)位移坐標(biāo)，方便了相鄰單元的連接和邊界條件的處理，并利用虛位移原理得到了BSWI層合板單元的單元?jiǎng)偠确匠?數(shù)值算例表明，文中構(gòu)造的BSWI層合板單元在薄板的靜態(tài)分析中，具有求解精度高、所用單元和自由度數(shù)量少的優(yōu)點(diǎn)，為進(jìn)行較復(fù)雜層合板結(jié)構(gòu)分析提供了高精度單元.

［1］Lee Sang-Ho，Yoon Young-Cheol.Numerical predication of crack propagation by an enhanced element-free Galerkin method ［J］.Nuclear Engineering and Design，2004(3):257-271.

［2］Metin Aydogdu.A new shear deformation theory for laminated composite plates ［J］.Composite Structures，2009(89):94-101.

［3］Ferreira A J M.A formulation of the multiquadric radial basis function method for the analysis of laminated composite plates［J］.Composite Structures，2003(59):385-392.

［4］王勖成.有限單元法［M］.北京:清華大學(xué)出版社，2002.

［5］何正嘉，陳雪峰.小波有限元理論研究與工程應(yīng)用的進(jìn)展 ［J］.機(jī)械工程學(xué)報(bào)，2005，41(3):1-11.He Zheng-jia，Chen Xue-feng.Advanced in theory study and engineering application of wavelet finite element［J］.Chinese Journal of Mechanical Engineering，2005，41(3):1-11.

［6］Kim Y Y，Jang G W.Hat interpolation wavelet-based multi-scale Galerkin method for thin-walled box beam analysis［J］.International Journal for Numerical Methods in Engineering，2002，53:1575-1592.

［7］周又和，王記增，鄭曉靜.小波伽遼金有限元法在梁板結(jié)構(gòu)中的應(yīng)用 ［J］.應(yīng)用數(shù)學(xué)和力學(xué)，1998，19(8):697-706.Zhou You-he，Wang Ji-zeng，Zheng Xiao-jing.Applications of wavelet Galerkin FEM to beam and plate structures［J］.Applied Mathematics and Mechanics，1998，19(8):697-706.

［8］向家偉，陳雪峰，何正嘉.區(qū)間B樣條小波平面彈性及Mindlin板單元構(gòu)造研究［J］.計(jì)算力學(xué)學(xué)報(bào)，2007，24(6):869-875.Xiang Jia-wei，Chen Xue-feng，He Zheng-jia.A study of the construction of wavelet-based plane elastomechanics and mindlin plate elements using B-spline wavelet on the interval［J］.Chinese Journal of Computational Mechanics，2007，24(6):869-875.

［9］Phoon K K，Huang S P，Quek S T.Implementation of Karhunen-Loeve expansion for simulation using a wavelet-Galerkin scheme［J］.Probabilistic Engineering Mechanics，2002(3):293-303.

［10］安寧剛，陳雪峰，曹巨江.樣條小波有限元的算法研究及其應(yīng)用［J］.陜西科技大學(xué)學(xué)報(bào)，2003，21(3):56-59.An Ning-gang，Chen Xue-feng，Cao Ju-jiang.The study and application of apline wavelets element methods［J］.Journal of Shaanxi University of Science ＆ Techonlogy，2003，21(3):56-59.

［11］舒小平，江永真，史林興.精確的復(fù)合材料層合板有限元模型［J］.淮海工學(xué)院學(xué)報(bào)，2003，12(1):8-11.Shu Xiao-ping，Jiang Yong-zhen，Shi Lin-xing.An accurate finite element of laminated composite plates［J］.Journal of Huaihai Institute of Technology，2003，12(1):8-11.

［12］鐘軼峰，余文斌.用漸近變分法對(duì)復(fù)合材料層合板簡(jiǎn)化數(shù)值建模及仿真［J］.重慶大學(xué)學(xué)報(bào)，2011，34(2):130-135.Zhong Yi-feng，Yu Wen-bin.Simplified numerical modeling and simulation of composite laminated plates by the variational asymptotic method ［J］.Journal of Chongqing University，2011，34(2):130-135.

［13］Shu X.A refined theor y of laminated shells with higherorder transverse shear deformation ［J］.Int J solids Struct，1997，34(6):673-683.

［14］Ghafoori E，Asghari M.Dynamic analysis of laminated composite plates traversed by a moving mass based on a first-order theory ［J］.Composite Structures，2010(92):1865-1876.

［15］沈觀林，胡更開.復(fù)合材料力學(xué)［M］.北京:清華大學(xué)出版社，2006.

［16］張萬平，張俊乾，吳堅(jiān).金屬基復(fù)合材料纖維斷裂脫落后的應(yīng)力分析［J］.重慶大學(xué)學(xué)報(bào)，2004，27(9):119-123.Zhang Wan-ping，Zhang Jun-qian，Wu Jian.Analysis for a broken and debonging fiber in ductile matrix composites［J］.Journal of Chongqing University，2004，27(9):119-123.