基于ABAQUS的固体火箭发动机复合材料壳体快速化建模方法及验证分析①

祖 磊，葛 庆，李德宝，张桂明，吴乔国，张 骞

(1.合肥工业大学 机械工程学院，合肥 230009；2.合肥工业大学 土木与水利工程学院，合肥 230009)

0 引言

碳纤维增强树脂基复合材料(CFRP)具有高比模量、高比强度、优异的耐疲劳耐腐蚀性等诸多优异特性。其中，纤维缠绕压力容器作为重要的复合材料能源贮存装置，在固体火箭、飞行器、海洋船舶、汽车等各个领域有着广泛的应用。国内外针对纤维缠绕压力容器研究已有较完备的分析设计流程，包括线型设计、失效分析、疲劳分析和爆压预测等。

近年来，海内外大量学者针对复合材料压力容器进行了大量的研究，JUAN通过纤维缠绕压力容器插件WCM对复合材料压力容器的缠绕层进行了精确建模，并对设计压力为70 MPa的高压储氢Ⅳ型气瓶进行有限元仿真和爆破实验验证。JOHANSEN等设计了纤维缠绕分析程序，通过CAD/CAE/CAM一体化设计方法，实现了对任意的轴对称回转体及其组合体进行缠绕分析。矫维成基于封头筒身的缠绕纱带总体积相等的原理，提出了一种基于三次样条函数预测复合材料压力容器封头厚度的方法。许家忠等通过多项式逼近算法实现压力容器缠绕层封头的厚度分布预测，解决了传统厚度预测公式在封头极孔壁厚预测值和实际值存在较大误差的问题。ZU等通过微分理论推导出了压力容器的非测地线方程，研究极孔切点数和纱线宽度对非测地线缠绕方式的影响。王华毕等基于网格理论设计缠绕层，通过三维扫描仪测量压力容器轮廓验证厚度预测，然后进行爆压预测和失效分析。

目前，针对复合材料压力容器的铺层设计方法主要包括网格理论和有限元方法。其中，建立准确的压力容器模型是有限元分析准确的前提。由于碳纤维复合材料压力容器在封头处的缠绕层厚度为连续变厚度、变角度。导致压力容器有限元模型建模难度大、建模工作繁琐。针对缠绕层数较多的模型，考虑时间成本只能对模型简化处理，影响最终的有限元计算精度。

国内外学者在纤维缠绕压力容器的设计及建模分析等方面进行了大量的研究，但针对壳体的快速化建模方法研究较少。综上，本文在前人研究的基础上，提出了一种固体火箭发动机复合材料壳体快速建模的方法。实现对商业缠绕软件信息文件进行数据处理，导入到ABAQUS中参数化生成有限元模型的功能。通过复合材料壳体的水压爆破实验验证，与生成的有限元模型计算结果进行对比，验证了建模方法的准确性。这种复合材料壳体快速化建模方法提高了建模效率，能较快速地实现从缠绕设计软件导入生成有限元几何模型，同时可以便捷地生成缠绕层封头补强，缩短了设计周期，该研究成果对壳体的设计分析有着重要的意义。

1 复合材料壳体快速建模方法

1.1 快速化建模流程

复合材料壳体有限元模型的准确性主要是由封头上纤维厚度和缠绕线型所决定的。传统的壳体设计及分析流程中，需要手动在ABAQUS中建立有限元模型。因壳体在封头处的角度和厚度是连续变化的，导致实际建模的过程中封头处的建模难度大，需要花费大量时间。考虑时间成本问题，人工对封头进行切分赋予角度时，无法保证铺层角度的连续性，影响模型计算精度。

本文基于ABAQUS二次开发功能，实现了从缠绕软件中快速导入的参数化建模功能，并可以实现包括缠绕补强、封头补强、螺旋缠绕等多种铺层方式的建模。可以将几天的建模工作缩短至几分钟，可以较大程度地缩短设计人员对有限元模型的前处理时间，建模时间成本降低，方便壳体的结构设计及优化分析。本插件除了实现标准椭圆封头结构的壳体建模功能，也可以实现如固体火箭发动机尾喷管等异形件的非标准结构建模功能。

缠绕软件设计线型建立模型后，可生成缠绕层信息文件，其中包括缠绕层落纱点坐标、落纱点缠绕角度和落纱点位置纱线厚度。图1为CADWIND软件的缠绕层信息文件部分截图。本文将纤维缠绕层数据离散化，通过PYTHON语言对缠绕软件生成的信息文件进行数据处理，在ABAQUS中参数化建模生成有限元模型，实现复合材料壳体缠绕层的快速化建模。

图1 数据文件Fig.1 Data file

