有限元理论与Abaqus实现研究：壳的应力方向探讨

==概述==

本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式。有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。

一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元求解器，通过自研求解器和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。

iSolver介绍视频：

http://www.jishulink.com/college/video/c12884

==壳的应力方向 ==

有限元中，物理量用的最多的是标量、矢量和二阶张量。其中位移、坐标等都是矢量，而应变、应力等都是二阶张量。矢量很容易理解，体的应力等二阶张量直接就采用了全局坐标系的也不会有方向理解问题，但梁壳的应力结果很容易搞错，后处理结果中的S11、S12等的方向有时会觉得和预想的不一致但又不明所以。同时，这个方向也是单元材料的方向，在自编程序时，如果一开始坐标系的定义就弄错了，那么将直接导致和材料相关的刚度矩阵的错误，所以弄清应力的方向定义对自编程序和理解有限元结果都相当重要。

本章将简单介绍一下数学上张量和Abaqus中壳的应力方向，并说明Abaqus这么选取的意义，最后通过自编程序iSolver来验证壳的应力方向的正确性。

1.1 数学上的张量方向

矢量的方向是一定的，但它的分量都是基于某个坐标系定义的，坐标系不同，那么分量结果也会不同。

矢量可以表示为：                    

显然，分量和坐标系的选取有关。譬如我们一般的直角全局坐标系如下，那么分量就是普通的x、y、z三个分量值。

和矢量类似，二阶张量可以表示如下，当然也可以用一个更简单的3X3的矩阵表示，显然，二阶张量的分量等也与坐标系的取值有关。

1.2 Abaqus壳的应力方向

Abaqus后处理中壳的应力会输出S11，S22，S12等分量，分别对应上面二阶张量a的a11、a22、a12等分量，其它分量不输出，这三个量与壳的坐标系的选取密切相关。Abaqus中，在《anlysis user manual 12.1》的29.6.7 Three-dimensional conventional shell element library的Element output说明了壳采用的是局部坐标系，具体定义如下：

S11：Local 11 direct stress.

S22：Local 22 direct stress.

S12：Local 12 shear stress.

1）   local 1（以下称为T1方向）默认情况下为x方向在表面上的投影。

如果x方向正好垂直于表面(Abaqus取cos夹角<0.1度）,那么取z方向投影，所以下方的圆柱面上当x轴和表面垂直时，S11的方向明显和上下不一致了。

2）n方向（即T3方向）垂直于壳的表面，至于是往上还是往下，由节点顺序决定。

3）第三个方向local 2（以下称为T2方向）和T1，T3满足右手定则。显然在壳的表面绕T3转动90度。

由上面的定义，可以看出应力的方向的S11，S22都是沿着表面的，而不是全局坐标系下的。

1.3 应力方向选取原因

壳的应力方向为何要取成上面的T1，T2和T3？有两个原因：

1.3.1 插值函数

T1，T2在面内才能使得坐标和位移可以使用二维平面内的插值公式。

位移的插值是有限元的基础，插值方法如下：

壳的插值方式有两种，一种是平面壳的理论，先按平面壳来计算K，然后再将K做坐标变换，此时单元内部任意点的坐标值只有x，y，没有z，插值为四个节点的坐标。这种方法就要求在壳的坐标系下，四个节点的z方向的坐标zi=0，只剩xi和yi，这种方法计算精度不是很准，因为一个四边形的四个点不能总是保证在一个平面上的，只要是曲面划分成四边形很有肯定就会是这种情况。另一种更精确的算法是当做曲面壳，abaqus就是这么做的。此时坐标的插值函数不变，但换为了三个坐标x、y、z的分别插值，四个节点的z方向的坐标zi不需要强制要求是0了。

不管是哪种方法，插值函数都是二维的，如下：

只有T1、T2在壳平面内实际坐标才能映射到等参上。

1.3.2 本构关系

壳的应变方向为何要取成T3和面垂直的另一个原因是因为只有这样用壳来简化分析对象时材料的本构关系才是最简单的。因为壳是对体的简化，当体的厚度远小于面内尺寸时，那么可以用壳的理论来近似，此时可以把壳的刚度分为薄膜效应刚度、面外弯曲刚度、面外横向剪切刚度等部分，每一部分本构关系由相应简单的材料相关的D矩阵决定，譬如薄膜效应和面外弯曲的D矩阵都是下方矩阵：

上面的本构关系是各向同性材料的，面内的T1和T2不影响D矩阵，但各向异性就不一样，此时D矩阵将与面内的T1，T2相关，这也是应力的坐标系也叫做材料坐标系的原因。

1.4 程序内部实现方法

其实，原理非常简单，但程序实现并不像看起来的那么简单。无论是Abaqus还是自编程序内部想要实现壳的应力的方向，和前面讨论的应力方向选取原因是一致的，主要是两点：

（1）坐标和位移用局部坐标系表示，也就是下面等式：

S1，S2，S3分别对应T1，T2，T3三个方向的坐标分量。

（2）本构关系用局部坐标系做旋转变化。

1.5 iSolver验证

为了证明壳的S11、S22是沿表面方向而不是全局坐标系方向，我们用一个球壳加均匀压力，可发现S11和S22在abaqus中是球对称的，只有沿着表面才会是应力相等，如果沿着全局的xyz，那么显然不应该球对称。

为了进一步验证，我们在iSolver中按照Abaqus的T1、T2、T3定义壳单元方向，并通过一个简单的例子来验证Abaqus壳的应力方向。

1.5.1 算例

只取一个简单的长方形单元。

参数如下：

尺寸：5X1，厚度0.1。

材料：Young’s Modulus 1e8， Poisson Ratio 0.3。

右侧两个节点固支。

左侧两个节点每个加集中力1e5，在壳的平面内往外拉伸。

一开始在XY平面内，长度5方向在X轴，宽度方向在Y轴上。此时将得到一个应力的二阶张量S。

1.5.2 Y轴转动45度

沿Y轴转动45度，按照Abaqus的定义方式，显然，二阶张量S完全不变。

1.5.3 Z轴转动45度

沿Z轴转动45时，显然，壳的应力将不同，需要坐标变换。

设R为从原始模型的到新模型的变换矩阵，那么从数学的角度如果是一个列矢量a，那么变换后的矢量为：

A=R*a

而如果是二阶张量S，在原始坐标系下二阶张量可以表示为

S=a*b’

b’为列矢量b的转置。

显然变换后的张量为：

SNew=(R*a)*(R*b)’=R*a*b’*R’=R*S*R’

其中R’表示转置。

1.5.4 iSolver验证演示

iSolver中详细的验证证明请参考演示视频（带配音）：

https://www.jishulink.com/college/video/c12884 20理论系列文章14：壳的应力方向

==总结==

本章简单介绍了一下数学上张量和Abaqus中壳的应力方向，并说明Abaqus这么选取的意义，最后通过自编程序iSolver来验证壳的应力方向的正确性。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删