计算固体课件

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,3.4,板壳单元,在前面所处理的二维或三维体元中，应变位移关系是精确满足的，当然，最后的数值解是近似的而在板壳理论中对几何关系引入了一些假设,，,使问题化为二维问题例如在前面已经介绍过的经典薄板理论中，采用了直法线的假设因此，薄板理论中的精确解只是在这些假设有效情况下而言,1,3.4 板壳单元 在前面所处理的二维或三维,在薄板理论中，位移未知数只有一个 广义应变和广义应力定义为：,2,在薄板理论中，位移未知数只有一个 ,现在，由于应变是由挠度的二阶导数表示，因此，在有限元分析中，单元之间界面上的连续性条件要求形状函数和它的一阶法向导数必须在界面上连续，也即这种单元为 类单元，以前所介绍的连续体单元只要求位移在单元的界面上连续，这种单元称为 类单元。,3,现在，由于应变是由挠度的二阶导数表示，因,由于 类单元的一阶导数连续的要求不容易做到因此，在薄板弯曲的有限元中，除了按位移求解的非协调单元外，杂交单元和混合单元也颇受重视，提出了许多种板壳单元。我们只介绍比较广泛采用的一种薄板单元和一种厚壳单元,4,由于 类单元的一阶导数连续的要求不容,3.4.1,三角形薄板单元和薄壳单元,用有限元法求解薄板弯曲问题时，用一 些离散的薄板单元代替原来的连续薄板，各单元只在节点互相连接，由于相邻单元之间要传递力矩，所以把节点当作刚性连接的每个节点有三个位移分量，即挠度 ，绕 轴的转角 和绕 轴的转角 ,5,3.4.1 三角形薄板单元和薄壳单元 用,在小变形的条件下，有几何关系：,因此，节点的位移可表示为：,6,在小变形的条件下，有几何关系：因此，节点的位移可表示为：6,三角形薄板单元由于可以较好地适应复杂的边界形状，因此实用价值较大,如图,3.8,所示的,三角形三节点单元中，每个节点有三个位移分量，三个节点共有,9,个节点位移，因此，所取的位移函数应包含,9,项，但一个完整的三次项包含,10,项如下：,7,三角形薄板单元由于可以较好地适应复杂的边界,图,3.8,三角形薄板单元,8,图3.8 三角形薄板单元8,因此必须从上式中删去一项，式中前三项代表刚体位移，随后三个二次项代表常量变形（曲率），为保证解答的收敛性都不能删去，只能从后面的四个三次项中删去一项，如从中任意删去一项，使表达式失去对称性,，,如令 ，可保持公式的对称性，但在某些情况下，例当三角形的二边平行于 轴 和 轴时，求解式中的系数将成为不可能如采用面积坐标，可解决这个矛盾,9,因此必须从上式中删去一项，式中前三项代表,面积坐标的一次，二次，三次式分别有以下各项：,一次式：,二次式：,三次式：,10,面积坐标的一次，二次，三次式分别有以下各项,从面积坐标的公式可知，,的完整一次多项式可用面积坐标表示为：,的完整二次多项式至少应包括二式中共六项，例如可取为：,11,从面积坐标的公式可知， 的完整一次多项式可,的完整三次多项式应至少包括三次式中的四项和一次,二次,三次三式中其他任取的六项，共十项，现取下列十项：,12,的完整三次多项式应至,上列十项的组合可作为位移函数，但节点位移只有,9,项，应设法减少一个独立项式中最后一项在三个节点都具有下列特点：,13,上列十项的组合可作为位移函数，但节点位移只,现把它归并入其他三次项中，取位移项数为：,14,现把它归并入其他三次项中，取位移项数为：1,经过一些代数运算后发现，为了满足薄板的变形条件，上式中只能取 ,在式中代入相应的节点位移,15,经过一些代数运算后发现，为了满足薄板的变形条件，上式中只能,利用节点位移可求出式中的系数,最后可得到位移函数如下：,16,利用节点位移可求出式中的系数最后可得到位移,17,17,18,18,经过适当整理后可得：,19,经过适当整理后可得：19,上式前六项是完整的二次多项式，因此位移函数可满足关于单元刚体位移和常应变的条件另外，沿着单元的三条边界，挠度将是边长的三次函数，在每条边界的两个端点（节点）上各有两个参数（ 及 ），每条边界共有四个相同参数，因此可保证相邻两单元沿着公共边界上具有相同的位移,20,上式前六项是完整的二次多项式，因此位移函数,在单元边界的法向斜率 也是边长的三次函数，但只在两个节点上有相同的斜率，只能部分地控制相邻单元的法向斜率，不能完全保证法向斜率的连续,从位移函数可以得到单元应变,21,在单元边界的法向斜率,其中 是 的矩阵：,22,其中 是 的矩阵：22,单元刚度矩阵为：,其中 为薄板弯曲理论中的弹性矩阵，对于各向同性材料的弹性矩阵为：,23,单元刚度矩阵为：其中 为薄板弯曲理论中的弹性矩,其中 为弹性薄板的弯曲刚度：,为板的厚度。,计算单元刚度矩阵时，利用面积坐标的积分公式可以得到刚度矩阵的解析表达式但公式很长，实际上在计算机中采用数值积分是很方便的由于被积函数是二次式，采用三点积分公式可得到较精确的结果,24,其中 为弹性薄板的弯曲刚度： 为板的厚度,单元节点的等效载荷为：,在均布载荷作用时可得到解析表达式：,式中 为三角形的面积,25,单元节点的等效载荷为：在均布载荷作用时可得到解析表达式：式,上述的三角形薄板单元是比较简单的单元但由于相邻单元的法向斜率不是完全协调的，只有计算网格比较规则时，计算结果能收敛于正确解答对于某些不规则网格，其计算精度就较差，为了克服这个缺点，有人建议对单元曲率进行修正，使之比较均匀也有人提出了各种混合单元，杂交单元等这里就不 一 一介绍了,26,上述的三角形薄板单元是比较简单的单元但,用有限元法分析弹性薄壳时，有两种不同的途径一条途径是用薄板单元组成的折板系统去代替原来的薄壳由平面应力状态和板弯曲应力状态加以组合而得到薄壳的应力状态另一途径是直接采曲面单元，根据壳体理论推导单元刚度矩阵由于每个单元的尺寸总是比较小的，因而有可能采用扁壳理论以简化单元刚度矩阵的计算,27,用有限元法分析弹性薄壳时，有两种不同的途径,弹性薄壳的应力状态可以认为是平面应力状态和弯曲应力状态的组合因此，薄壳单元的刚度矩阵也可以由这两种应力状态的刚度矩阵加以组合而得到参阅图,3.9.,把局部坐标系的轴和轴取在单元所在平面内,我们介绍用折板代替薄壳的计算方法,28,弹性薄壳的应力状态可以认为是平面应力状态和,图,3.9,薄壳单元的节点力和节点位移,(a),平面应力状态,(b),弯曲应力状态,29,图3.9 薄壳单元的节点力和节点位移29,对于平面应力状态，单元的应变状态完全取决于各节点的位移 ，以三角形单元为例，单元的节点力与节点位移的关系是：,30,对于平面应力状态，单元的应变状态完全取决,其中：,式中上标 表示平面问题， 表示平面节点力,节点的转角 不影响节点力，可以不考虑，相应的节点力 也不存在,31,其中：式中上标 表示平面问题，,对于弯曲应力状态，单元应变状态取决于节点在 方向的线位移 ，绕 轴的转角 及绕 轴的转角 ，节点位移和节点力的关系是：,32,对于弯曲应力状态，单元应变状态取决于节点,其中：,式中上标 表示弯曲问题 ， 表示弯曲时的节点力和弯矩,33,其中：式中上标 表示弯曲问题 ，,把平面应力和弯曲应力加以组合后，单元的节点位移和节点力是：,虽然转角 不影响单元的应力状态，为了便于以后把局部坐标系刚度矩阵转化为整体坐标系的刚度矩阵并集合成整体刚度矩阵，这里特地把 也包括在节点位移中，并在节点力中相应地包括一个虚拟弯矩 ,34,把平面应力和弯曲应力加以组合后，单元的节点,单元的节点位移和节点力的关系可写成：,式中 是在组合应力状态下的单元刚度矩阵由于平面应力状态下的节点力与弯曲应力状态下的节点位移互不影响，弯曲应力状态下的节点力与平面应力状态下的节点位移也互不影响，所以组合应力状态下的单元刚度矩阵的子矩阵可以写成如下形式：,35,单元的节点位移和节点力的关系可写成：式中,36,36,式中,引入标记,37,式中引入标记37,子矩阵 和 分别是平面应力问题和薄板弯曲问题的相应子矩阵因此，把平面应力问题和薄板弯曲问题中三角形刚度矩阵结合起来，即可得到组合应力状态下的三角形单元的刚度矩阵,38,子矩阵 和,前面是在单元的局部坐标系中推导单元刚度矩阵的，即以单元的中面为 面， 轴垂直于单元的中面，为了建立薄壳的整体刚度矩阵，必须把不同平面内的单元刚度系数在节点处加以集合，因此，需要一个统一的整体坐标系，并把各单元在局部坐标系中的刚度矩阵转换到整体坐标系中去,39,前面是在单元的局部坐标系中推导单元刚度矩,现在用 表示整体坐标，局部坐标仍用 表示在局部坐标系中，节点 的位移和节点力仍用 和,表示在整体坐标系中，节点处的位移和节点力用 和 表示，即：,40,现在用 表示,图,3.10,局部坐标系和整体坐标系,41,图3.10 局部坐标系和整体坐标系41,节点位移和节点力在两个坐标系中的变换可按下式进行（参看图,3.10,）：,式中：,42,节点位移和节点力在两个坐标系中的变换可按下,而 为 轴与 轴的夹角，等等,43,而 为 轴与 轴的夹,把整个单元在整体坐标中的节点位移和节点力分别用,表示，可得到：,式中：,44,把整个单元在整体坐标中的节点位移和节点力分,由前面公式可知：,将局部坐标中单元节点力与节点位移的关系式代入可得：,45,由前面公式可知： 将局部坐标中单元节点力与节,由此可见此单元在整体坐标系中的刚度矩阵为：,因 是正交矩阵，故 ，由此可得：,46,由此可见此单元在整体坐标系中的刚度矩阵为：,这就是刚度矩阵的变换公式。利用这个公式，可以由局部坐标系中的单元刚度矩阵 计算整体坐标系中的单元刚度矩阵 当整体坐标系中的单元刚度矩阵计算出来以后，就可以加以集合得到壳体的整体刚度矩阵 ,47,这就是刚度矩阵的变换公式。利用这个公式,为了建立节点平衡方程，还须计算节点载荷令单元载荷所做的功等于节点载荷所做的功，不难求出节点载荷一般，在分布载荷作用下，单元节点载荷包括力和力矩分量 实际上，当单元比较小时，常可采用比较简单的方法计算单元节点载荷，即忽略力矩分量，并把单元承受的外载荷的合力平均分摊到单元的各节点上去,48,为了建立节点平衡方程，还须计算节点载荷,例如，当壳面承受着法向分布载荷,时，在三角形单元三个节点上的载荷强度分别是 ，那么作用于单元上的载荷合力是：,其中 是单元面积,49,例如，当壳面承受着法向分布载荷,假定载荷方向与局部坐标轴同时为正把这个合力平均分摊到三个节点上，并转换到整体坐标中去，得到在整体坐标中的节点载荷为：,其中 是整体坐标系,50,假定载荷方向与局部坐标轴同时为正把这个合,把各单元的节点载荷加以集合，即可得到壳体的载荷列阵 壳体的节点平衡方程是：,从上式可解得各节点在整体坐标系中的位移利用变换式求出各节点在局部坐标系中的位移，再利用局部坐标系中的应力矩阵，即可求得各单元的应力,51,把各单元的节点载荷加以集合，即可得到壳体的,有一个特殊情况必须注意如果交会于 一个节点的各单元都在一个平面内，由于在前面已令 方向的刚度系数等于零，在局部坐标系中，这个节点的第六个平衡方程（相当于 方向）将是,52,有一个特殊情况必须注意如果交会于 一个节,如果整体坐标与这一局部坐标一致，显然，整体刚度矩阵的行列式 ，因而前面的平衡方程不满足有唯一解的条件如果整体坐标与局部坐标不一致，经变换后，在这个节点上得到了表面上正确的平衡方程，行列式,的行向量是线性相关的，从而仍有,的情况发生，所以平衡方程仍不满足有唯一解的条件,53,如果整体坐标与这一局部坐标一致，显然，整体刚,为了去掉这一困难，对于这种各有关单元位于同一平面内的节点，可以在局部坐标系内建立节点平衡方程，并删去 方向的平衡方程，于是剩下的方程组满足有唯一解的条件，这个方法在程序设计上比较麻烦另一个办法是，在这些特殊节点上，给以任意的刚度系数 ，因此，在局部坐标系中，这个节点的平衡方程是：,54,为了去掉这一困难，对于这种各有关单元位于同,经过坐标变换，在整体坐标中的节点平衡方程将满足有唯一解的条件在解出节点位移中包括 由于,不影响单元应力，并与其他节点平衡方程无关，所以实际上可以给定任意的 ，而不会影响计算结果,55,经过坐标变换，在整体坐标中的节点平衡方程,图,3.11,三角形单元的局部坐标,56,图3.11 三角形单元的局部坐标56,在前面所述中，为了建立节点平衡方程，需要用到各单元局部坐标的方向余弦矩阵 ， 下面介绍该矩阵的计算方法,如图,3.11,所示的三角形 ，三个节点的整体坐标为,今取节点 为局部坐标的原点， 轴为 轴， 轴在单元平面内并垂直于 边， 轴垂直于单元平面,57,在前面所述中，为了建立节点平衡方程，需要用,如果把边作为一个矢量 ，可用该矢量在整体坐标系中的三个分量表示：,其中,58,如果把边作为一个矢量 ，可用该矢量在,用该矢量的长度 除它的三个分量，即得到该矢量的方向余弦，也就是局部坐标轴在整体坐标系中的方向余弦，从而得到一个单位矢量：,59,用该矢量的长度 除它的三个分量，即得,轴是垂直于三角形平面的，利用两个矢量的矢积垂直于两个矢量所在平面的性质，可得到垂直于三角形平面的矢量 ：,式中,60,轴是垂直于三角形平面的，利用两个矢,矢量 的长度：,61,矢量 的长度：61,2024/8/28,62,2023/9/462,故 轴的方向余弦是：,63,故 轴的方向余弦是：63,由于 轴垂直于 平面，因此单位矢量 和 的矢积等于 轴方向单位矢量：,这样，局部坐标的方向余弦矩阵的全部元素都求出来了,64,由于 轴垂直于 平面，因此,3.4.2,厚板单元和厚壳单元,在厚板和厚壳理论中，采用如下三 个基本假设：,(,一,),板的挠度是小的；,(,二,),变形前垂直于板中面的法线在变形后仍旧是直线，但不一定垂直于中面；,(,三,),忽略垂直于中面的法向应变,65,3.4.2 厚板单元和厚壳单元 在厚板和厚,这三个假设中的第一和第三两个假设是与,Kirchhoff,的薄板理论中的假设相同，第二个假设有所不同，这个假设最早是,Mindlin,提出来的,我们介绍的厚板和厚壳单元是由三维体元退化而成的板壳单元在退化过程中利用了上述的假设二和三在这种单元中把板壳的中面位移和转角作为自 由度,66,这三个假设中的第一和第三两个假设是与 K,图,3.12,曲边厚壳单元,67,图3.12 曲边厚壳单元67,现在先介绍厚板单元，考虑图,3.12,所示的厚板曲边单元，单元中面是平面,是单元中面内的曲线坐标，而 是在厚度方向的直线坐标,及 代表单元的,上、下表面单元中面内任一点的直角坐标可用曲线坐标表示如下：,68,现在先介绍厚板单元，考虑图3.12所示的,式中 为节点 的坐标， 为形函数，例如，对于图中所示,8,节点四边形曲边单元，形函数如下：,对于角点：,对于边中点：,69,式中 为节点 的坐标，,式中 为节点 的曲线坐标在用三角形单元时，则应采用面积坐标单元内任一点 的 坐标为：,形函数 在节点 处取值,1,，在其他节点处取值为,0,并可保证在单元边界上坐标连续采用二次或三次形函数，可使单元边界适应任意的二次或三次曲线 一般用二次形函数,70,式中 为节点 的曲线坐标在用三角形,图,3.13,厚板单元的节点位移和节点力,71,图3.13 厚板单元的节点位移和节点力71,下面定义单元的位移场由于我们是分析厚板的弯曲问题，在中面内的线位移也予以忽略，根据厚板理论的基本假设，对于每个中面上的节点 ，我们只须确定 方向的线位移 和中面法线绕 轴的转角 和绕 轴的转角 ，如图,3.13,所示即：,72,下面定义单元的位移场由于我们是分析厚板,在节点 ，由于转角 ，离中面距离为 的一点在 方向的线位移为 ，由于转角 ，离中面距离为 点在 方向的线位移为 因此，在单元内任一点 的三个位移分量分别为：,73,73,将上式对 求微商，得到：,74,将上式对 求微商，得到：74,75,75,76,76,77,77,因此，单元的应变可以写成：,78,因此，单元的应变可以写成：78,在厚板理论中忽略了厚度方向的正应变 ，因而也忽略厚度方向的正应力 ，其他应力由下式决定：,其中 为应力应变矩阵：,79,在厚板理论中忽略了厚度方向的正应变,80,80,式中 的 和是材料的弹性模量和泊松比，最后两个剪应力项中所包含的系数 可取为,1.20,引入系数 是为了考虑剪应力分布不均匀的影响在前面的位移函数中，剪应力沿厚度方向接近均匀分布，实际上更接近于抛物线分布，数值 就是两种应变能的比值。,81,式中 的 和是材料的弹性模量和,有了 和 后就可以计算单元的刚度矩阵：,当然，单元刚度矩阵是转换到坐标 中用数值积分来计算的,82,有了 和 后就可以计,现在将上述分析厚板的基本概念应用于厚壳的分析,采用厚壳曲面单元如图,3.14,所示单元的局部坐标为 为单元的中面， 为单元的上下表面，均为曲面， 和 是由直线产生的四个截面单元的节点取在中面上我们仍以,8,结点为例,83,现在将上述分析厚板的基本概念应用于厚壳的分,图,3.14,厚壳曲面单元,84,图3.14 厚壳曲面单元84,厚板单元的中面是平面，而厚壳单元的中面是曲面设节点 的直角坐标为 ，则中面,( ),上任一点,的整体坐标 可表示如下,(,图,3.15),：,85,厚板单元的中面是平面，而厚壳单元的中面是曲,图,3.15,壳体中面,86,图3.15 壳体中面86,式中 为单元的形函数，八节点的形函数就是前面的形函数相同,单元内任一点的整体坐标可用局部坐标表示如下：,87,式中 为单元的形函数，八节点的形,88,88,式中下标 分别代表下表面和上表面令 为节点 在厚度方向的坐标差：,89,式中下标 分别代表下表面和上表面令,由前面两式可得到：,式中下标 代表中面，,这样定义的局部坐标只是近似地垂直于中面,90,由前面两式可得到：式中下标 代表中面，这样定义的局部坐标,现在定义位移函数在每一个节点 有三个线位移和两个角位移为了定义角位移，如下图所示，作三个正交向量 和 首先作垂直于中面的向量,91,现在定义位移函数在每一个节点 有,的方向余弦为：,其中 为节点 处的壳体厚度,92,的方向余弦为：其中 为节点 处,在节点 处再作两个向量 ，它们均与中面相切，并与 正交为了使 与 方向大体一致（这对计算结果的整理及边界条件处理带来一些方便），令：,93,在节点 处再作两个向量 ，,即 正交于 和 轴，则 与和 正 交根据向量运算法则，不难由上式求出 的方向余弦 和 的 方向余弦,图,3.16,中给出了节点 在 方向的线位移分别为 ,角位移如图,3.17,所示,为壳体中面法线向量绕 向量的转角， 为绕 向量的转角,94,即 正交于 和 轴，则 与和 正 交根,图,3.16,节点的线位移,95,图3.16 节点的线位移 95,图,3.17,节点的角位移,96,图3.17 节点的角位移 96,由于转角 ，在 方向产生的线位移为,，它在 方向的位移分量为：,由于转角 ，在 方向产生的线位移为,，它在 方向的位移分量为：,97,由于转角 ，在 方向产生的线位移为,利用形函数 ，单元内任一点的位移可用中面节点位移 表示：,98,利用形函数 ，单元内任一点的位移可用中,比较单元几何和单元位移两式可见，定义单元几何尺寸的式比定义单元位移的式具有较多的自由度，因此这种单元属于超参数单元，不会自动满足常应变条件,但从后面对应变的定义中可以看出，对于刚体位移和常应变的条件都是满足的在每个中面节点，共有五个自由度,99,比较单元几何和单元位移两式可见，定义单元几何尺寸的,整体坐标 中的应变取决于下列矩阵：,100,整体坐标 中的应变取决于下列矩阵： 10,由位移函数式分别对 求微商，可得到上述矩阵中的各元素，例如：,101,由位移函数式分别对 求微商，可得到上述矩,式中 ，由偏微分法则可知：,102,式中 ，由偏微分法则可知： 102,其中 为雅可比矩阵：,103,其中 为雅可比矩阵： 103,同理可求得：,104,同理可求得： 104,有了这些公式以后，前面矩阵中的元素均可以计算,因此,单元的总体应变可以计算了在厚壳分析中，由于壳体中面的法线方向与整体坐标的轴 不一致，而且法线方向随着点的位置不同而变化，既然我们假设垂直于壳体中面的正应力等于零，就必须求出局部坐标系中的应力和应变今假定局部坐标系为 ，其中 轴正交于壳体中面,105,有了这些公式以后，前面矩阵中的元素均可以计算,图,3.18,任意点的局部正交坐标系,106,图3.18 任意点的局部正交坐标系 106,前面已求出了节点 的局部正交坐标,现在建立单元内任意一点的局部正交坐标 ,如图,3.18,所示，在点 作曲面,，在此曲面上有两条空间曲线 和 ，再作曲面 的两个切向量 和 ，其中 切于曲线 ， 切于曲线 ，由此可知：,107,前面已求出了节点 的局部正交坐标,式中 为 方向的单位向量今作 垂直于 和 ，则正交于 曲面：,108,式中 为 方向的单位向量今作 垂直于,由此可求出 的方向余弦 ,109,由此可求出 的方向余弦 109,为使局部坐标 与整体坐标 大体一致，以利于计算结果的整理和边界条件的处理，对于其他两个局部坐标 和 用如下方法选取：,110,为使局部坐标 与整体,由此可求出 的方向余弦 ，如 与 轴平行则改用 最后作 正交于 和 即,由此可求出 的方向余弦 现得到局部坐标系的方向余弦矩阵如下：,111,由此可求出 的方向余弦 ，如,利用矩阵 可将整体坐标系中的矩阵转换到局部坐标系中去：,112,利用矩阵 可将整体坐标系中的矩阵转换到局部坐,由此得到局部坐标中的应变如下：,113,由此得到局部坐标中的应变如下： 113,矩阵 可以分块，,式中：,114,矩阵 可以分块，式中： 114,115,115,116,116,117,117,由广义虎克定律，局部坐标中的单元应力为：,为弹性矩阵，以前已经给出在输出计算结果时，除局部坐标中的应力外，还可输出整体坐标中的应力：,118,由广义虎克定律，局部坐标中的单元应力为：,厚壳单元的刚度矩可直接在局部坐标系中用下式计算：,设壳体中面承受了分布载荷 ，在结点 分布载荷在 方向的分量分别为 ，中面上任一点的分布载荷可表示为：,119,厚壳单元的刚度矩可直接在局部坐标系中用下式计算：,参阅图,3.18,，在曲面 上取向量 和 所构成的微小面积 ，由向量运算法则可知：,120,参阅图3.18，在曲面 上取向量 和 所构成,由虚功原理，可得到分布载荷 产生的节点载荷如下：,121,由虚功原理，可得到分布载荷 产生的节点载荷如下：,厚壳单元也能用于薄壳和薄板在求单元刚度矩阵时，对于八节点的二次单元，在,方向均取,2,个高斯点能取得很好的结果,122,厚壳单元也能用于薄壳和薄板在求单元刚度矩阵时，,2024/8/28,123,2023/9/4123,