基于缠绕层信息文件和PYTHON语言，编写了数据处理程序。信息文件中的数据一般根据纤维实际缠绕的位置生成对应参数，缠绕的坐标变化顺序是跟着丝嘴轨迹而变化的，故信息文件中缠绕坐标不是线性变化的。数据处理程序先读取数据文件中的落纱点坐标，根据坐标相对位置进行排序，剔除重复数据，绘制出缠绕层内轮廓。然后，读取各落纱点位置上的纱线厚度，绘制出缠绕层外轮廓，依次往复完成壳体有限元模型建立。

部分缠绕软件如CADWIND基于面片法原理，会根据芯模轮廓打断成直线进行求解建模，曲线打断的直线疏密会影响缠绕层信息文件，采样的数据往往会产生噪点，甚至会偏离模型轮廓，如图2所示。导致生成的有限元模型表面会产生褶皱，不利于模型的分析。这里采取B样条函数法对缠绕层外轮廓数据进行光滑处理，保证生成的模型表面光滑连续。

图2 原始数据Fig.2 Original data

由于部分商用缠绕软件只能生成螺旋层缠绕数据信息。为满足更通用性的建模工作，本方法中增加了环向层及补强层建模功能。在缠绕层内部中出现了环向层或补强层时，可以根据相邻螺旋缠绕层位置判断生成环向层和补强层数据。如果缠绕层最外层为环向层或补强层时，信息文件只能获取螺旋层位置，无法求解环向层或补强层位置，需要在插件界面输入最外层铺层信息，辅助生成对应铺层的有限元模型。

基于ABAQUS二次开发方法，编写了缠绕软件导入的快速化壳体建模脚本程序。将数据处理程序中采集到的数据轮廓点导入，建立几何模型。筒身段，根据缠绕角进行材料属性的赋予。封头段，通过切片法对封头模型沿母线方向进行切分，将数据处理程序中得到的缠绕角信息依次对封头进行材料属性的赋予，通过柱坐标系离散化的方法完成封头模型的材料属性赋予。然后，根据需求对壳体模型施加对应的边界条件、工作载荷以及最终的网格划分，完成整个壳体模型的前处理工作。

1.2 缠绕软件导入数据的处理方法

为解决数据处理程序中波动的外轮廓问题，本文采用了数值分析中的B样条方法，对实际的波动数据进行处理。样条线的概念是过一组指定点集而生成的柔性光滑曲线，而B 样条曲线法的原理即通过控制点集中局部位置，从而达到控制一条曲线形状的样条函数。

B样条定义方法较多，一般采取截尾幂函数差商定义，这里采用德布尔-考克斯递推公式。B样条曲线的方程为

(1)

其中，为控制顶点，，()为基函数，而基函数表示如下：

，()的第一个下标表示序号，第二个下标表示阶次。如确定第个次B样条，()，需要确定，+1，…，++1一共+2个节点。这里把[，++1]定义为，()的支撑区间。

B样条光滑曲线原理如图3所示。

图3 B样条数据处理Fig.3 B-spline data processing

当两节点为同一位置时，递推公式中组合系数在特殊情况时会出现分子分母均为0的现象，这种情况直接令系数为0即可。

该曲线在区间[，++1]中某一点()是由+1个顶点(=-，-+1，…，)和基函数决定，与其他顶点无关，如图4所示。

图4 B样条顶点控制原理Fig.4 B-spline vertices approach

由图4可见，移动第个控制顶点，最多只会对区间(，++1)上的局部曲线产生影响，而其他区间的B样条曲线不会受到干预。因此，阶次越高，B样条函数会拟合出的数据点越光滑，但递推函数计算量就越大。

1.3 可视化界面

基于PYTHON语言和ABAQUS二次开发功能，为建模方法设计了可视化界面。可以在插件中界面输入相关信息进行模型的前处理工作，方便设计人员操作，界面主要分为两部分，如图5所示。

(a)Interface A (b)Interface B图5 插件界面Fig.5 Plug-in interface

针对任意形状的构件，可以在图5(a)中输入参数实现模型生成。在文件指派模块中，选择缠绕软件生成的缠绕层信息文件，然后在缠绕层参数指派模块中，输入设计爆压、环向单层厚度、三维模型旋转角度以及筒身长度。在设计参数指派模块中，输入工程常数和缠绕层形式，其中主要缠绕形式为螺旋缠绕、环向缠绕等。在数据处理模块中，可以将每次使用到的参数进行保存。再次使用时，相同的参数选中文件可以直接导入，节约输入参数的时间。

