工程壳有限元软件研发进展：复合材料壳单元求解器

关键词：壳单元，复合材料，有限元，铺层，求解器

在我为什么研发一个工程壳有限元软件一文中，我们实现了常规壳单元的研发、集成，并开发了界面和视口 交互功能。

自研工程壳有限元软件

我最初的目的就是要搞一套自己的复合材料壳有限元软件，这两天完成了第一版复合材料壳单元求解器的开发和初步验证，这里做一个情况总结。

壳单元

最顶级的工程力学思维，就是把工程问题，进行既简单又不失关键特征的简化。尤其在固体力学方面，前辈们建立了杆、梁、膜、板、壳等经典的力学模型，在实际工程中发挥着巨大的作用。

从有限元的角度看，壳单元实际上是由膜单元和板单元组合而成的单元，其中膜负责抗拉，板负责抗弯。编程的时候，需要分别完成膜和板的刚度矩阵，再把它们的刚度矩阵的子项填入壳单元刚度矩阵。

壳从几何上来说，是没有厚度的，厚度只是作为参数体现在刚度矩阵中，由此获得的比实体单元更少的单元和计算量。

最常用的壳单元有三角形和四边形两种。考虑到三角形划分网格的便利性，我从规划软件的时候，就定下来开发三维的三角形壳单元。

膜部分整体来说较为简单，它的刚度矩阵很多文献上都有现成结果，直接拿来用就行。麻烦的是三角形板单元，其几何矩阵非常复杂，有些文献推导出来的显式结果甚至是错的，我花了很大力气才完成这部分的开发。

复合材料壳单元

最常用的复合材料结构就是层合板，层合板的力学分析天然适合用壳单元来做。原因是：

（1） 实际的结构铺层多，用实体单元计算量太大；

（2） 壳没有实际的几何厚度，相应的复合材料的铺层角度、厚度等等信息都可作为参数体现在刚度矩阵，和复合材料力学理论内核一致；

（3）大部分的层合板结构厚度远远小于面内尺寸，适配壳的应用场景。

这里面的第二条尤为重要。根据复合材料力学，求解层合板刚度特性的基本思路是：

（1）认为单层板是横观各向同性的，并建立它的刚度矩阵；

（2）根据铺层角度，得到转换矩阵，计算当前层在整体坐标系下的刚度矩阵：

（3）根据每一层的刚度矩阵、厚度，积分得到层合板的载荷、变形关系：

以上是复合材料力学的内容，而这个过程也是我们开发复合材料壳单元的过程。和各向同性壳相比，复合材料壳除了拉伸（A矩阵）、弯曲（D矩阵）部分，还存在拉弯耦合（B）部分。

在编程的时候，就需要按照复合材料力学的过程，先完成A、B、D矩阵的编写，然后A传给膜部分，D传给板部分，最后组集壳的时候，把拉弯耦合B再加进去。

为什么ABAQUS铺层模块和UMAT无法同时使用

在完成复合材料壳单元的开发后，我较为充分理解了为什么ABAQUS铺层模块和UMAT无法同时使用。

UMAT的逻辑是，输入材料参数（单向），然后计算刚度矩阵，更新应力。铺层模块的逻辑是，定义铺层角度、厚度、材料（单向），自动计算出A、B、D矩阵，然后根据壳模型，得到层合壳的刚度矩阵。

我们自始至终输入的都是单向板的材料参数，而单元在几何上只有一层，因此在UAMT中，我们无法获得层合板的参数。

即便未来UMAT开发了获取层合后材料参数的接口，我们还得自己写一个膜、板、壳的代码，这已经不是自定义材料（UMAT），而是自定义单元技术了。

测试验证：对称铺层

材料参数：E1=150GPa；E2=10GPa；u12=0.3；G12=7.9GPa。

对称铺层：[0/90/90/0]

单元数量：2。

算例设置

位移、应力结果如下，从结果可以看出自研复合材料壳结果与商用软件结果一致。

自研求解器结果：位移

某商用软件结果：位移

自研求解器结果：应力S11包络

某商用软件结果：应力S11包络

测试验证：非对称铺层，拉弯耦合效应验证

非对称铺层：[0/90/0/90]

位移、应力结果如下，从结果可以看出自研复合材料壳求解器结果与商用软件结果一致。

非对称铺层存在拉弯耦合效应，从结果可以看出，在面内仅有拉伸载荷情况下，在面外方向出现了远大于面内方向的弯曲变形。自研复合材料壳求解器可以精准模拟出这种效应。

自研求解器结果：合位移

某商用软件结果：合位移

自研求解器结果：拉伸方向位移

某商用软件结果：拉伸方向位移

自研求解器结果：厚度方向位移

某商用软件结果：厚度方向位移

后续工作

后续我们将嵌入常用的强度准则（蔡吴、哈辛），并提供刚度折减定义，实现复合材料渐进损伤模拟。

然后将复合材料壳单元求解器，集成进我们的交互界面，并增加更复杂工况（离心力、压力、温度应力）的定义功能。进一步的，还可以和优化算法集 合，对工程实际结构进行铺层优化设计。