清华大学计算固体力学第九次课件梁和壳.ppt

非线性有限元第9章 梁和壳,计算固体力学,第9章 梁和壳,引言梁理论基于连续体的梁CB梁CB梁的分析基于连续体的壳CB壳CB壳理论剪切和膜自锁假设应变单元一点积分单元,1 引言,第8章介绍了平面单元(二维)和实体单元(三维),在二维问题中，最经常应用的低阶单元是3节点三角形和4节点四边形。在三维单元中，是4节点四面体和8节点六面体单元。,结构单元可以分类为：梁，运动由仅含一个独立变量的函数描述；壳，运动由包含两个独立变量的函数描述；板，即平面的壳，沿其表面法线方向加载；膜，面内刚度很大，面外刚度很小的薄壳。,1 引言,1 引言,在工程构件和结构的模拟中，梁和壳及其他结构单元是极为有用的。应用薄壳，如汽车中的金属薄板，飞机的机舱、机翼和风向舵；以及某些产品的外壳，如手机、洗衣机和计算机。用连续体单元模拟这些构件需要大量的单元，如采用六面体单元模拟一根梁沿厚度方向至少需要5个单元，而既便采用低阶的壳单元也能够代替5个或者更多个连续体单元，极大地改善了运算效率。应用连续体单元模拟薄壁结构常常导致较高的宽厚比，从而降低了方程的适应条件和解答的精度。在显式方法中，根据稳定性的要求，采用连续体单元的薄壁结构被限制在非常小的时间步。,1 引言,通过两种途径建立壳体有限元：1 应用经典壳方程的动量平衡(或平衡)的弱形式；2 结构的假设直接由连续体单元建立基于连续体(CB)方法。,第一种途径是困难的，尤其是对于非线性壳，因为对于非线性壳的控制方程是非常复杂的，处理起来相当不方便；它们的公式通常由张量的曲线分量来表示，并且其特征，诸如厚度、连接件和加强件的变量一般也是难以组合。而且对于什么是最佳的非线性壳方程的观点也不一致。第二种CB方法(基于连续体)是直观的，得到非常好的解答，它适用于任意的大变形问题并被广泛地应用于商业软件和研究中。因此，我们将关注CB方法。这种方法也称为退化的连续体方法。,1 引言,在大多数板壳理论中，通过强制引入运动假设建立平衡或者动量方程，然后应用虚功原理推导偏微分方程。在CB方法中，在连续体弱形式的变分和试函数中强制引入运动假设。因此，对于获得壳和其它结构的离散方程，CB壳方法更加直观。在关于壳的CB方法中，由两种途径强化运动假设：1）在连续体运动的弱形式中，或者2）在连续体的离散方程。,由二维梁描述CB方法编程特点，应用第一种途径的理论，检验CB梁单元。建立CB壳单元，编程，发展CB壳理论，结合由于大变形在厚度上变化的处理方法，给出在三维问题中描述大转动的方法。CB壳单元的两点不足：剪切和膜自锁。将描述假设应变场的方法防止发生自锁，给出了缓和剪切和膜自锁的单元例子。描述应用在显式程序中的4节点四边形壳单元一点积分单元。这些单元是快速和强健的，并且适用于大规模问题的计算。,当结构一个方向的尺度(长度)明显大于其它两个方向的尺度，并且沿长度方向的应力最重要时，可以用梁单元模拟。梁理论的基本假设是：由一组变量可以完全确定结构的变形，而这组变量只是沿着结构长度方向位置的函数。应用梁理论获得可接受的结果，横截面尺度必须小于结构典型轴向尺度的1/10。典型的轴向尺度为：支承点之间的距离；横截面发生显著变化部分之间的距离；所关注的最高阶振型的波长。梁单元假设在变形中垂直于梁轴线的横截面保持平面。不要误解横截面的尺度必须小于典型单元长度1/10的提法。高度精细的网格中可能包含长度小于其横截面尺寸的梁单元(尽管一般不建议这样做)，在这种情况下实体单元可能更适合。,2 梁理论,2 梁理论,梁理论的假设,运动学假设关注梁的中线(也称为参考线)的运动。由垂直于中线定义的平面称之为法平面。,梁横截面几何形状,广泛应用的梁理论有两种：其运动学假设是：Euler-Bernoulli梁：假设中线的法平面保持平面和法向；称为工程梁理论，而相应的壳理论称为Kirchhoff-Love壳理论。Timoshenko梁：假设中线的法平面保持平面，但不一定是法向；称为剪切梁理论，相应的壳理论称为Mindlin-Reissner壳理论。,2 梁理论,梁理论的假设,考虑一点P的运动，它在中线上的正交投影为点C。如果法平面转动视为一个刚体，则P点的速度相对于C点的速度给出为,2 梁理论,Timoshenko梁理论,在二维问题中，角速度的非零分量是z 分量，所以,法线的角速率,相对速度为,中线上任何一点的速度是 x 和时间 t 的函数，因此有,即梁上任何一点的速度是相对速度和中线速度之和,2 梁理论,Timoshenko梁理论,应用变形率的定义,变形率的非零分量只有轴向分量和剪切分量，后者为横行剪切。,由于梁内的变形率是有限的，非独立变量 和 只要求C0 连续，位移（挠度）和截面转动各自独立，使截面发生剪切变形后保持平面。,Euler-Bernoulli理论,运动学假设是法平面保持平面和法向，因此，法线的角速度是由中线的斜率的变化率给出,上式等价于要求剪切分量为零，表示在法线和中线之间的夹角没有变化，即法线保持法向。轴向速度则给出为,变形率给出为,2 梁理论,注意在上式中的两个特征：1）横向剪切为零；2）在变形率的表达式中出现了速度的二阶导数，梁内的变形率是有限的，即非独立变量的速度场必须为C1连续。,Euler-Bernoulli理论,2 梁理论,E-B梁理论常称为C1 理论，因为它要求C1 近似。转角由位移对坐标的导数给出(区别于Tim梁位移与转角相对独立)。梁单元常常是基于E-B理论，在一维情况下，C1 插值是很容易构造的。,Timoshenko梁有两个非独立变量(未知)，在E-B梁中只有一个非独立变量。类似的简化发生在相应的壳理论中：在Kirchhoff-Love壳理论中只有3个非独立变量，而在Mindlin-Reissner壳理论中有5个非独立变量（经常应用6个）。,E-B梁理论要求C1 近似是E-B和Kirchhoff-Love理论的最大缺陷，在多维空间中C1 近似是很难构造的。由于这个原因，在软件中除了针对梁之外很少应用C1 构造理论。,2 梁理论,横向剪切在厚梁中是明显的，在Timoshenko梁和Mindlin壳中常常应用。当梁趋于薄梁时，Timoshenko梁中的横向剪切在理想性能单元情况将趋于零。因此，在数值结果中也观察到了垂直假设，它隐含着对于薄梁横向剪切为零。这些假设主要是以实验为依据的：这一理论预测与实验测量相吻合。对于弹性材料，梁的闭合形式解析解也支持这一理论。它带来的好处是在有限元程序中，用中厚壳代替薄壳，用铁摩钦柯梁代替伯努利梁。,3 基于连续体的梁CB梁,为什么要建立CB梁和CB壳：1 梁与板壳组合的偏置(offset)2 接触问题的处理3 边界条件的处理,通过指定一个偏置量，可以引入偏置。偏置量定义为从壳的中面到壳的参考表面的距离与壳体厚度的比值。,梁作为壳单元的加强部件：（a）梁截面无偏置（b）梁截面有偏置,3 基于连续体的梁CB梁,建立CB二维梁的公式，结构的控制方程与连续体的控制方程是一致的：质量守恒线动量和角动量守恒能量守恒本构方程应变-位移方程,右为CB梁单元，左为母单元。连续体单元的节点仅在顶部和底部，,在 方向的运动一定是线性的。这些节点称为从属节点,3 基于连续体的梁CB梁,为常数的线称为纤维，沿着纤维的单位矢量称为方向矢量,为常数的线称为迭层,主控节点,3 基于连续体的梁CB梁,在纤维将从属节点与参考线连接的内部截面上，引入主控节点，其自由度描述了梁的运动。以主控节点的广义力和速度建立运动方程。在一条纤维上，每一主控节点联系一对从属节点，三点共线。,3 基于连续体的梁CB梁,假设：1 纤维保持直线；2 横向正应力忽略不计，即平面应力条件；3 纤维不伸缩。,第一个假设与经典的Mindlin-Reissner假设中要求法线保持直线是不同的，纤维可以不垂直于中线，称其为修正的M-R假设。,如果CB梁单元近似地为Timoshenko梁，其纤维方向尽可能地接近中线的法线方向是必要的，通过指定从属节点的初始位置可以实现这一点。否则，CB梁单元的行为将从根本上偏离Timoshenko梁，并且可能与所观察到的梁的行为不一致。,3 基于连续体的梁CB梁,注意到纤维的不可伸缩仅适用于运动学描述，不适用于动力学描述。不可伸缩性与平面应力的假设相矛盾：纤维通常接近于y方向，如果，则必须考虑速度应变。,通过不使用运动，而是由本构方程来计算，消除了这种矛盾。令，由 计算沿厚度方向的变化。这等价于由物质守恒获得厚度，因为平面应力的本构方程与物质守恒有关。然后修正节点内力以反映沿厚度方向的变化。这样，不可伸缩性的假设仅仅适用于运动。,假设：1 纤维保持直线；2 横向正应力忽略不计，即平面应力条件；3 纤维不伸缩。,3 基于连续体的梁CB梁,运动：通过主控节点的平移x(t),y(t)和节点方向矢量的旋转描述运动,从x 轴逆时针旋转的转角为正,通过对连续体单元的标准等参映射，由从属节点运动给出梁的运动,连续体的标准形函数（在节点指标中*代表上节点或下节点）,为了使上面的运动与修正的M-R假设相一致，基本连续体单元的形函数在 方向必须是线性的。因此，母单元在该方向只有两个节点，沿着纤维方向只能有两个节点。速度场为,3 基于连续体的梁CB梁,运动,在从属节点的运动中，现在强制引入不可伸缩条件和修正的M-R假设,pI为主控节点的方向矢量，h0 是伪厚度(初始厚度)，因为它是沿着纤维方向在单元的顶部与底部之间的距离。这是连续体单元向CB梁单元转化的关键一步。,当前节点的方向矢量给出为,总体基矢量,3 基于连续体的梁CB梁,从属节点的速度是坐标的材料时间导数，服从,由每个节点的三个自由度描述主控节点的运动,运动,写出矩阵的形式为,上标slave和mast强调连续体节点是从属节点，梁节点为主控节点。,3 基于连续体的梁CB梁,由于从属节点速度是与主控节点速度相关的，所以节点力的关系为,节点力与主控节点的速度是功率耦合的，在I 处,节点力,可以看出,为了将标准连续体单元转化为CB梁单元，必须强化平面应力假设。采用应力和速度应变的层间分量是方便的。构造每层的基矢量为,3 基于连续体的梁CB梁,与迭层正切,垂直于迭层,本构更新,3 基于连续体的梁CB梁,在迭层分量上加“帽子”，它们随着材料转动，因此考虑是共旋的。变形率的迭层分量给出：,在应力计算中，必须观察平面应力约束,如果本构方程是率的形式，则约束是,例如，对于各向同性次弹性材料，应力率分量给出为,求解得到,3 基于连续体的梁CB梁,求解得到,对于更为一般的材料(包括模型中缺少对称性的定律，诸如非关联塑性的材料)，本构率关系可以写成为,为切线模量，最后一个方程强调了平面应力条件，求解 Dyy,3 基于连续体的梁CB梁,节点内力,除了强制平面应力条件外，从属节点内力采用与连续体单元节点内力相同的计算。积分由数值积分求得。在CB梁中既不应用完全积分公式，也不应用选择减缩积分公式(4-5-33)。这两种方法都会导致剪切自锁(见第7节)。,在2节点单元中，在,处采用单一束积分点，可以避免剪切自锁。,这种积分方法也称为选择减缩积分。它能精确地积分求得轴向正应力，但是不能准确地积分求得横向切应力,3 基于连续体的梁CB梁,沿方向的积分点数目依赖于材料定律和对精度的要求。1 平滑的超弹性材料定律，3个积分点是足够的。2 弹-塑性材料，应力分布不是连续可导至少需要5个积分点。对于弹-塑性材料定律，沿方向的Gauss积分并不是最佳选择，因为这些积分方法是基于高阶多项式的插值，其默认假设数据是平滑的。所以，对于非光滑函数，常常采用梯形规则，其运算效率更高。,3 基于连续体的梁CB梁,为了说明在剪切自锁情况下，选择减缩积分的过程，考虑一个基于4节点四边形连续体单元的2节点梁单元。通过对在,处一串积分点的积分得到节点力：,质量矩阵,CB梁单元的质量矩阵可以由转换公式得到：,运动方程,对于对角化质量矩阵，在一个主控节点的运动方程为,4 CB梁的分析,例9.1 2节点梁单元,应用CB梁理论建立基于4节点四边形连续体的2节点CB梁单元，将参考线(中线)置于上下表面的中间位置，将主控节点放置在参考线与单元边界的交点处，从属节点是角点。,主控节点,从属节点,4节点连续体单元的运动,上述的运动得到,令,也给出4节点连续体单元的运动,例9.1 2节点梁单元,4 CB梁的分析,所有纤维的不可伸长性 尽管节点的纤维是不可伸长的，在一个单元中的其它纤维可能发生长度的变化。通过图示的指定条件，不用任何方程就可以看到：在中点处的纤维明显变短了。,节点力：主控节点力由从属节点力给出,4 CB梁的分析,计算上式得到,这个变换给出了平衡的结果：主控节点力是从属节点力的合力；主控节点力矩是从属节点力绕主控节点的力矩。,Green应变 以PK2应力和Green应变的形式应用于本构方程。,从属节点的位置,4 CB梁的分析,通过取节点坐标的差值得到从属节点的位移。通过,给出的连续体位移场，则可以得到任何点的位移。通过本构关系和PK2应力则可以计算Green应变。,初始和当前构形的方向矢量给出为,矩形单元的速度应变 当基本连续体单元为矩形，并且梁的中线沿着x轴时，由于方向矢量是沿着y方向,4 CB梁的分析,用一维形式写出上式的分量，线性形状函数给出,其中,速度应变分量则给出为：,由平面应力条件计算分量,剪切自锁,4 CB梁的分析,为了检验剪切自锁的原因，考虑例9.1的2节点梁单元。令单元位于x 轴方向，线性响应，因此在运动学中，用线性应变,代替,用位移代替速度。由公式的对应部分给出横向剪切应变：,考虑在纯弯状态下的单元有,对于这些节点位移，给出,由平衡方程，当力矩为常数时，剪力为零。,但是在大多数单元中，横向剪切应变和剪切应力不为零,事实上，除了在,外它们处处不为零。出现附加剪切(见上图)。,剪切自锁,4 CB梁的分析,这附加的横向剪切对于单元的性能具有很大的影响。为了解释这种影响的严重性，对于一个单位宽度矩形横截面的线弹性梁，检验与弯曲和剪切应变有关的能量。上面节点位移的弯曲能量给出为,关于梁的剪切能量给出为,这两种能量的比值为,当,剪切能量是显著地大于弯曲能量。由于在纯弯曲中剪切能量应该为零，这附加剪切能量吸收了大部分的能量。其结果是明显地低估了总体位移。采用单元细划使得剪切自锁的单元可以收敛于精确解，只是非常慢。而体积自锁根本得不到收敛结果。,剪切自锁,4 CB梁的分析,由方程立即提出问题，在这些单元中为什么不采用不完全积分消除剪切自锁：注意到在 处横向剪力为零，这对应于在一点积分中的积分点。因此，通过对剪切相关项的不完全积分消除了伪横向剪切。,相对于2节点梁，在采用二次插值的3节点梁中，剪切自锁是很不明显的。考虑一个长为 l 的3节点梁单元，采用母坐标,在该单元中的剪切应变给出为,剪切自锁,4 CB梁的分析,考虑纯弯曲的状态，,将这些节点位移代入公式，证明在单元中横向剪切为零。基于这个结果，这里没有理由出现自锁。但是，考虑另外一种变形,由于法线保持法向，剪力应该为零。但是，公式却给出了横向剪切，相应于这种变形的节点位移是,所以，除了，有限元数值近似处处给出了非零剪力。因此，对于自由剪切(纯弯)模式，横向剪切将发生在这种单元中，而在模拟薄梁时，它将是无效的。,当结构一个方向的尺度(厚度)远小于其它方向的尺度，并忽略沿厚度方向的应力时，可以用壳单元模拟。例如，压力容器结构的壁厚小于典型整体结构尺寸的1/10，一般用壳单元进行模拟。以下尺寸可以作为典型整体结构的尺寸：支撑点之间的距离；加强件之间的距离或截面厚度有很大变化部分之间的距离；曲率半径；所关注的最高阶振动模态的波长。壳单元假设垂直于壳面的横截面保持为平面。请不要误解为在壳单元中也要求厚度必须小于单元尺寸的1/10，高度精细的网格可能包含厚度尺寸大于平面内尺寸的壳单元（尽管一般不推荐这样做），实体单元可能更适合这种情况。,5 基于连续体的壳CB壳,5 基于连续体的壳CB壳,两种壳单元：1 常规壳单元2 基于连续体的壳单元,通过定义单元的平面尺寸、表面法向和初始曲率，常规的壳单元对参考面进行离散。但是，常规壳单元的节点不能定义壳的厚度；还需要通过截面性质定义壳的厚度。基于连续体的壳单元类似于三维实体单元，它们对整个三维物体进行离散和建立数学描述，其动力学和本构行为类似于常规壳单元。对于模拟接触问题，基于连续体的壳单元与常规的壳单元相比更加精确，因为它可以在双面接触中考虑厚度的变化。然而，对于薄壳问题，常规的壳单元提供更优良的性能。,壳体公式-厚壳或薄壳,壳体问题一般归结为两类之一：薄壳和厚壳。厚壳假设横向剪切变形对计算结果有重要的影响；薄壳假设横向剪切变形小到足以忽略。图(a)描述了薄壳的横向剪切行为：初始垂直于壳面的材料线在整个变形过程中保持直线和垂直。因此，横向剪切应变假设为零。图(b)描述了厚壳的横向剪切行为：初始垂直于壳面的材料线在整个变形过程中并不要求保持垂直于壳面，因此，发生了横向剪切变形。,5 基于连续体的壳CB壳,5 基于连续体的壳CB壳,经典壳理论中的假设,两种壳理论中，运动学假设为：Kirchhoff-Love理论：中面的法线保持直线和法向。Mindlin-Reissner理论：中面的法线保持直线，但不是法向。实验结果表明，薄壳满足Kirchhoff-Love假设。较厚的壳或组合壳体，Mindlin-Reissner假设是更为合适的，横向剪切的效果特别重要。Mindlin-Reissner理论也可以应用于薄壳中：在这种情况下，法线将近似地保持法向，并且横向剪切将几乎为零。,厚度与曲率半径的比值的条件是对于壳理论适用性的重要要求。当它不满足时，壳理论将是不适用的。,5 基于连续体的壳CB壳,1 纤维保持直线（修正的M-R假设）；2 垂直于中面的应力为零（也称为平面应力条件）；3 动量源于纤维的伸长，和沿纤维方向忽略动量平衡。假设1与经典的Mindlin理论不同之处在于约束纤维保持直线，而不是法向。必须布置节点，使纤维方向尽可能地接近于法线。,CB壳理论中的假设,在CB壳理论中常常认为纤维是不可伸长的，这是与平面应力条件矛盾的，当应用不可伸长的条件时，是为了忽略在p方向的相关运动的动量平衡，在计算节点内力时，厚度的改变是不能忽略的。,定义三种坐标系统：1.总体Cartesian坐标系统（x,y,z），应用基矢量,坐标系统,5 基于连续体的壳CB壳,2.旋转的层坐标系统,应用基矢量,3.与主控节点相关的节点坐标系统,表示在下表面和参考面之间沿着纤维方向的距离,表示在上表面和参考面之间沿着纤维方向的距离,表示厚度的改变，壳的厚度一般定义为在上下两个表面之间沿着法线的距离。,壳的材料方向,应用局部的直角、圆柱或者球坐标系，可以代替整体的笛卡尔坐标系,例如，如果在图中的圆柱中心线与整体坐标3轴一致，局部材料方向可以这样定义，使局部材料1方向总是沿着圆环方向，并使相应的局部材料2方向总是沿着轴方向。,5 基于连续体的壳CB壳,5 基于连续体的壳CB壳,在主控节点上的节点速度和力为,运动的有限元近似,以从属节点的形式，描述运动的有限元近似为,主控节点转换速度和方向矢量角速度,表示从属节点的速度,厚度的变化率,当前节点的方向矢量给出为,5 基于连续体的壳CB壳,在主控节点上的节点速度和力为,运动的有限元近似,以从属节点的形式，描述运动的有限元近似为,主控节点转换速度和方向矢量角速度,表示从属节点的速度,厚度的变化率,当前节点的方向矢量给出为,5 基于连续体的壳CB壳,运动的有限元近似,将从属节点速度联系到主控节点速度为,在主控节点力的计算中采用了当前厚度，因此考虑了纤维的伸长,5 基于连续体的壳CB壳,本构方程,连续介质材料的所有本构都可以应用于CB壳。但是，必须引入平面应力条件。可以应用强制引入约束的方法，如Lagrange乘子法和罚方法。以Voigt形式写出的率更新方程为：,切线模量矩阵是55和51子矩阵,通过消去第6个方程，可以获得与非零应力增量相关的修正矩阵,从第6个方程得到变形率分量，计算厚度的变化。,5 基于连续体的壳CB壳,厚度,可以直接或者由率形式得到厚度。在任意时刻的厚度给出为,这里给出的更新厚度提供了关于厚度的双参数近似。在等参CB单元中，由于变形梯度在厚度方向近似为线性，这通常是足够的。单参数形式经常仅用于说明厚度的平均变化。双参数形式是更精确的，因为当伸长时叠加弯曲，在压缩边和拉伸边的厚度改变是不同的。另一种更加精确的方法是计算所有积分点的新的位置，但通常是不必要的。,壳体厚度和截面点（section points）描述壳体的横截面必须定义壳体的厚度。此外，还要选择是在分析过程中还是在分析开始时计算横截面的刚度。如果选择在分析过程中计算刚度，采用数值积分法，在沿厚度方向的每一个截面点（section points）（积分点）上独立地计算应力和应变值，这样就允许了非线性的材料行为。例如，弹塑性材料的壳在内部截面点还保持弹性时，其外部截面点可能已经达到了屈服，4节点减缩积分壳单元中唯一的积分点位置和沿壳厚度方向上截面点的分布如图示。,5 基于连续体的壳CB壳,壳体厚度和截面点 当在分析过程中积分单元特性时，可指定壳厚度方向的截面点数目为任意奇数。对性质均匀的壳单元，一般在厚度方向上取5个截面点，对于大多数非线性设计问题是足够了。但是，对于一些复杂的模拟必须采用更多的截面点，尤其是当预测会出现反向的塑性弯曲时（在这种情况下一般采用9个截面点）。对于线性问题，3个截面点已经提供了沿厚度方向的精确积分。当然，对于线弹性材料壳，选择在分析开始时计算材料刚度更为有效。如果选择仅在分析开始时计算横截面刚度，材料行为必须是线弹性的。在这种情况下，所有的计算都是以整个横截面上的合力和合力矩的形式进行。如果要求输出应力或应变，输出在壳底面、中面和顶面的值。,5 基于连续体的壳CB壳,5 基于连续体的壳CB壳,主控节点力,在主控节点的内力和外力可以由从属节点力得到，即,由连续体单元的程序计算从属节点力,质量矩阵,利用基本连续体单元的质量矩阵,通过转换公式获得CB壳单元的质量矩阵。则66子矩阵给出为,转动惯量,平移质量,5 基于连续体的壳CB壳,离散动量方程,对于上面给出的对角化质量矩阵，在节点处的3个平动方程为,节点力和节点速度为,以节点坐标系统表示的3个转动方程为,上式就是著名的Euler运动方程。它们对于角速度是非线性的，但是对于一个各向同性转动质量矩阵，二次项将消失。,No sum on I,5 基于连续体的壳CB壳,切线刚度,由基本连续体单元刚度矩阵的变换，得到切线刚度和荷载刚度矩阵：,连续体单元的切线刚度矩阵,5个自由度公式,如果没有扭转，不考虑，或刚体转动对变形没有影响，每个节点处的壳的运动用5个自由度描述。主控节点的节点速度为,对于CB壳理论，5个比6个自由度的描述更加合适。当壳为平坦时，对于6个自由度的描述其刚度为奇异的。另一方面，5个自由度的描述必须在角点处进行修正，使其符合结构特点。对于在节点处采用可变化自由度数目的软件，仅在需要增加附加自由度的那些节点处应用6个自由度可能是最合适的，如连接处。,在Mindlin-Reisssner理论中，横向剪应力,沿着壳的厚度方向为常数。由于应力张量的对称性，在这些表面的横向剪力必须为零，除非一个剪切面力施加在上表面或者下表面。对于平衡状态下弹性梁的分析表明，沿梁的厚度方向，横向剪应力应该为二次的，在上下表面处为零。因此，常值剪切应力分布高估了剪切能量。通常采用一个修正因数，如矩形截面为5/6，已知为剪切修正，以减少与横向剪切相关的能量，并且对于弹性梁和壳可以做出关于这个因数的精确估计。然而，对于非线性材料，估计一个剪切修正因数是非常困难的。,结构理论的非协调性和特殊性,6 CB壳理论,结构理论的非协调性和特殊性,6 CB壳理论,在Kirchhoff-Love理论中的非协调性甚至更加严重，由于运动学的假设(平截面假设并垂直中面)，导致了横向剪力为零。在结构理论中，如果力矩不是常数，在梁中的剪力必须非零。因此，Kirchhoff-Love的运动学假设与平衡方程是矛盾的。但是，与实验结果比较证明它是相当精确的，并且对于薄的均匀壳，它恰与Mindlin-Reissner理论同样精确。在薄壁结构的变形中，横向剪力并没有起到重要的作用，是否考虑它的作用几乎没有影响。M-R单元简单，甚至当横向剪力的影响可以忽略时，应用M-R单元也满足精度。,修正的Mindlin-Reissner CB模型提供了产生误差的附加可能性。如果方向矢量不垂直于中面，则运动与实验观察到的运动将会有明显的偏差。当一个法向面力施加在壳的任何面上时，零法向应力的假设是矛盾的，这里是(这是学生经常问到的问题)。为了平衡，法向应力必须等于所施加的法向面力。然而，在结构理论中它们被忽略了，因为与轴向应力相比它们是非常小的；法向应力仅仅吸收了很小部分能量，对变形几乎没有影响。,结构理论的非协调性和特殊性,6 CB壳理论,sm=?,st=?,结构理论的非协调性和特殊性,6 CB壳理论,在壳体的分析中，要注意边界效应。某些边界条件导致了边界效应，在较窄的边界层处性能发生了剧烈的变化。对于某些边界条件，在边界的角点处可能发生奇异性（管道交接处）。应用结构运动学假设的一个原因是它们改善了离散方程的适应性。如果一个壳体由三维的连续体单元模拟，自由度是所有节点的平动，与厚度方向应变相关的自然模态具有非常大的特征值。其结果，对于一个隐式更新算法，线性化平衡方程或者线性化方程的适应性可能是非常差的。壳方程的适应性也不如标准连续体模型那样好，但是比薄壳的连续体模型的适应性要明显强一些。在显式方法中，由于沿厚度方向模态的较大特征值，导致薄壁结构的连续体模型需要非常小的临界时间步长。CB壳模型可以提供更大的临界时间步长。,7 剪切和膜自锁,壳单元的最大困难特性是剪切和薄膜自锁。剪切自锁源于出现了伪横向剪切。更确切地说，它源于许多单元没有能力表现弯曲变形，剪切刚度通常远远大于弯曲刚度，伪剪切吸收了大部分由外力产生的能量，而预计的挠度和应变成为非常小的量值，这就是所谓的剪切自锁。,薄膜自锁的出现是源于在壳单元中没有能力表现变形的不可伸长模式。壳弯曲而没有伸长：一张纸，能够很容易地将其弯曲，称为不可伸长弯曲。然而，用手拉伸一张纸几乎是不可能的。壳的行为类似：弯曲刚度很小，而薄膜刚度很大。当有限元没有伸长又不能弯曲时，能量是不准确地转换成为薄膜能，于是导致低估了位移和应变。在屈曲模拟中，薄膜自锁是尤为重要，因为许多屈曲模态是完全的或者接近于不可伸缩的。,7 剪切和膜自锁,剪切和薄膜自锁与体积自锁是相似的：当有限元近似的运动不能满足约束时，约束模式比正确运动的刚度表现的更为刚硬。在体积自锁的情况中，约束是不可压缩，体积刚度过大，而对于剪切和薄膜自锁，在弯曲中的约束为Kirchhoff-Love正常状态约束和不可伸长约束(与纤维的不可伸长没有任何关系)。,应该注意的是薄壳的自由剪切行为不是一个精确的约束。对于较厚的壳和梁，希望有某些横向剪切，但是，就像对于几乎不可压缩材料，体积自锁的单元表现很差一样，对于厚度适中的壳，即使当横向剪切出现时，壳单元的剪切表现是很差的。,7 剪切和膜自锁,自锁现象比较,不可压缩，等体积运动，J=常数,自锁：单元没有能力锁住不该有的变形 或单元没有能力表现该有的变形,Kirchhoff-Love约束,7 剪切和膜自锁,薄膜自锁,为了说明薄膜自锁，利用Maguerre浅梁方程,考虑一个长为 l 的3节点梁单元，采用母单元坐标为,在一个不可伸缩的模式中，薄膜应变,必须为零。,对公式中的表达式,进行积分，对于y0 令,考虑一个梁的纯弯模式，有,在没有横向剪切时，由公式可以得到,7 剪切和膜自锁,薄膜自锁,在没有横向剪切时，证明：,证明：,积分：,由边界条件：,得到：,7 剪切和膜自锁,薄膜自锁,设在一个初始对称构形中，,如果,要求满足公式。由公式计算薄膜应变给出为,则在这种变形的不可伸缩模式下，除了在,外,伸缩应变处处不为零。如果单元包括伸缩应变不为零的积分点，单元将展示薄膜自锁。,对于3节点的CB梁，在 中给出的剪切和在上面给出的薄膜应变，它们都在点 处为零，,即两点积分的Gauss积分点。这些点常称为Barlow点。,7 剪切和膜自锁,消除自锁,通过在积分点 处对剪切能进行不完全积分，限制对其他点的剪切能取值可以避免附加剪切，从而使单元不会自锁。多场方法，通过设计合适的应变场，也可以回避自锁。例如，如果应用Hu-Washizu弱形式，通过使得横向剪切为常数可以回避剪切自锁。横向剪切速度和剪切应力场为（Hybrid法）：,系数由离散协调方程和本构方程确定。,应用假设应变方法也可以避免自锁，实质是设计横向剪切场和薄膜应变场，从而使得附加剪切和薄膜自锁为最小。必须设计假设应变场，以保证刚度矩阵是正确的秩。对于2节点的梁，假设的应变场必须是常数，并且在纯弯时必须为零。实现这些目的，如果,7 剪切和膜自锁,消除自锁,是在中点处等于 的常数场，对于这个场，在纯弯时整个单元中的假设剪切应变率将为零。,如果应用单元表示一般的三次位移场，在22 Gauss点处的应力将与节点位移具有相同的精度。在设计壳单元中，这一发现证明是非常有用的。例如，由Zienkiewicz、Taylor等提出的减缩积分的成果。因此，关于应用两点Gauss积分（Barlow点）的剪切和薄膜功率，减缩积分将消除剪切和薄膜自锁。,8 假设应变单元,假设应变4节点四边形,例如横向剪切场的构造，,在壳单元中，通过假设应变方法和选择减缩积分也可以避免剪切和薄膜自锁。然而，这些方法的设计对于壳比对于梁或者连续体更加困难。例如，Hughes等(1978)应用选择减缩积分处理4节点四边形板单元，单元始终存在一个伪奇异模式，即w 沙漏模式。因此，选择减缩积分为连续体提供了强健的单元，而对于壳体是不成功的。,可以推论，对于一个纯弯曲的梁，如果横向剪切分布是线性的，并且在中间点为零，则它在常数场中的映射也为零。,除了 的点，它们处处不为零。在纯弯状态时出现的横向剪切常常称为附加剪切。,8 假设应变单元,假设应变4节点四边形,为任意参数。,考虑平面矩形壳单元，类似于梁，当弯矩施加到两端时，横向剪力为零（纯弯曲）；当材料为各向同性时，横向剪切应变为零。通过使剪切应变为常数可以满足这些条件，即令，。但是，一个常数的横向剪切将会导致秩缺乏，从而出现不稳定单元，为了保存稳定性，增加一个取决于y的项，给出横向剪切为，线性项对于弯曲力学性能没有影响，不会影响不自锁行为。,对于，也有类似的讨论。,扩展到四边形的母单元上，有,8 假设应变单元,单元的秩,上面的单元是否满秩？利用x-y面的平面壳元(板单元)来说明。我们仅考虑弯曲性能，每个节点则有3个相关的自由度：1个位移，2个转角()。单元有4个节点共12个自由度，考虑3个刚体运动，单元适当的秩是9。由第8章平面Q4单元的完全积分，8个自由度3个刚体位移5个线性独立场，另外，这里2个横向剪切有4个线性独立场，因此，总的线性独立场有5+4=9个，提供了适当的秩。,4节点四边形壳元关于剪切的插值点,8 假设应变单元,9节点四边形,9节点壳避免薄膜和剪切自锁的假设应变场，由内部点插值的假设速度应变，如图a,为a次一维Lagrange插值,8 假设应变单元,9节点四边形,在Gauss点，弯曲时横向剪切为零，而在不可伸缩的弯曲时薄膜应变为零。因此，单元不会展示附加的横向剪切或者薄膜应变，即不会产生自锁行为。,在假设速度应变场中的高次项,提供了稳定性。,由点插值的速度应变如图c 所示。,由点插值的剪切分量如图b 所示。,9 一点积分单元,在显式软件中，最常用的壳单元是一点积分的4节点四边形。一点积分是指在参考面上积分点的数目(高斯点)，在厚度方向的任何位置可以采用3到30或者更多的积分点，取决于非线性材料响应的程度。这些单元一般是应用于大规模的分析中，它们采用对角化质量矩阵，极为强健。高阶单元，诸如基于二次等参插值的，很快收敛到平滑的结果。但是，大多数大模型问题包括不平滑的现象，诸如弹-塑性和接触-碰撞，因此，在这些问题中，高阶单元的更高次近似的功能是不能实现的。,表9.2中列出了在计算软件中最频繁应用的单元，及它们的一些特点和缺陷。,9 一点积分单元,表9.2 4节点四边形壳单元,YASE单元（yet another shell element的缩称）,9 一点积分单元,最早的一点积分壳单元是BT单元。通过组合一个平板4节点单元和一个平面四边形4节点薄膜单元，构造了壳单元。当它出现翘曲构形时，不能够正确地作出反应（当应用单元的一条或者两条线模拟扭曲梁时，这个缺点就自己暴露了）。,HL单元，是基于CB壳理论。在显式程序中，它采用单一系列的积分点，因此需要沙漏控制；应用由Belytschko、Lin和Tsay发展的技术。在运算中它明显地慢于BT单元。,BWC单元修正了扭曲，即在BT单元中翘曲构形的缺陷。在BL单元中，引入了物理沙漏控制。这个沙漏控制是基于多场变分原理和Dvorkin-Bathe应变近似公式（见9.8.4-5，假设横向剪切速率）。在实际中，应变和应力状态的非均匀性限制了精确的物理沙漏控制。沙漏控制的这种形式提供了实质上的优点；它可以增加到适度大小的值而没有自锁，而在BT单元中，沙漏控制参数较高时将导致剪切自锁。,9 一点积分单元,BL单元和任何完全积分单元都被另一个缺陷困扰着。在具有大扭曲的问题中，这些单元会突然地和戏剧性地失效，并终止模拟。另一方面，BT单元在严重扭曲下是非常强健的，并且很少终止运算。在工业应用中，这具有很高的实用价值。因此，单一积分点单元的优点，不仅仅归于它们的高速度，在出现严重扭曲的问题中，它们趋于更加强健。诸如汽车碰撞模拟。,YASE单元改善了在梁弯曲时薄膜响应的薄膜场，即改进的弯曲运算。否则，它与BT单元是一致的。,BT、BWC和BL单元都是基于离散的Mindlin-Reissner理论；它们不是基于连续体单元。离散指这样的事实，即仅在积分点处，将Mindlin假设应用于运动。通过要求当前的法线保持直线，对运动施加约束。这可以看做是对Mindlin假设的另一种修正；不是要求初始法线保持直线，而是要求当前法线保持直线。,9 一点积分单元,速度场给出为,上波浪是指当前参考面的法线。关于运动的有限元近似为,将叉乘转换为矩阵相乘，上式可以写作,NI 为4节点等参形函数，,在积分点,关于转动方向的偏斜对称张量，定义在公式(9.5.42)中,处的转动变形率给出为,曲率,9 一点积分单元,在转动坐标系中计算薄膜应变和薄膜沙漏控制。在积分点处的曲率给出为,在一个任意的坐标系中，对于刚体转动，在曲率表达式中的最后一项不为零。而在转动坐标系中，刚体转动的节点速度正比于zyh，可以证明曲率为零，满足框架客观性。,其中，,9 一点积分单元,由于U1仅利用了一个系列的积分点，缺乏稳定性，单元是秩缺少的。单元有3个刚体位移模式：平动w，绕x和y的转动。有4个运动模式：3个沙漏模式和1个扭曲模式，3个沙漏模式是可以相互表示的，而1个平面内的扭曲模式是不能相互表示的。,其中，E 和G 为杨氏模量和剪切模量，A为单元的面积，在公式中定义的材料参数b，和 由用户自己设定，其范围必须在0.010.05之间。,由于U1仅利用了一个系列的积分点，单元是秩缺少的。对于一点积分，弯曲部分的秩是5：变形率场包含三个常数力矩和两个常数剪切。通过秩的分析，在横向剪切和曲率中，由于单元缺少线性项，所以，弯曲部分的秩缺乏为4。Hughes证明了伪奇异模式，3个沙漏模式是可以相互表示的(1个挠度和2个转角)，而1个平面内的扭曲模式是不能相互表示的。对3个相关模式应用沙漏控制：,9 一点积分单元,由于仅利用了一个系列的积分点，缺乏稳定性，单元是秩缺少的。Hughes证明了伪奇异模式。模式中的3个是可以相互表示的，而1个平面内的扭曲模式是不能相互表示的。,9 一点积分单元,9 一点积分单元,四边自由板的沙漏模式,9 一点积分单元,尚留下一个非传播的奇异模式扭曲模式,Ci是因数，相