复合材料层合结构应力分析

复合材料层合结构应力分析

随着工业和技术的进步，一体化板的外壳结构在航空航天等领域得到了越来越多的应用。关于复合材料层合板壳的有限元分析一直是一个令人关注的问题。为了克服传统的退化壳元难以处理大的转动增量的问题,最近几年人们把研究目光转向了固体壳元,但随之而来的是固体壳元的厚度自锁问题。现有的解决方法有采用平面应力假设的本构方程、加强假定应变列式和杂交应力元列式等,但这些方法都存在着各种各样的问题。采用平面应力假设的本构方程在处理厚度变化较大的层合结构时会导致厚度方向正应力的明显不连续分布,加强假定应变列式同样不能得到连续的厚度方向的应力分布,且它要在单元一级上消除内部参数,降低了计算效率,杂交应力元列式虽可得到连续的厚度方向的应力分布,但有时会影响位移的计算精度。本文中则引入广义应力的概念,建议了一种新的改进本构阵的方法,克服了厚度自锁问题,尤其适用于复合材料层合板壳结构的数值分析。正交化方法是近年来多变量有限元研究领域取得的一个主要成果,它能明显地提高低阶单元的计算效率。在高阶单元的研究方面,S.W.Lee等建议了一种杂交稳定列式,将应变插值函数分为低阶和高阶两部分,与低阶项对应的刚度阵和相应的减缩积分单元等价,与高阶项对应的部分用来克服单元的零能模式,但其建议的列式中对应高阶项的部分还要采用高阶的数值积分,且要求一满阵的逆。本文列式中借鉴了低阶杂交元列式中的正交化方法,消除了低阶和高阶项的交叉项,并提出一个精化措施,只用低阶积分即可求出高阶稳定阵,且避免了矩阵求逆运算,显著提高了计算效率。1根据ddlagrangaan的位移法，有限列方程1.1增量应变变分的计算九节点固体壳元(每个节点六个自由度)的坐标、位移采用下列的标准形式:x(ξ,η,ζ)=9∑i=1Νi(ξ,η)(xi+ζvi)=Ν(ξ,η)ˉx+ζΝ(ξ,η)ˉv=x0(ξ,η)+ζxn(ξ,η)(1)U(ξ,η,ζ)=9∑i=1Νi(ξ,η)[Ui+ζ(v´i-vi)]=U0(ξ,η)+ζUn(ξ,η)(2)在传统的退化壳元中,假定变形前后壳体的厚度保持不变,故v′i-vi(v′i、vi分别为变形前、后的中面法线)由两个转动矢量来表示,目前的列式中抛弃了这一假设,因此v′i-vi为一有三个分量的矢量。ΔU(ξ,η,ζ)=9∑i=1Νi(ξ,η)(ΔUi+ζΔVi)=Ν(ξ,η)ΔˉU+ζΝ(ξ,η)ΔˉV=ΔU0(ξ,η)+ζΔUn(ξ,η)(3)其中:xi、Ui、ΔUi分别是节点坐标、节点位移和增量节点位移矢量,Ni是二维拉格朗日插值函数。ˉx={x1⋯x9}={X1+U1⋯X9+U9},ˉV={V1⋯V9},ΔˉU={ΔU1⋯ΔU9},ΔˉV={ΔV1⋯ΔV9},Ν=[Ν1,⋯,Ν9]由应变定义,t时刻的Almansi逆变应变和增量Almansi逆变应变为εi=εii=xT′iU′i-UT′iU′i/2,2εij=xΤ´iU′j+xΤ´jU′i-(UΤ´iU′j+UΤ´jU′i)/2(4)Δεi=Δεii=XΤ´iΔU′i+ΔUΤ´iΔU′i/2,2Δεij=XΤ´iΔU′j+XΤ´jΔU′i-(ΔUΤ´iΔU′j+ΔUΤ´jΔU′i)/2(5)其中:i,j=ξ,η,ζ,且i≠j。ε=={εx′εy′2εx′y′}=Τ={εξεη2εξη},εt={2εz′x′2εz′y′}=Τt{2εζξ2εζη},εz=Τzεζ(6)Δε=={Δεx′Δεy′2Δεx′y′}=Τ={ΔεξΔεη2Δεξη},Δεt={2Δεz′x′2Δεz′y′}=Τt{2Δεζξ2Δεζη},Δεz=ΤzΔεζ(7)在增量应变中忽略关于ζ的二次项,得Δε==Δεm+ζΔεb=(BLm+BΝm/2)ΔˉU+ζ(BLb+BΝb/2)ΔˉU,Δεz=(BLz+BΝz/2)ΔˉU,Δεt=Δεo+ζΔεn=(BLo+BΝo/2)ΔˉU+ζ(BLn+BΝn/2)ΔˉU(8)显然增量应变的变分可写成δΔελ=(BLλ+BΝλ)δΔˉUδΔεz=(BLz+BΝz)δΔˉU(9)其中:λ=m,b,o,n。当前时刻的Cauchy应力和增量Kirchhoff应力可写成{σ=σz}=[C=C×CΤ×Cz]{ε=εz},σt=Ctεt,{ΔσΔσz}=[C=C×CΤ×Cz]{Δε=Δεz},Δσt=CtΔεt(10)每个单元的最小势能原理可写为ΠeΡ=12n∑l=11∫-11∫-11∫-1[{ε=+Δε=εz+Δεz}Τ[C=C×CΤ×Cz]⋅{ε=+Δε=εz+Δεz}+(εt+Δεt)ΤCt(εt+Δεt)]⋅Jh1hdξdηdζl-(ˉU+ΔˉU)ΤF(11)其中:n是复合材料的总铺层数,h是所有层的总厚度,hl是各层的厚度,J是当前构形与等参构形间的雅可比转换矩阵,F是外部载荷形成的等效节点力。为提高计算效率,近似取ˉJ=J|ζl=0,同时代入式(10),有ΠΡe=12〈(ε⊥+Δε⊥)ΤC˜⊥(ε⊥+Δε⊥)+(εΤ+ΔεΤ)ΤCΤ(εΤ+ΔεΤ)〉-(U¯+ΔU¯)ΤF(12)其中:ε⊥={εmεzεb},Δε⊥={ΔεmΔεzΔεb},C˜⊥=12∑l=1n∫-11[C=C×ζlC=C×ΤCzζlC×ΤζlC=ζlC×ζl2Cz]hlhdζl,εt={εoεn},Δεt={ΔεoΔεn},CΤ=12∑l=1n∫-11[CtζlCtζlCtζl2Ct]hlhdζl,〈o〉=2∫-11∫-11(o)J¯dξdη此外,我们定义当前时刻的广义应力如下(增量广义应力类似):Ν⊥={ΝΤΜ}=12∑l=1n∫-11{σ=σzζlσ=}hlhdζl=12∫-11{C=(εm+ζlεb)+C×εzC×Τ(εm+ζlεb)+CzεzζlC=(εm+ζlεb)+ζlC×εz}hlhdζl=C˜⊥ε⊥QΤ={QoQn}=12∑l=1n∫-11{σtζlσt}hlhdζl=12∫-11{Ct(εo+ζlεn)ζlCt(εo+ζlεn)}hlhdζl=CΤεΤ将广义应力和增量广义应力代入变分泛函(12),有δΠΡe=〈δ(Δε⊥)Τ(Ν⊥+C˜⊥Δε⊥)+δ(ΔεΤ)Τ(ΝΤ+CΤΔεΤ)〉-δ(ΔU¯)ΤF(13)经线性化处理,有δΠΡe=δ(ΔU¯)Τ[(ΚL+ΚΝ)ΔU¯+R-F](14)其中单元刚度矩阵为ΚL=〈B⊥LΤC˜⊥B⊥L〉,ΚΝΔU¯=〈B⊥ΝΤΝ⊥+BΤΝΤQΤ〉(15)等效节点力为R=〈B⊥LΤΝ⊥+BΤLΤQΤ〉(16)其中:B⊥L[BmLBzLBbL],B⊥L=[BmΝBzΝBbΝ],BΤL=[BoLBnL],BΤΝ=[BoΝBnΝ]1.2改进的本构阵将式(10)改写为{σ=εz}=[AB-BΤD]{ε=σz}=[S=-1-S=-1S×-S×ΤS=-1Sz-S×ΤS=-1S×]{ε=σz}(17)为了克服厚度自锁,定义εz=12∑l=1n∫-11([-BΤD]{ε=σz})hlhdζl=12∑l=1n∫-11([-BΤD]{εm+ζlεbσz})hlhdζl=[-B0ΤD0-B1Τ]{εmσzεb}(18){ΝΜ}=12∑l=1n∫-11{σ=ζlσ=}hlhdζl=12∑l=1n∫-11[ABζlAζlB]{εm+ζlεbσz}hlhdζl=[A0B0A1A1B1A2]{εmσzεb}(19)其中:[A0A1A2B0B1D0]=12∑l=1n∫-11[AζlAζl2ABζlBD]hlhdζl由式(10),(18),(19)得改进的本构阵为C⊥=[A0+B0B0Τ/D0B0/D0A1+B0B1Τ/D0B0Τ/D01/D0B0Τ/D0A1+B1B0Τ/D0B1/D0A2+B1B1Τ/D0](20)2单刚单干单阶项与减缩积分元等价的计算为了表述的简洁性,引入记号ue001φ,ue001φ代表⊥和T(⊥和T所代表的意义与1.1中所述相同)。Hellinger-Reissner变分原理可写成ΠΗRe=∑φ=⊥,Τ〈-(σφΤ+ΔσφΤ)Cφ-1(σφ+Δσφ)/2+(σφΤ+ΔσφΤ)(εφ+Δεφ)〉-(U¯+ΔU¯)ΤFδΠΗRe=∑φ=⊥,Τ〈δ(ΔσφΤ)[Cφ-1(σφ+Δσφ)-(εφ+Δεφ)]+(σφΤ+ΔσφΤ)δ(Δεφ)〉-δ(ΔU¯)ΤF(21)广义应力分为低、高阶两部分,即:Δσue001φ=Pue001φLΔβue001φL+Pue001φHΔβue001φH,其中应力插值函数满足下列正交性条件:〈PTue001φLC-1ue001φPue001φH〉=0,引入记号φ,φ代表L和H(分别代表低、高阶项,以下各式中同),则应力插值函数可简化为:Δσφ=∑φΡφφΔβφφ,将应力插值函数和式(9)、(20)代入式(21)有δΠΗRe=∑φφ〈-δ(ΔβφφΤ)[ΗφφΔβφφΤ-GφφΔU¯-Ζφφ]+δ(ΔU¯Τ)GφφΔβφφΤ〉+δ(ΔU¯)Τ(ΚσΔU¯+Q-F)(22)其中:Ηφφ=〈ΡφφΤCφ-1Ρφφ〉,Gφφ=〈ΡφφΤBφ〉,Ζφφ=〈ΡφφΤ(εφ-Cφ-1σφ)〉(23)由式(22),得独立假定的增量应力参数为Δβφφ=Ηφφ-1(GφφΔU¯+Ζφφ)(24)将式(24)代入式(22),得δΠΗRe=δ(ΔU¯Τ){(∑φ,φGφφΤΗφφ-1Gφφ+Κσ)ΔU¯-F+Q+∑φ,φGφφΤΗφφ-1Ζφφ}(25)为了提高计算效率,低阶项部分选择如下的插值函数,使本杂交元列式的单刚低阶项部分与减缩积分元等价Ρ⊥L=[Ρ1Ι7Ρ2Ι7Ρ3Ι7Ρ4Ι7],ΡΤL=[Ρ1Ι4Ρ2Ι4Ρ3Ι4Ρ4Ι4]其中:Ρ1=(1-ξ3)(1-η3)/4,Ρ2=(1+ξ3)(1-η3)/4,Ρ3=(1+ξ3)(1+η3)/4,Ρ4=(1-ξ3)(1+η3)/4则式(23)可简化为GφL=〈ΡφLΤBφ〉L=[J¯1Bφ1⋯J¯4Bφ4]‚ΖφL=[J¯1(εφ-Cφ-1σφ)1⋯J¯4(εφ-Cφ-1σφ)4]‚ΗφL=〈ΡφLΤCφ-1ΡφL〉L=diag[J¯1Cφ1-1,⋯,J¯4Cφ4-1]其中∶〈o〉L=∑i=14J¯i(o)i‚(o)i为矩阵或变量‘o’在第i个2×2高斯积分点的值。所以有GφLΤΗφL-1GφL=〈BφΤCφBφ〉L,GφLΤΗφL-1ΖφL=〈BφΤ(Cφεφ-σφ)〉L故式(25)可变为δΠΗRe=δ(ΔU¯Τ){[∑φ〈BφΤCφBφ〉L+GφΗΤΗφΗ-1GφΗΤ)+Κσ]ΔU¯-F+Q+∑φ[〈BφΤ(Cφεφ-σφ)〉L+GφΗΤΗφΗ-1ΖφΗ]}(26)为克服减缩积分元的零能模式,选取下列高阶应力插值函数:Ρ⊥Η=1J[Τ=-1Ρ⊥03×201×201×203×2Τ=-1Ρ⊥],ΡΤΗ=1J[Τt-1Ρ⊥02×102×1Τt-1ΡΤ],Ρ⊥=[Ρξ00Ρη00],ΡΤ=[ΡξΡη],Ρξ=ξ(η2-1/3)‚Ρη=η(ξ2-1/3)则G⊥Η=〈Ρ⊥ΗΤB⊥〉L=[GmξGmηGbξGbη]=[x¯Μmξx¯Μmηx¯Μbξx¯Μbη]‚Η⊥Η=〈Ρ⊥ΗΤC⊥-1Ρ⊥Η〉=〈1J02[Ρ⊥03×201×201×203×2Ρ⊥]ΤS⊥[Ρ⊥03×201×201×203×2Ρ⊥]〉GΤΗ=〈ΡΤΗΤBΤ〉=[GtoGtn]=x¯Τ[ΜtoΜtn],ΗΤΗ=〈ΡΤΗΤCΤ-1ΡΤΗ〉=〈1J02[Ρ⊥02×102×1ΡΤ]ΤSΤ[ΡΤ02×102×1ΡΤ]〉其中:S⊥=diag{Τ=-1,0,Τ=-1}⋅C⊥-1⋅diag{Τ=-Τ,0,Τ=-Τ},SΤ=diag{Τt-1,Τt-1}⋅CΤ-1⋅diag{Τt-Τ,Τt-Τ}其中:Mmξ、Mmη、Mbξ、Mbη、Mto、Mtn可利用Maple(WaterlooMapleInc.)等软件推得显式。近似取:J¯≈J0=J¯|ξ=η=0,S⊥≈S¯⊥=S¯⊥|ξ=η=0‚SΤ≈S¯Τ=S¯Τ|ξ=η=0,则H⊥H,HTH均可由2×2的高斯积分求出。代入式Q=∑φ〈BφΤσφ〉L‚ΖφΗ=0,式(26)变为δΠΗRe=δ(ΔU¯Τ){[∑φ〈BφΤCφBφ〉L+GφΗΤΗφΗ-1GφΗ+Κσ]ΔU¯-F+∑φ〈BφΤCφεφ〉L}(27)3计算值的示例3.1材料及结构几何四边固支圆柱壳受如图1所示均匀分布径向力作用,其几何尺寸为:R=2540mm,L=508mm,θ=0.1°,壳体总厚度为h=2.54mm,考虑[0°/90°]铺层,所使用的材料常数为:E1=25×106N/mm2,E2=E3=106N/mm2,G12=G31=0.5×106N/mm2,G23=0.2×106N/mm2,ν12=ν23=ν13=0.25。由于结构和载荷的对称性,仅将四分之一壳体离散为2×2网格。图2反映了使用十个加载步时,变形随载荷变化的情况,并与文献中的结果进行了比较。3.2流固体的结构和载荷一端固支圆柱壳,如图3所示自由端曲边中点受一集中力作用,其几何尺寸为:R=101.6mm,h=3mm,L=304.8mm,材料为横向各向同性,弹性常数为:EL=2068.5N/mm2,ET=517.125N/mm2,νLT=νTL=0.3,GLT=795.6N/mm2,考虑90°/0°/90°和0°/90°/0°两种铺层。壳体所受外力为Pmax=36N,由于结构和载荷的对称性,仅将二分之一半圆柱壳体离散为8×8网格。图4反映了A点的挠度随载