针对任意比例的椭圆壳体，可在图5(b)中实现参数化建模。首先，在ABAQUS中建立芯模模型；之后，在参数栏中输入相关缠绕和建模参数，实现包括螺旋缠绕、环向缠绕、扩孔和补强等多种缠绕形式的壳体有限元模型建立。

2 壳体缠绕有限元模型建立

文中插件实现的功能是对已设计完成的壳体结构进行快速化建模功能，侧重体现在建模功能，线型设计部分并非此插件功能，故直接取网格理论设计好的壳体铺层方案，使用合肥工业大学研发的缠绕软件生成对应模型的缠绕层信息文件。使用本插件将信息文件导入到ABAQUS中生成壳体的有限元模型。根据有限元软件中仿真结果与对壳体的水压爆破实验进行对比，观察壳体在不同压力阶段的应变、位移和失效形式等。

2.1 复合材料壳体结构

本文使用的壳体芯模主要由橡胶内衬、金属接头和砂芯组成。芯模的具体尺寸如图6所示，直径为150 mm的壳体，左右极孔尺寸分别为30 mm和50 mm。

图6 芯模尺寸Fig.6 Mandrel dimension

通过网格理论对壳体铺层进行设计，设计爆压为30 MPa。设计考虑封头承压能力弱，适当增加缠绕层，减少环向层。由于本插件可以实现补强建模功能，故增加了补强层辅以验证。其中，缠绕补强是通过将筒身段的复合材料层剪断，只保留封头段复合材料层作为补强层，在有限元模型中表现形式如图7所示，设计的缠绕方案如表1所示。

图7 补强区域模型图Fig.7 Model diagram of the reinforce part

表1 铺层方案Table 1 Layup scheme

2.2 材料参数

采用干法缠绕进行实验，选用T700-12k纤维预浸带，中温固化为120 ℃，预浸带厚度0.15 mm，树脂含量32%，纤维含量150 g，带宽6 mm。

材料参数见表2。

表2 T700-12k预浸带力学性能参数Table 2 Mechanical properties of T700-12k prepreg

2.3 有限元模型导入

在缠绕软件中输入芯模尺寸参数和对应的铺层方案，生成缠绕层模型及导出的缠绕层封头角度如图8所示。

图8 复合材料壳体线型仿真图Fig.8 Composite case line type simulation

快速化建模功能是导入缠绕软件生成的数据文件进行数据处理生成模型。基于本课题组自研的缠绕软件导入生成的模型，自研缠绕软件中纤维缠绕层厚度分布采取的三次样条函数法进行预测，矫维成等对三次样条函数法预测厚度进行了理论推导和实验验证，证明了三次样条函数厚度预测的准确性，在纤维不滑纱的情况下，实际纤维最大厚度在一个带宽附近。

为节约计算时间，建立复合材料壳体的1/72模型。将缠绕层信息文件导入到插件中，壳体缠绕层采用C3D8R单元，生成的壳体模型已经完成了材料属性赋予和网格划分工作，完成了壳体的前处理工作。本建模流程一般只需几分钟即可完成，实现了壳体的快速化建模功能。

因缠绕软件中只能生成缠绕层信息文件，不包含金属接头部分。在ABAQUS中建立金属接头模型，将壳体与金属接头装配。对壳体中间部分施加轴向位移约束，避免有限元计算中发生刚体位移，壳体和接头模型侧边施加循环对称条件。在壳体内表面施加均布压力，从0～30 MPa依次增加。

边界及载荷施加如图9所示。

图9 载荷与边界条件Fig.9 Load and boundary conditions

2.4 壳体损伤模型

渐进损伤模型(Progressive Damage Model)已被广泛应用于复合材料结构的渐进破坏和失效分析中。首先，选取适合的失效准则表征材料失效状态，根据不同失效形式的退化模型，对材料进行刚度折减，当损伤到达一定程度后，整体结构失效，完成计算分析。渐进损伤计算流程如图10所示。

图10 渐进损伤流程图Fig.10 Flow chart of progressive damage process

针对壳体结构失效分析，主要考虑的失效模式有纤维拉伸压缩破坏、基体拉伸压缩破坏、纤维基体剪切破坏和分层破坏。

纤维拉伸破坏：

(3)

纤维压缩破坏：

(4)

基体拉伸破坏：

(5)

基体压缩破坏：

(6)

纤维基体剪切破坏：

(7)

分层破坏：

(8)

式中 下标11、22、33表示正轴坐标系下1、2、3方向的正应力；、、为剪切力；、分别为表单层板1方向的拉伸、压缩强度；、分别为单层板2方向的拉伸、压缩强度；、、分别为单层板各向的剪切强度；、分别为3方向的拉伸和压缩强度；α为经验系数，=1。

