基于20节点辛元的复合材料层合板应力分析*

童 瑶1,姚玉喆2

(1.上海飞机设计研究院, 上海 201206;2.中国民航大学 航空工程学院, 天津 300300)

摘要: 通常情况下,常规位移有限元法获得的应力结果比位移精度低一阶次,且面外应力难以满足连续性要求联合最小势能原理和H-R变分原理,构造出包含位移和3个面外应力两类变量的20节点六面体辛元由于两类变量采用高阶插值函数近似,无需引入单元内部的非协调位移项,因此相关理论的推导过程非常简单与Hamilton部分混合元不同,该辛元涉及的变量沿3个坐标方向均做离散处理,不受单元厚度和结构几何形状的限制数值实例表明20节点辛元的数值结果收敛稳定在粗糙网格的情况下,与20节点位移元相比,该文单元的面外应力更接近精确解

关 键 词: 复合材料;面外应力;H-R变分原理;20节点辛元;位移法;混合法

引 言

复合材料层合板壳结构的面外应力对于层间分层的起始和扩展具有重要影响为了模拟面外应力沿厚度方向的非线性分布,前人提出了许多解析方法和数值方法由于理论公式比较简单、计算资源消耗低且适应任何复杂的工程结构,位移有限元法具有吸引力但是,若采用较粗的网格,常规的位移法难以给出复合材料结构面外应力的理想结果一方面,因为位移法的应力结果是基于本构关系通过假设的位移插值函数的微分得到的,应力结果的精度比位移结果的精度低一阶次另一方面,由于位移法只有位移是独立假设的,当板厚变得很薄时,会出现剪切锁定现象,即使应用高阶层合板理论[1-2]也很难得到比较准确的沿厚度方向的面外剪切应力特别指出:如果不对相邻单元于同一节点上的应力结果进行磨平处理,位移法的结果很难满足复合材料层合结构面外应力的连续性要求

假设应力和位移两类变量的杂交应力元法[3-5]的应力结果是根据光滑的应力模式得到的这是杂交应力元法的应力结果比较准确的主要原因之一专著[3]中根据状态空间理论建立的部分杂交应力元法完全满足层间面外应力的连续性要求有关位移法和杂交应力元法求解复合材料层合板/壳应力问题的调查分析可以参考文献[1-6]

将空间坐标z类比为时间度量,H-R变分原理的Reissner函数等价于Hamilton原理[7]这正是半解析的Hamilton部分混合元[8-10]的理论基础Hamilton部分混合元包含了3个位移分量和全部面外应力分量,均采用相同的插值多项式该混合元的主要优点包括:1) 满足面外应力的层间连续性要求和面外应力的边界条件,同时满足面内应力和面外应变的层间不连续性要求; 2) 可同时获得高精度的位移结果和面外应力结果但是,就目前的发展情况来说,这类模型要求网格模型中同一子层内所有单元的厚度相同因此,不适合厚度不均匀或几何形式比较复杂的复合材料层合结构另一方面,对应Hamilton部分混合元的状态空间模型是线性微分方程系统,求解位移和应力之前,首先需要进行比较耗时的矩阵指数运算,不适应大规模的有限元分析

当然,包含位移变量和应力变量的经典混合有限元法[11-12]为同时引入应力边界条件和位移边界条件提供了方便,并可同时获得精度比较高的位移和应力数值结果因为混合有限元法中位移变量和应力变量的插值多项式相同,所以其位移和应力解的精度是一致的这也是该类方法的应力结果比较准确的原因之一一般情况下,经典混合有限元法的线性系统是非正定的(因为主对角线上存在零元素)因此,若不考虑稳定混合元技术方案[11-12],数值结果是不稳定的附加稳定技术方案混合元的数学模型比较复杂,不易被工程师们所接受这可能是经典混合法未能得到充分发展和广泛应用的原因之一另一方面,通过全混合元(包括位移、面外应力和面内应力的所有分量)对应的有限元线性方程系统求得的全部应力分量结果与位移结果相同,于材料界面间均是连续的然而,由于复合材料层合板壳结构各层材料不同,面内法向应力分量于层间通常是不连续的因此,全混合法不适应复合材料层合板壳结构的应力分析

为了更合理、更准确地分析复合材料层合板的层间面外剪切应力,Liao和Tsai[13]提出了部分混合有限元法,其中假定了全部位移分量和两个面外剪切应力分量,因而降低了混合有限元线性系统的自由度,提高了计算效率但是,这种模型忽略了面外法向应力于层间的连续性要求

