1、3.4 板壳单元板壳单元 在前面所处理的二维或三维体元中，应变位移关系是精确满足的，当然，最后的数值解是近似的而在板壳理论中对几何关系引入了一些假设， 使问题化为二维问题例如在前面已经介绍过的经典薄板理论中，采用了直法线的假设因此，薄板理论中的精确解只是在这些假设有效情况下而言1 在薄板理论中，位移未知数只有一个 广义应变和广义应力定义为：w yxwywxw/2/22222 xyyxMMM2 现在，由于应变是由挠度的二阶导数表示，因此，在有限元分析中，单元之间界面上的连续性条件要求形状函数和它的一阶法向导数必须在界面上连续，也即这种单元为 类单元，以前所介绍的连续体单元只要求位移在单元的界面上
2、连续，这种单元称为 类单元。1C0C3 由于 类单元的一阶导数连续的要求不容易做到因此，在薄板弯曲的有限元中，除了按位移求解的非协调单元外，杂交单元和混合单元也颇受重视，提出了许多种板壳单元。我们只介绍比较广泛采用的一种薄板单元和一种厚壳单元1C43.4.1 三角形薄板单元和薄壳单元三角形薄板单元和薄壳单元 用有限元法求解薄板弯曲问题时，用一 些离散的薄板单元代替原来的连续薄板，各单元只在节点互相连接，由于相邻单元之间要传递力矩，所以把节点当作刚性连接的每个节点有三个位移分量，即挠度 ，绕 轴的转角 和绕 轴的转角 wxxyy5在小变形的条件下，有几何关系：xwx/ywy/因此，节点的位移可表
3、示为：iiwxiTyiiwixw)(Tiyw)/(6 三角形薄板单元由于可以较好地适应复杂的边界形状，因此实用价值较大 如图 3.8 所示的三角形三节点单元中，每个节点有三个位移分量，三个节点共有 9 个节点位移，因此，所取的位移函数应包含 9 项，但一个完整的三次项包含10项如下：31029283726524321yxyyxxyxyxyx7图3.8 三角形薄板单元8 因此必须从上式中删去一项，式中前三项代表刚体位移，随后三个二次项代表常量变形（曲率），为保证解答的收敛性都不能删去，只能从后面的四个三次项中删去一项，如从中任意删去一项，使表达式失去对称性，如令 ，可保持公式的对称性，但在某些情
4、况下，例当三角形的二边平行于 轴 和 轴时，求解式中的系数将成为不可能如采用面积坐标，可解决这个矛盾98xy9 面积坐标的一次，二次，三次式分别有以下各项：一次式：mjiLLL,二次式：222,mjiimmjjiLLLLLLLLL三次式：333222222,mjiimmjjiimmjjimjiLLLLLLLLLLLLLLLLLL10 从面积坐标的公式可知， 的完整一次多项式可用面积坐标表示为：yx,mjiLLL321 的完整二次多项式至少应包括二式中共六项，例如可取为：yx,immjjimjiLLLLLLLLL65432111 的完整三次多项式应至少包括三次式中的四项和一次, 二次, 三次三式
5、中其他任取的六项，共十项，现取下列十项：yx,mjiimmjjiimmjjimjiLLLLLLLLLLLLLLLLLL,22222212 上列十项的组合可作为位移函数，但节点位移只有 9项，应设法减少一个独立项式中最后一项在三个节点都具有下列特点：0/ )(/ )(yLLLxLLLLLLmjimjimji13 现把它归并入其他三次项中，取位移项数为：)()(2524321mjimjmjijimjiLLcLLLLLcLLLLLLw)()(2726mjijimjiimLLcLLLLLcLLL)()(2928mjiimmjimjLLcLLLLLcLLL14经过一些代数运算后发现，为了满足薄板的变形条
6、件，上式中只能取 在式中代入相应的节点位移2/1ciwixixw)/(iyiyw)/(),(mji15 利用节点位移可求出式中的系数921,最后可得到位移函数如下：jjyiyixixiiiwNNNwNweymymxmxmmmyjyjxjxjNNNwNNN16iNNxiNyiNjNxjNyjNyjNxmNymNiewxiyijwxjyjmwymTym17)()(222mjimjiiiLLLLLLLN2/)()(2mjimjjmmjixiLLLbbLbLbLN2/)()(2mjimjjmmjiyiLLLccLcLcLNmjiyybjmixxc),(mji18经过适当整理后可得：miyijxijmm
7、jjiiLLcbLwLwLww)(mjymixmijiyjmxjmLLcbLLcb)()(19 上式前六项是完整的二次多项式，因此位移函数可满足关于单元刚体位移和常应变的条件另外，沿着单元的三条边界，挠度将是边长的三次函数，在每条边界的两个端点（节点）上各有两个参数（ 及 ），每条边界共有四个相同参数，因此可保证相邻两单元沿着公共边界上具有相同的位移wsw/20 在单元边界的法向斜率 也是边长的三次函数，但只在两个节点上有相同的斜率，只能部分地控制相邻单元的法向斜率，不能完全保证法向斜率的连续 从位移函数可以得到单元应变eBnw/21其中 是 的矩阵： B93 yxNyNxNBiii22222
8、2yxNyNxNymymym22222222单元刚度矩阵为：ATedxdyBDBk其中 为薄板弯曲理论中的弹性矩阵，对于各向同性材料的弹性矩阵为：01bDD012/ )1 (00D23其中 为弹性薄板的弯曲刚度：bD)1 (1223EhDb 为板的厚度。h 计算单元刚度矩阵时，利用面积坐标的积分公式可以得到刚度矩阵的解析表达式但公式很长，实际上在计算机中采用数值积分是很方便的由于被积函数是二次式，采用三点积分公式可得到较精确的结果24单元节点的等效载荷为：ATeqdxdyNp在均布载荷作用时可得到解析表达式： ieppjpTmp31qA24mjbb Tjicc24式中 为三角形的面积A25 上
9、述的三角形薄板单元是比较简单的单元但由于相邻单元的法向斜率不是完全协调的，只有计算网格比较规则时，计算结果能收敛于正确解答对于某些不规则网格，其计算精度就较差，为了克服这个缺点，有人建议对单元曲率进行修正，使之比较均匀也有人提出了各种混合单元，杂交单元等这里就不 一 一介绍了26 用有限元法分析弹性薄壳时，有两种不同的途径一条途径是用薄板单元组成的折板系统去代替原来的薄壳由平面应力状态和板弯曲应力状态加以组合而得到薄壳的应力状态另一途径是直接采曲面单元，根据壳体理论推导单元刚度矩阵由于每个单元的尺寸总是比较小的，因而有可能采用扁壳理论以简化单元刚度矩阵的计算27 弹性薄壳的应力状态可以认为是平
10、面应力状态和弯曲应力状态的组合因此，薄壳单元的刚度矩阵也可以由这两种应力状态的刚度矩阵加以组合而得到参阅图3.9. 把局部坐标系的轴和轴取在单元所在平面内 我们介绍用折板代替薄壳的计算方法28图3.9 薄壳单元的节点力和节点位移(a)平面应力状态 (b)弯曲应力状态29 对于平面应力状态，单元的应变状态完全取决于各节点的位移 ，以三角形单元为例，单元的节点力与节点位移的关系是：vu, pmjipmjipeFFFk30其中： iipivu iipiVUF式中上标 表示平面问题， 表示平面节点力节点的转角 不影响节点力，可以不考虑，相应的节点力 也不存在piiVU ,zmzjzi,zmzjziMM
11、M,31 对于弯曲应力状态，单元应变状态取决于节点在 方向的线位移 ，绕 轴的转角 及绕 轴的转角 ，节点位移和节点力的关系是：zwxxyy bmjibmjibeFFFk32其中： yixiibiw yixiibiMMWF式中上标 表示弯曲问题 ， 表示弯曲时的节点力和弯矩byixiiMMW,33 把平面应力和弯曲应力加以组合后，单元的节点位移和节点力是：iiuiviwxiyiTziiiUFiViWxiMyiMTziM 虽然转角 不影响单元的应力状态，为了便于以后把局部坐标系刚度矩阵转化为整体坐标系的刚度矩阵并集合成整体刚度矩阵，这里特地把 也包括在节点位移中，并在节点力中相应地包括一个虚拟弯
12、矩 ziziziM34 单元的节点位移和节点力的关系可写成：eeeFk式中 是在组合应力状态下的单元刚度矩阵由于平面应力状态下的节点力与弯曲应力状态下的节点位移互不影响，弯曲应力状态下的节点力与平面应力状态下的节点位移也互不影响，所以组合应力状态下的单元刚度矩阵的子矩阵可以写成如下形式：ek35 00002111pperskkk00002212ppkk000534333bbbkkk000544434bbbkkk000554535bbbkkk00000036式中引入标记 ppperskkk2111ppkk2212 bbbberskkkk534333bbbkkk544434bbbkkk554535
13、37 子矩阵 和 分别是平面应力问题和薄板弯曲问题的相应子矩阵因此，把平面应力问题和薄板弯曲问题中三角形刚度矩阵结合起来，即可得到组合应力状态下的三角形单元的刚度矩阵persk bersk 38 前面是在单元的局部坐标系中推导单元刚度矩阵的，即以单元的中面为 面， 轴垂直于单元的中面，为了建立薄壳的整体刚度矩阵，必须把不同平面内的单元刚度系数在节点处加以集合，因此，需要一个统一的整体坐标系，并把各单元在局部坐标系中的刚度矩阵转换到整体坐标系中去yx,z39 现在用 表示整体坐标，局部坐标仍用 表示在局部坐标系中，节点 的位移和节点力仍用 和 表示在整体坐标系中，节点处的位移和节点力用 和 表示
14、，即：zyx,zyx,iiiFtitiFtitiuiviwxiyiTtzititiUFiViWxiMyiMTtziM40图3.10 局部坐标系和整体坐标系41 节点位移和节点力在两个坐标系中的变换可按下式进行（参看图3.10）：tiitiiFF式中：42 000),cos(),cos(),cos(xzxyxx),cos(),cos(),cos(000zzzyzx而 为 轴与 轴的夹角，等等),(xxxx43 把整个单元在整体坐标中的节点位移和节点力分别用 , 表示，可得到：teteFteeLteeFLF式中： 00L000044由前面公式可知：1eteFLF 将局部坐标中单元节点力与节点位移的
15、关系式代入可得：11teeeeteLkLkLF45 由此可见此单元在整体坐标系中的刚度矩阵为：1LkLkete 因 是正交矩阵，故 ，由此可得：LTLL1LkLkeTte46 这就是刚度矩阵的变换公式。利用这个公式，可以由局部坐标系中的单元刚度矩阵 计算整体坐标系中的单元刚度矩阵 当整体坐标系中的单元刚度矩阵计算出来以后，就可以加以集合得到壳体的整体刚度矩阵 ekekK47 为了建立节点平衡方程，还须计算节点载荷令单元载荷所做的功等于节点载荷所做的功，不难求出节点载荷一般，在分布载荷作用下，单元节点载荷包括力和力矩分量 实际上，当单元比较小时，常可采用比较简单的方法计算单元节点载荷，即忽略力矩
16、分量，并把单元承受的外载荷的合力平均分摊到单元的各节点上去48 例如，当壳面承受着法向分布载荷 时，在三角形单元三个节点上的载荷强度分别是 ，那么作用于单元上的载荷合力是：pmjippp,3/ )(mjipppAR其中 是单元面积A49 假定载荷方向与局部坐标轴同时为正把这个合力平均分摊到三个节点上，并转换到整体坐标中去，得到在整体坐标中的节点载荷为：3),cos(),cos(),cos(RzzzyzxZYXiii其中 是整体坐标系zyx,50 把各单元的节点载荷加以集合，即可得到壳体的载荷列阵 壳体的节点平衡方程是：PttPK 从上式可解得各节点在整体坐标系中的位移利用变换式求出各节点在局部
17、坐标系中的位移，再利用局部坐标系中的应力矩阵，即可求得各单元的应力51 有一个特殊情况必须注意如果交会于 一个节点的各单元都在一个平面内，由于在前面已令 方向的刚度系数等于零，在局部坐标系中，这个节点的第六个平衡方程（相当于 方向）将是ziz00 52 如果整体坐标与这一局部坐标一致，显然，整体刚度矩阵的行列式 ，因而前面的平衡方程不满足有唯一解的条件如果整体坐标与局部坐标不一致，经变换后，在这个节点上得到了表面上正确的平衡方程，行列式 的行向量是线性相关的，从而仍有 的情况发生，所以平衡方程仍不满足有唯一解的条件0KK0K53 为了去掉这一困难，对于这种各有关单元位于同一平面内的节点，可以在
18、局部坐标系内建立节点平衡方程，并删去 方向的平衡方程，于是剩下的方程组满足有唯一解的条件，这个方法在程序设计上比较麻烦另一个办法是，在这些特殊节点上，给以任意的刚度系数 ，因此，在局部坐标系中，这个节点的平衡方程是：zzk0zizk54 经过坐标变换，在整体坐标中的节点平衡方程将满足有唯一解的条件在解出节点位移中包括 由于 不影响单元应力，并与其他节点平衡方程无关，所以实际上可以给定任意的 ，而不会影响计算结果zizizk55图3.11 三角形单元的局部坐标56 在前面所述中，为了建立节点平衡方程，需要用到各单元局部坐标的方向余弦矩阵 ， 下面介绍该矩阵的计算方法 如图3.11所示的三角形 ，
19、三个节点的整体坐标为 , 今取节点 为局部坐标的原点， 轴为 轴， 轴在单元平面内并垂直于 边， 轴垂直于单元平面ijm,iiizyx,jjjzyxmmmzyx,iijxyijz57 如果把边作为一个矢量 ，可用该矢量在整体坐标系中的三个分量表示：ijVjiijxVjiyTjiz 其中ijjixxxijjiyyyijjizzz58 用该矢量的长度 除它的三个分量，即得到该矢量的方向余弦，也就是局部坐标轴在整体坐标系中的方向余弦，从而得到一个单位矢量：xljijijixxzyxlzxyxxxv1),cos(),cos(),cos(592/1222)(jijijixzyxl 轴是垂直于三角形平面的
20、，利用两个矢量的矢积垂直于两个矢量所在平面的性质，可得到垂直于三角形平面的矢量 ： 式中zzV60mijimijimijimijimijimijiimijzxyyxzxxzyzzyVVV矢量 的长度：zV2/1222)()()(mijimijimijimijimijimijizxyyxzxxzyzzyl612022-5-1562故 轴的方向余弦是：zzzzlVzzyzxzv),cos(),cos(),cos(63 由于 轴垂直于 平面，因此单位矢量 和 的矢积等于 轴方向单位矢量：yxzxvzvyxzyvvzyyyxyv),cos(),cos(),cos( 这样，局部坐标的方向余弦矩阵的全部元
21、素都求出来了643.4.2 厚板单元和厚壳单元 在厚板和厚壳理论中，采用如下三 个基本假设：(一) 板的挠度是小的；(二) 变形前垂直于板中面的法线在变形后仍旧是直线，但不一定垂直于中面；(三) 忽略垂直于中面的法向应变65 这三个假设中的第一和第三两个假设是与 Kirchhoff的薄板理论中的假设相同，第二个假设有所不同，这个假设最早是Mindlin提出来的 我们介绍的厚板和厚壳单元是由三维体元退化而成的板壳单元在退化过程中利用了上述的假设二和三在这种单元中把板壳的中面位移和转角作为自 由度66图3.12 曲边厚壳单元67 现在先介绍厚板单元，考虑图3.12所示的厚板曲边单元，单元中面是平面
22、 , 是单元中面内的曲线坐标，而 是在厚度方向的直线坐标 及 代表单元的上、下表面单元中面内任一点的直角坐标可用曲线坐标表示如下：0z,11iixNx),(iiyNy),(68式中 为节点 的坐标， 为形函数，例如，对于图中所示8节点四边形曲边单元，形函数如下：iiyx ,iiN对于角点：) 1)(1)(1 (iiiiiN对于边中点：, 0, 0ii, 2/ )1)(1 (, 2/ )1)(1 (22iiiiNN69式中 为节点 的曲线坐标在用三角形单元时，则应采用面积坐标单元内任一点 的 坐标为：ii,i,z2/iitNz 形函数 在节点 处取值1，在其他节点处取值为0, 并可保证在单元边界
23、上坐标连续采用二次或三次形函数，可使单元边界适应任意的二次或三次曲线 一般用二次形函数iNi70图3.13 厚板单元的节点位移和节点力71 下面定义单元的位移场由于我们是分析厚板的弯曲问题，在中面内的线位移也予以忽略，根据厚板理论的基本假设，对于每个中面上的节点 ，我们只须确定 方向的线位移 和中面法线绕 轴的转角 和绕 轴的转角 ，如图3.13所示即：iziwyyixxi iiwyiTxi72 在节点 ，由于转角 ，离中面距离为 的一点在 方向的线位移为 ，由于转角 ，离中面距离为 点在 方向的线位移为 因此，在单元内任一点 的三个位移分量分别为：iyix2/yitxi2/xiit,y732
24、/ixiiNtv2/iyiiNtuiiwNw将上式对 求微商，得到：zyx,74yiiitxNxu)(21xiiitxNxv)(21iiwxNxw75yiiityNyu)(21iiwyNywxiiityNyv)(2176yiiitzNzu)(21xiiitzNzv)(21iiwzNzw77 因此，单元的应变可以写成： /Bzuxwywzvxvyuyvxuzxyzxyyx78 在厚板理论中忽略了厚度方向的正应变 ，因而也忽略厚度方向的正应力 ，其他应力由下式决定：zz xyxyyz DTzx其中 为应力应变矩阵：D79 0001)1 (2ED000100002100002121000080 式中
25、 的 和是材料的弹性模量和泊松比，最后两个剪应力项中所包含的系数 可取为1.20引入系数 是为了考虑剪应力分布不均匀的影响在前面的位移函数中，剪应力沿厚度方向接近均匀分布，实际上更接近于抛物线分布，数值 就是两种应变能的比值。E20. 181 有了 和 后就可以计算单元的刚度矩阵：BD 111111dddJBDBdxdydzBDBkTVTe 当然，单元刚度矩阵是转换到坐标 中用数值积分来计算的,82 现在将上述分析厚板的基本概念应用于厚壳的分析 采用厚壳曲面单元如图3.14所示单元的局部坐标为 为单元的中面， 为单元的上下表面，均为曲面， 和 是由直线产生的四个截面单元的节点取在中面上我们仍以
26、8结点为例.,011183图3.14 厚壳曲面单元84 厚板单元的中面是平面，而厚壳单元的中面是曲面设节点 的直角坐标为 ，则中面( )上任一点 的整体坐标 可表示如下(图3.15)：iiiizyx,0,iixNx),(iiyNy),(iizNz),(zyx,85图3.15 壳体中面86式中 为单元的形函数，八节点的形函数就是前面的形函数相同 单元内任一点的整体坐标可用局部坐标表示如下：),(iN87biiiizyxNzyx21),(tiiiizyxN21),(88式中下标 分别代表下表面和上表面令 为节点 在厚度方向的坐标差：tb,iiizyx,itiiibiiiiiizyxzyxzyx89
27、由前面两式可得到：iiiimiiiizyxNzyxNzyx2式中下标 代表中面，这样定义的局部坐标只是近似地垂直于中面m90 现在定义位移函数在每一个节点 有三个线位移和两个角位移为了定义角位移，如下图所示，作三个正交向量 和 首先作垂直于中面的向量iiiHH21,iH3iH3iixH3iyTiz 91 的方向余弦为：iH3iiitxl/3iiitym/3iiitzn/3其中 为节点 处的壳体厚度iti2/1)222(iiiizyxt92在节点 处再作两个向量 ，它们均与中面相切，并与 正交为了使 与 方向大体一致（这对计算结果的整理及边界条件处理带来一些方便），令：iiiHH21,iH3ii
28、iHHH321,zyx,iHHii32iiiHHH32193即 正交于 和 轴，则 与和 正 交根据向量运算法则，不难由上式求出 的方向余弦 和 的 方向余弦 图3.16中给出了节点 在 方向的线位移分别为 角位移如图3.17 所示, 为壳体中面法线向量绕 向量的转角， 为绕 向量的转角 iH2Hi 3xHi 1Hi2iH3Hi 1Hi2lmniii222,ix y z, ,u v wiii,iHi 3iHi2iiinml111,94图3.16 节点的线位移 95图3.17 节点的角位移 96 由于转角 ，在 方向产生的线位移为 ，它在 方向的位移分量为：i2/iittliii12/tmiii
29、12/tniii12/ 由于转角 ，在 方向产生的线位移为 ，它在 方向的位移分量为： iHi22/iitx y z, ,tliii22/tliii22/tniii22/zyx,iH197 利用形函数 ，单元内任一点的位移可用中面节点位移 表示： Ni( , ) i iiu viwiiiTuvwNuvwN tlmniiiii tiii2111lmniiiii22298 比较单元几何和单元位移两式可见，定义单元几何尺寸的式比定义单元位移的式具有较多的自由度，因此这种单元属于超参数单元，不会自动满足常应变条件 但从后面对应变的定义中可以看出，对于刚体位移和常应变的条件都是满足的在每个中面节点，共有
30、五个自由度 99 整体坐标 中的应变取决于下列矩阵： x y z, ,zuyuxu/zvyvxv/zwywxw/100 由位移函数式分别对 求微商，可得到上述矩阵中的各元素，例如： zyx,uxNxuMxt lliiiiiiii1212(),),(2121iiiiiiiilltxMuxNxv),(2121iiiiiiiilltxMuxNxw101式中 ，由偏微分法则可知： MNii 0/1iiiiiNNJzNyNxN102其中 为雅可比矩阵： J /xxxJ/yyy/zzz103同理可求得： MxMyMzJMMMJNNNiiiiiiiii11104 有了这些公式以后，前面矩阵中的元素均可以计算
31、, 因此, 单元的总体应变可以计算了在厚壳分析中，由于壳体中面的法线方向与整体坐标的轴 不一致，而且法线方向随着点的位置不同而变化，既然我们假设垂直于壳体中面的正应力等于零，就必须求出局部坐标系中的应力和应变今假定局部坐标系为 ，其中 轴正交于壳体中面 zzyx, z105图3.18 任意点的局部正交坐标系 106 前面已求出了节点 的局部正交坐标 现在建立单元内任意一点的局部正交坐标 如图3.18所示，在点 作曲面 ，在此曲面上有两条空间曲线 和 ，再作曲面 的两个切向量 和 ，其中 切于曲线 ， 切于曲线 ，由此可知： ixy z, 000,0000ddd0d0iiiHHH321,107k
32、ji,dixjykzddixjykzd()()式中 为 方向的单位向量今作 垂直于 和 ，则正交于 曲面： x y z, ,0z108/xxiddz/yyj/zzk z由此可求出 的方向余弦 l m n333,109 为使局部坐标 与整体坐标 大体一致，以利于计算结果的整理和边界条件的处理，对于其他两个局部坐标 和 用如下方法选取：xy z,zyx,xy yzxil31jm30kn30110 y zx由此可求出 的方向余弦 ，如 与 轴平行则改用 最后作 正交于 和 即lm n222, yzyxy zzyx由此可求出 的方向余弦 现得到局部坐标系的方向余弦矩阵如下： xl m n111, lm
33、n111lmn222lmn333111 利用矩阵 可将整体坐标系中的矩阵转换到局部坐标系中去： uxuyuzvxvyvz wxwywzuxuyuzTvxvyvz wxwywz112 由此得到局部坐标中的应变如下： xyxyyzzxuxvyuyvxvzwywxuzB113矩阵 可以分块， , B B1 iinnBBB212式中： iiuviwiiiT114 3113233212212211llllllllBimmmmmmmm1122122123323113115 1122122123323113nnnnnnnn1122122123323113 1122122123323113116ssisisi
34、ssisisiissisisissisisilNxmNynNzlMxmMynMztl lm mn nl lm mn n,(),2111222117 由广义虎克定律，局部坐标中的单元应力为： D 为弹性矩阵，以前已经给出在输出计算结果时，除局部坐标中的应力外，还可输出整体坐标中的应力： Dxzxyxxyyyz xzyzzxxyxzxyyyz Tzyzxz118 厚壳单元的刚度矩可直接在局部坐标系中用下式计算： 111111dddJBDBkTe设壳体中面承受了分布载荷 ，在结点 分布载荷在 方向的分量分别为 ，中面上任一点的分布载荷可表示为： ix y z, ,ziyixippp,xiixpNpyiiypNppN pziizp119参阅图3.18，在曲面 上取向量 和 所构成的微小面积 ，由向量运算法则可知： 0ddddddAd d 2)(yxyxA2/ 122)()(xzxzzyzy120 由虚功原理，可得到分布载荷 产生的节点载荷如下： pXNN pAd drpriix() 1111YNN pAd drpriiy() 1111ZNN pAd drpriiz() 1111MMrprp00,121 厚壳单元也能用于薄壳和薄板在求单元刚度矩阵时，对于八节点的二次单元，在 方向均取2个高斯点能取得很好的结果 , ,1222022-5-15123