应力满足破坏条件时，会对材料损伤部分进行刚度折减。本文采用CAMANHO材料退化模型进行渐进失效分析。材料退化参数如表3所示。

表3 材料性能退化参数Table 3 Degradation parameters of material properties

3 结果与讨论

3.1 实验结果分析

为壳体设计水压爆破实验，如图11所示，监测爆破中壳体的应变位移变化。在壳体上均匀设置应变监测点，具体位置如图12所示。

图11 实验测量示意图Fig.11 Schematic diagram of experimental measurement

图12 应变测量点分布Fig.12 Distribution of strain measurement points

由于升压的过程中，壳体变形量较大，壳体伴随有树脂开裂，导致有一些应变片损坏，应变值无法显示，取数据完整的测试点1、3、5、7、9、13、19和21的结果，如图13(a)所示，其中每个应变测点监测两个方向的应变变化，取纤维方向为a，垂直于纤维方向为b，应变单位为με。在小端极孔处的13号测点应变产生了突变。位移变化监测设置两点，与应变监测点3、8重合。提取两点位移结果如图13(b)所示。两点均在27 MPa后产生了位移突变。结合实验中应变位移测试结果可知，此壳体的实际爆破压强为27 MPa。

(a)Tendency of strain with pressure (b)Tendency of displacement with pressure图13 应变与位移曲线Fig.13 Strain and displacement curves

壳体爆破实验结果如图14所示。可见，小端极孔接头的纤维从接头的肩部断裂，接头飞出。因为封头肩部具有较大的应力集中现象，导致纤维方向上的应力过大。由于封头肩部轴向承载能力大于环向承载能力，在均布载荷下，环向承受应力大于轴向，从而极易导致基体的损伤。随着环向应力的增大，基体损伤加剧，而基体损伤也会引起纤维的损伤。直到纤维损伤到达其强度极限后，纤维最终会断裂，从而导致封头肩部的纤维全部剪断，金属接头飞出，轴向位移突变。

图14 实验结果Fig.14 Experimental result

3.2 有限元结果分析

通过观察壳体有限元计算结果对比可知，壳体的主要失效原因为纤维失效，失效位置发生在壳体小端封头处，这里取壳体有限元计算中纤维拉伸失效结果进行分析，失效部位与未发生失效部位取高对比度的色阶便于观察。当内压为22.8 MPa时，整个复合材料缠绕层并未发生失效，如图15(a)所示。随着内压增加到23.2 MPa时，小端极孔封头肩部附近纤维发生失效，如图15(b)所示。随着内压逐渐增加，复合材料层逐层的纤维不断发生失效。当内压最终达到26.02 MPa时，金属接头肩部附近复合材料层纤维发生贯穿性失效，最终导致封头在该位置发生爆破，如图15(c)所示。

(a)p=22.8 MPa (b)p=23.2 MPa (c)p=26.02 MPa图15 纤维破坏云图Fig.15 Damage contours of fiber

除封头肩部位置，在筒身赤道处发生的损伤较多，但没有形成致命的失效。壳体的轴向位移云图如图16所示，提取小端封头段某节点的径向位移和轴向位移，如图17所示。在内压为26.02 MPa时，壳体小端极孔封头轴向位移发生突变，而径向位移并未发生突变，表明此时壳体发生轴向爆破，即外在表现形式为金属接头轴向飞出。

图16 轴向位移云图Fig.16 Axial displacement contours

图17 封头位移曲线Fig.17 Head displacement curves

通过实例分析对比可知，本建模方法下建立的有限元模型可以较准确地对复合材料壳体失效形式进行分析。只需从缠绕软件中生成线型文件，然后在插件界面中配置相关参数，即可快速生成壳体的有限元模型。如果计算结果不满足设计要求，只需更改铺层参数等相关信息，或者在缠绕软件中设计线型，重新计算，直到满足设计要求即可。在满足壳体强度刚度的要求下，也可对壳体结构进行优化设计，从而达到节约成本的目的。

4 结论

(1)针对固体火箭发动机复合材料壳体，提出了一种基于ABAQUS二次开发的快速化建模方法。通过对缠绕软件生成的缠绕层信息文件进行数据处理，导入到ABAQUS中，实现从缠绕设计软件到有限元分析软件的快速化建模方法。

(2)在数据处理过程中，针对波动性较大的数据，通过B样条函数对厚度进行光滑处理，保证生成模型的平滑性。

(3)通过壳体水压爆破实验结果的位移和失效形式分析，对比本方法生成的有限元仿真计算结果，对建模方法的准确性进行了验证。

(4)对上述建模方法设计了可视化界面，方便设计人员使用，大大缩短了设计的建模时间。在实际工程中，固体火箭发动机复合材料壳体设计与分析有着重要的意义。