最近,对文献[14]基于修正的H-R变分原理和势能原理建立了包含3个位移和全部面外应力分量的8节点非协调六面体辛元,其有限元线性系统对称保辛,并且系数矩阵主对角线上不存在零元素与经典混合法和Hamilton部分混合元法相同,该辛元满足层合结构层间面外应力的连续性要求,通过其有限元线性系统求得的面外应力结果于节点处自然连续与Hamilton部分混合元不同的是:该辛元所涉及的变量在3个坐标方向均被离散处理,对单元的厚度和几何形状没有限制文献[15]从修正的广义H-R 变分原理出发,建立了求解复合材料层合板稳态温度场静力学问题的非协调8节点辛元,是文献[14]中辛元的扩展和应用这两类辛元的共同缺点是:由于单元内非协调位移项的引入,需要对每个单元非协调部分的“刚度矩阵”求逆后方可消去内部节点的位移参数,该操作过程消耗大量的时间

本文的主要目的是建立形式比较简单、效率较高的20节点六面体等参辛元,并用其分析复合材料层合板的应力

1 基 本 理 论

1.1 最小势能原理

若不考虑体积力,最小势能原理可表示为

u)TC(

(1)

式中,V为弹性体的体积,S为结构的表面积,C表示材料的刚度矩阵,u=[u v w]T是关于3个方向(x1,x2,x3)的位移向量,为表面Sσ上的外载荷,是阶数为6×3的微分算子矩阵

1.2 20节点位移单元的刚度矩阵

三维六面体位移单元的位移变量u可表示为

u=Nqe,

(2)

式中,qe=[u1e u2e u3e]T表示单元节点的位移向量;N为对角矩阵,diag N=[Ne Ne Ne]T, Ne=[N1 N2N20]为等参形式的20节点插值多项式

将式(2)代入式(1),考虑位移向量qe的变分,则单元的刚度方程为

Dqqqe=fq,

(3)

其中N)TC(

1.3 修正的H-R变分原理

假设边界的位移事先满足,则修正的H-R变分原理[7-10,14](忽略体积力)为

(4)

其中LMHR是Reissner能密度函数:

2u)TΦ222u

(5)

式中,σo=[σxz σyz σzz]T是面外应力向量;123是微分算子矩阵;Φ11Φ21Φ22是材料参数矩阵

面内应力σi=[σxx σyy σxy]T可由位移和面外应力表示:

σi=φ21σo+φ22(2u)

(6)

1.4 20节点辛元

与文献[14-15]相同,面外应力可表示为

σo=Npe,

(7)

其中,pe=[σ13e σ23e σ33e]T为单元节点的面外应力向量

将式(2)和(7)代入式(4)可得修正的H-R变分原理的离散形式:

(8)

其中

1N+Φ212N+3N)dV,

2N)TΦ22(2N)dV,

fq的积分式与式(3)中fq的积分式相同

peqe视为相互独立的变量,由δΠMHR(pe,qe)=0可得

(9)

这里必需说明:式(8)和(9)中的Kqq的主对角线上含有零元素(因为微分算子矩阵2的主对角线上有零元素)实践表明:与经典的混合法一样,如果简单地考虑式(9)取和后的控制方程来求解具体问题,位移和应力结果并不随着网格的增加稳定地收敛

为了消去式(9)中主对角线上的零元素,将其与式(3)合并,则有

(10a)

其中

交换上式的行,可得

(10b)

显然,式(10b)与文献[14]中非协调辛元在形式上完全相同,满足辛空间Hamilton矩阵的定义因此,这里称式(10b)为20节点等参辛元

2 数 值 算 例

算例1 考虑三层的矩形层合梁[0°/90°/0°]层合板总厚度为H=0.1,上中下三层厚度均为H/3,长a=4H各层的材料参数[16]为:E11=25E22=25E33G12=G13=0.5E33,G23=0.2E33,μ12=μ13=μ23=0.25该梁的两端简支约束(当x1=0和x1=a时,σ11=u2=u3=0)上表面(x3=0.1处)垂直方向施加了单向正弦载荷q=q0sin(πx1/a)

因为几何模型、位移边界条件和载荷均是对称的,因此考虑1/4对称模型:x1方向划分为9个单元;x2方向划分为1个单元;每层厚度方向(x3)划分为2个单元,共54个单元20节点位移元(CDE20)和20节点辛元(CSE20)采用相同的有限元模型

图1和图2分别为面外应力沿厚度方向的分布图1~8图例中的CDE20代表20节点位移元,其应力结果是根据材料的本构关系计算Gauss点的应力, 然后用外推方法得到节点上的应力, 最后对相邻单元同一个节点应力值取平均值; CSE20代表本文20节点辛元, 面外应力结果和位移结果通过式(10)求得图1~3中EXACT(精确解)取自文献[16]

