工学第八章 复合材料结构分析有限元法

第八章复合材料构造分析有限元法§8.1引言§8.2考虑剪切效应的层合梁单元§8.3具有离散克希霍夫假定的层合板单元〔DKT〕§8.4考虑剪切效应的分项插值单元§8.5以一阶剪切理论为根底的Panda单元§8.6层合板的三维单元§8.7分层分析理论为根底的有限层单元§8.1引言根据板的理论,复合材料层合板的有限元，可分为:〔1〕以Kirchhoff假定的一阶理论层合板单元〔2〕以Reissner-Mindlin假定的一阶剪切层合板单元〔FOSDLT〕〔3〕以高阶理论为根底的层合板单元〔HOSDLT〕〔4〕以分层转化理论为根底的层合板单元〔5〕以三维理论为根底的层合板单元假设以变分原理分类又可分为：〔1〕以势能原理为根底的位移元〔2〕以余能原理为根底的应力元〔3〕以广义变分原理为根底的混合元或杂交元二.考虑一阶剪切变形Mindlin理论的势能和余能表示式(1)势能表达式：（1）extensionalstiffnesscouplingstiffnessbendingstiffnessshear（outofplane）stiffness(2)余能表达式〔2〕extensionalcompliancecouplingcompliancebendingcomplianceshear（outofplane）compliance三.广义变分原理(1)HellingerReissner原理〔二变量〕(2)胡海昌-鹫津广义变分原理〔三变量〕符号规定采用右手法则----向量符号XYZO

考虑剪切效应的层合梁单元以直角梁断面为例一维问题，设〔1〕位移场表达式〔3〕〔2〕应变表达式〔a〕面内〔b〕出平面将〔3〕代入上式，得〔3〕单层板的本构关系对第K层板整体坐标应力应变关系，并注意〔a〕面内〔4〕经约化可得（5）其中即（6）〔b〕出平面上式可约化为〔7〕其中〔8〕〔4〕层合梁本构关系〔a〕面内（9）其中（b）出平面（10）〔5〕势能表达式〔11〕〔12〕〔13〕〔14〕（6）刚度阵（a）设位移插值函数为线性插值其中（b）为二次插值∵对应弯曲部分将（13）（14）代入（11）得132L〔15〕〔16〕〔c〕设W用Hermet插值应变—节点位移阵刚度阵节点位移向量其中（17）将（17）和（14）式代入（12）式得（18）（19）刚度阵应变—节点位移阵〔20〕〔d〕在单刚装配时必须先凝聚掉

那么最后节点位移向量为:节点位移向量123456§7.3具有离散克希霍夫假定的层合板单元（DKT元）DiscreteKirchhoffTriangularElement

该单元是法国贡比涅大学Batoz教授提出的，该单元特点仅在三角形单元边中点满足Kirchhoff假定，故放松了原位移的约束，使此单元不存在剪切闭锁（即ShearLocking），且具有收敛好的特点，该单元还能保证各单元相邻边位移的协调性。但是该单元原来适用于各向同性板。我们将该单元发展成为复合材料层合板单元。

已知三角形单元其节点为i、j、k，定义为1、2、3，边中点为4、5、6，位移和转角均定义为向量，按右手定律规定其正负。O、X、Y为整体坐标。（1）内力势能—应变能

∵不考虑剪切效应，故（21）Z设节点位移向量分成两局部假设面内应变—节点位移阵为弯曲应变—节点位移阵为那么〔22〕将〔22〕式代入〔21〕式〔23〕由最小势能原理可得其中（24）（25）问题如何确定阵（2）阵的确定采用常应变插值函数，则（26）为面积坐标L1=1L3=0L1=1L1=0L3=1L2=0L1=constL2=constL3=constA1A2A3o〔27〕将（26）代入上式得其中该单元能保证各边连续〔3〕阵确实定〔28〕弯曲能量表达式中包含的导数项，故其插值函数要求满足沿边界上曲率，位移及转角（0，1，2）连续，其次,必须避免ShearLocking现象.令0为形心，边中点分别为4，5，6。设仅在节点1，2，3，4，5，6上满足Kirchhoff假定，即:〔29〕〔b〕建立边中点的转角与两端点转角的关系假设任选一个边，其两端节点为i，j，边中点为m，引入局部坐标ξοη，令w在ij边按Hermite插值，即（a）对分项插值其中为辅助自由度，最后需消去。〔30〕其中L为ij边长，S=ξ/l，无量纲则（31）当S=0.5为m点（边中点）（32）ijmxy0假设设沿ij边线性插值〔33〕在S=0.5处为m点〔34〕由〔32〕与〔34〕可合并写为:〔35〕〔36〕将转换到整体坐标上去