图1~8中各应力分量的结果采用了如下的无量纲公式:

(11a)

(11b)

式中,q0表示载荷,S表示长厚比算例1中的S=a/H=4

在相同的有限元网格条件下,图1表明本文辛元的面外剪切应力比位移元的面外剪切应力更贴近精确解事实上,CSE20的结果与精确解的最大误差小于1%图2、3表明辛元和位移元的结果差别甚微不难看出,面外法向应力和面内法向应力沿厚度方向的非线性分布特性不显著

算例2 考虑三层的矩形层合板[0°/90°/0°]层合板总厚度为H=0.1,上中下三层厚度均为H/3,长a=10H,宽b=a各层材料参数[17]为:E11=E22=4×104E33=5×105,G13=G23=6×104,G12=1.6×104μ12=μ13=μ23=0.25层合板为四边简支约束:当x1=0和x1=a时,σ11=u2=u3=0;当x2=0和x2=b时,σ22=u1=u3=0上表面(x3=0.1处)垂直方向施加均布载荷q0本例中用于无量纲化的S=a/H=10

考虑1/4对称的有限元模型:x1x2方向各6个单元;每层厚度方向(x3)2个单元图4~10中EXACT(精确解)采用了文献[17]中的方法求得

沿厚度的分布 图沿厚度的分布
Fig.1 Distribution of along Fig.2 Distribution of along
the thickness the thickness

沿厚度的分布 图沿厚度的分布
Fig.3 Distribution of along Fig.4 Distribution of
the thickness the thickness

沿厚度的分布 图 沿厚度的分布
Fig.5 Distribution of 6 Distribution of
along the thickness along the thickness

沿厚度的分布 图沿厚度的分布
Fig.7 Distribution of 8 Distribution of
along the thickness along the thickness

由图4和图6可以看出,算例2进一步表明辛元(CSE20)的面外剪切应力的结果分别比20节点位移元(CDE20)的结果更接近精确解

同时,该例题表明,层合板面外法向应力面内法向应力和面内剪切应力沿厚度方向的非线性分布特性不显著20节点辛元的这些应力结果比20节点位移元的结果相对理想,参见图6~820节点辛元的面外应力与精确解的误差为1.27%,是面外应力中的最大误差

图9和10中的数据结果采用以下公式进行了无量纲化处理:

(12)

沿厚度的分布 图沿厚度的分布
Fig.9 Distribution of 10 Distribution of
along the thickness along the thickness

20节点辛元的位移与精确解的误差为0.82%,是位移的最大误差

3 结 论

1) 本文推导20节点六面体辛元时已说明面外应力变量和位移变量采用了相同的高次插值函数因此,没有必要再引入单元内部非协调位移正因如此,与文献[14-15]中建立8节点六面体非协调辛元的过程比较,本文的推导过程非常简便另一方面,与基于最小余能原理的经典应力元比较,本文应力变量的插值多项式不需要满足平衡方程,形式简单

2) 20节点辛元避免了传统混合法中系数矩阵主对角线上含有零元素的弊端,因而保证了其有限元线性系统数值结果的精度和稳定性

3) 由于3个坐标方向均被离散处理,因此本文的20节点辛元对分析对象的几何形状没有限制

4) 辛元的有限元线性系统包含了位移和应力两类变量,为面外应力边界条件的引入提供了方便数值实例表明,20节点辛元面外应力结果的精度和位移结果的精度比较接近

参考文献(References):

[1] REDDY J N.A simple higher-order theory for laminated composite plates[J].Journal of Applied Mechanics, 1984, 51(4): 745-746.

[2] NOOR A K, BURTON W S.Assessment of shear deformation theories for multilayered composite plates[J].Applied Mechanics Reviews, 1989,42(1): 1-13.

[3] PIAN T H H, WU C C.Hybrid and Incompatible Finite Element Methods[M].Chapman and Hall CRC, 2005.

[4] JING H S, LIAO M L.Partial hybrid stress element for the analysis of thick laminated composite plates[J].Computers & Structures, 1990, 36(1): 57-64.

[5] 田宗漱, 卞学鐄.多变量变分原理与多变量有限元方法[M].北京: 科学出版社, 2011.(TIAN Zhongshu, PIAN T H H.Variational Principle With Multi-Variables and Finite Element With Multi-Variables[M].Beijing: Science Press, 2011.(in Chinese))

[6] REDDY J N, ROBBINS D H.Theories and computational models for composite laminates[J].Applied Mechanics Reviews, 1994, 47(6): 147-169.

[7] 钟万勰.应用力学的辛数学方法[M].北京: 高等教育出版社, 2006.(ZHONG Wanxie.Symplectic Solution Methodology in Applied Mechanics[M].Beijing: Higher Education Press, 2006.(in Chinese))

[8] 唐立民, 褚致中, 邹贵平, 等.混合状态Hamilton元的半解析解和叠层板的计算[J].计算力学学报, 1992, 9(4): 347-360.(TANG Limin, CHU Zhizhong, ZOU Guiping, et al.The semi-analytical solution of mixed state Hamilton element and the computation of laminated plates[J].Chinese Journal of Computational Mechanics, 1992, 9(4): 347-360.(in Chinese))

[9] ZOU G P, TANG L M.A semi-analytical solution for laminated composite plates in Hamilton system[J].Computer Methods in Applied Mechanics and Engineering, 1995, 128(3/4): 395-404.

[10] QING G H, QIU J J, LIU Y H.Free vibration analysis of stiffened laminated plates[J].International Journal of Solids and Structures, 2006, 43(6): 1357-1371.

[11] ATLURI S N, GALLAGHER R H, ZIENKIEWICZ O C.Hybrid and Mixed Finite Element Methods[M].Chichester: Wiley, 1983.

[12] ARNOLD D N, FALK R S, WINTHER R.Mixed finite element methods for linear elasticity with weakly imposed symmetry[J].Mathematics of Computation, 2007, 76: 1699-1724.

[13] LIAO C L, TSAI J S.Partial mixed 3-D element for the analysis of thick laminated composite structures[J].International Journal For Numerical Methods in Engineering, 1992, 35(7): 1521-1539.

[14] QING G, TIAN J.Highly accurate symplectic element based on two variational principles[J].Acta Mechanica Sinica, 2018, 34(1): 151-161.

[15] 刘艳红, 李锐.含参数辛元与热弹性复合材料层合板分析[J].复合材料学报, 2019, 36(5): 1306-1312.(LIU Yanhong, LI Rui.Parametered symplectic element and analysis of thermoelastic composite laminates[J].Acta Materiae Compositae Sinica, 2019, 36(5): 1306-1312.(in Chinese))

[16] PAGANO N.Exact solutionsfor composite laminates in cylindrical bending[J].Journal of Composite Materials, 1969, 3(3): 398-411.

[17] 范家让.强厚叠层板壳的精确理论[M].北京: 科学出版社, 1996.(FAN Jiarang.Exact Theory of Laminated Thick Plates and Shells[M].Beijing: Science Press, 1996.(in Chinese))

20-Node Hexahedron Symplectic Elements for Stress Analysis of Composite Laminates

TONG Yao1, YAO Yuzhe2

(1.Shanghai Aircraft Design & Research InstituteShanghai 201206, P.R.China; 2.College of Aeronautical Engineering, Civil Aviation University of China, Tianjin 300300, P.R.China)

Abstract: Usually, the accuracy of the stresses obtained with the conventional displacement finite element method is one-order lower than that of the displacements, and the out-of-plane stresses can hardly meet the continuity requirements.Then, combined with the minimum potential energy principle and the H-R variational principle, a 20-node hexahedral symplectic element involving displacement and out-of-plane stress variables was established.Incompatible displacement terms are needless in the element since the 2 kinds of variables are approximated with higher-order interpolation functions.Hence, the derivation process of the theory is very simple.Unlike in the partially mixed Hamiltonian element, the variables involved in the symplectic element are discretized in 3 coordinate directions without restriction of the element thickness and the structure geometry.Numerical examples show that, the 20-node symplectic elements exhibit stable convergence.Under the coarse mesh, the out-of-plane stresses obtained with the proposed element are closer to the exact solution than those by the incompatible 8-node symplectic element.

Key words: composite material; out-of-plane stress; H-R variational principle; 20-node symplectic element; displacement method; mixed method

*收稿日期: 2019-09-20;

修订日期:2019-11-22

基金项目: 国家自然科学基金(11502286)

作者简介:

童瑶(1980—),女,高级工程师, 硕士(E-mail: tongyao@comac.cc);

姚玉喆(1993—), 男,硕士生(通讯作者.E-mail: 727394427@qq.com).

ⓒ 应用数学和力学编委会,ISSN 1000-0887

http://www.applmathmech.cn

引用格式: 童瑶, 姚玉喆.基于20节点辛元的复合材料层合板应力分析[J].应用数学和力学, 2020, 41(5): 509-516.

中图分类号: O342; O343

文献标志码:A

DOI: 10.21656/1000-0887.400283

Foundation item: The National Natural Science Foundation of China(11502286)