设在端点与边中点满足克希霍夫假定,即:〔37〕同理可得〔38〕将〔38〕代入〔29〕式〔39〕〔40〕将（39）代入（28），得将（27）（40）代入（25）可得单元刚度阵〔4〕几何刚度阵〔41〕代入〔41〕可得〔42〕则几何刚度阵（5）数例受均压四边简支边界为三角形板设层数精确解DKT误差456.0957.101.8%662.4663.060.9%§7.4考虑剪切效应的分项插值单元（1）内力势能→应变能→刚度阵（43）若板是薄板或忽略剪切效应，则阵必须是奇异阵，即剪切能量为零。而面内，弯矩耦合刚度阵为非奇异阵，则总单元刚度阵也必须是非奇异阵。设为广义位移,在单元内位移采用分项插值，且等参插值。（44）其中为i节点位移向量，即设应变位移阵为〔45〕其中将〔45〕代入〔43〕，并对节点位移向量变分得(46)其中（2）刚度阵奇异性的判别a、秩的判别规则在刚度阵求解时必须采用数值积分，一般采用高斯积分方法。然而采用不同的高斯积分方式所得到的刚度阵，是具有不同性质的特征。考虑剪切效应的板单元，要求剪切刚度阵为奇异，而其它刚度阵必须非奇异，总单元刚度阵必须非奇异。判别刚度阵奇异与非奇异，首先可判断它是缺秩或满秩。若缺秩为奇异。令K为n阶方阵，其满秩，则其秩＞n。①矩阵相乘的秩

若

则秩②矩阵相加的秩

若

则秩b、单刚的秩其中为权系数,D为弹性阵，其d×d，秩为d。2Dd=33Dd=6轴对称d=5B应变—节点位移向量，其尺寸,为节点未知数〔自由度〕故秩为d为高斯积分阶数。故,秩为假设构造由M个单元组成那么总刚度阵的秩为:〔47〕为总刚度非奇异性判别准那么。〔M为单元总数〕〔3〕数例以四与八节点平面等参元为例d=3每节点两个自由度(黑字为非奇异,红字为奇异)线性单元（四节点）二次单元（八节点）nMnMa4×2-3=5113奇异8×2-3=131212奇异b6×2-3=9216奇异13×2-3=232224非奇异c25×2-18=3216148非奇异[(5×9)+(4×5)]×2-34=96162192非奇异

四节点(hg=1)八节点(hg=2)

(a)(b)(c)奇异奇异奇异非奇异非奇异非奇异〔4〕简易刚度阵秩的判别准那么构造考虑剪切效应板单元，假设仅判别单刚奇异性，可采用如下简易判别条件，即：奇异〔48〕J

为增加一个单元形成的自由度的增加d为弹性阵的秩为高斯积分阶数其假定原构造单元总刚是满秩，在此根底上再增加一单元，去判别此单元增加后引起秩的变化。该方法计算秩可能偏大，因为在边界上增加一单元，其节点自由度的增加。将小于j〔5〕例构造考虑剪切效应单元单元形式j总单刚LR311×1×3=311×1×2=2522×2×3=1211×1×2=214QS922×2×3=1222×2×2=820QLR1222×2×3=1222×2×2=820CSR1533×3×3=2733×3×2=1845CLR2733×3×3=2733×3×2=1845黑字为非奇异,红字为奇异LRQSRQLRCSRCLR结论〔1〕三种Lagrange〔LR，QLR，CLR〕减缩积分比Serendipity单元〔QSR，CSR〕有更好的性能。而QSR剪切刚度阵为奇异。总刚度阵为非奇异，而CSR也经常剪切刚度阵为非奇异，故经常出现剪切Locking现象。〔2〕对厚板这两种积分都将得到较好的结果。§7.5以一阶剪切理论为根底的Panda单元本单元是在采用等参分项插值根底上，再引入厚度概念〔ThicknessConcept〕是一种退化单元。设一八节点壳单元，每节点仍为五个自由度。〔1〕由等参定义设（49）7O234568中面为自然坐标下ξ，η中i点的形函数。T为板总厚度。

为i节点处厚度，-1<ζ<1，为第j层厚度则位移场插值:（50）由应变位移关系:〔51〕〔2〕由三维弹性力学本构关系〔令，对本构方程约化〕，即〔51〕〔3〕势能原理与单元刚度阵〔52〕Jocobi阵

由于每层材料层内相同，而沿厚度方向各层不同，故引入ζ（局部局部坐标（每一层一个）

，ζ不能是连续函数，第K层为，

其原点在每层的中面。

由几何关系可知〔53〕那么〔54〕将〔53〕〔54〕代入〔52〕得〔55〕n为层板总层数。〔4〕例0/90/90/0，等厚简支方板，外载:力学性能如下:结果见图。102030402468a/tCLT51015201.52.02.53.0a/tCLTPANDA元解§7.6层合板的三维单元〔20节点块体元〕（1）位移与几何插值

等参插值（56）0其中在边中点：（2）本构关系各层采用横观各向同性假定（材料主轴）（57）在统一坐标下，（58）注意：z轴不转动为铺设角。（3）单刚（59）∵每层材料性质不同，与panda元相同，引入厚度概念，即令（60）t为单元总厚度，为第k层厚度（60）式代入（59）得，（61）§7.7分层理论为基础的有限层元法

若边界形状为规则层板，且边界条件比较简单（例如四边简支，固定等），则我们可采用分层理论为基础的有限层法。此方法特征，沿X—Y平面以符合边界条件位移函数展开，而在厚度方向采用有限元法，形成有限层法。（1）设位移函数，第

i层（62）其中为每层厚度，而下标t、c、b表示第i层上表面、中面与下表面。分别为满足边界条件的位移函数。（2）本构关系（第i层）（63）（3）刚度阵其中dξdη部分可积分积出，而沿厚度z