尤军峰，刘 浩，王春光

(1.中国航天科技集团公司四院四十一所 固体火箭发动机燃烧、热结构与内流场国防科技重点实验室，西安 710025；2.西安近代化学研究所，西安 710065)

0 引言

固体火箭发动机壳体主要用来承受药柱燃烧产生的内压载荷[1]，对于各向同性材料的金属壳体，其承压能力及其结构分析，已有成熟的理论及较为简单的有限元方法可利用，而对于复合材料壳体，用有限元法来分析，也要困难和复杂得多，其原因主要是除了需要用到各向异性层合薄壳理论，还要考虑固体火箭发动机壳体封头上每个缠绕层纤维的缠绕方向(缠绕角)和层厚从赤道到极孔是不断变化的。目前，工程上常用的简化方法是将封头沿母线方向分成若干段，计算每一段的平均半径，将每一段认为正交各向异性材料，用层合板理论等效出9个弹性常数用宏观的方法进行处理，忽略封头处缠绕层厚度及缠绕角度不断变化，这种方法只能粗略地计算变形，而对于需要获得每个缠绕层纤维的应力、应变以及层间的剪切应力就无能为力了，更无法模拟封头部位每个纵向缠绕层间经纬双向织物的补强，无法对复合材料壳体封头部位进行精确计算和受力分析，使得封头部位的型面设计和补强结构设计更多依赖工程经验，无法实现精细化设计。

本文通过编制计算程序，精确计算出封头部位每个纵向缠绕层的缠绕角和层厚，并考虑封头补强，采用ANSYS ACP模块铺层建模，较真实地模拟了固体火箭发动机壳体的实际结构。

1 理论及方法

1.1 复合材料有限元数值计算方法

复合材料是由两种或两种以上不同性质的材料，通过物理或化学的方法，在宏观上组成具有新性能的材料。各种材料在性能上互相取长补短，产生协同效应，使复合材料的综合性能优于原组成材料而满足各种不同的要求。复材可通过加强材料的类型进行分类：颗粒增强型复合材料和纤维增强型复合材料，纤维增强型复合材料又分为短纤维增强复合材料和长纤维增强复合材料。固体火箭发动机复合材料壳体属于长纤维增强型复合材料，它是将炭纤维(或玻璃纤维)预浸胶后通过缠绕设备逐层缠绕到芯模上，在一定温度条件下固化而成。目前，复合材料数值仿真方法分为三种：

(1)微观方法(Fiber Level)。微观方法考虑了纤维和基体的微观结构、形状、位置以及材料属性,是最详细、真实地体现了复合材料的细微结构，但由于纤维的直径往往只有复合结构尺寸的百万分之一，采用微观方法建立有限元模型，将会导致节点与单元数量非常多，无法计算。

(2)宏观方法(Laminate Level)。宏观方法是将复合材料假设为正交各向异性材料(即假定沿铺层方向等效为各向同性材料)，采用静力等效的方法获得正交各向异性材料的9个性能参数：Ex，Ey，Ez，νxy，νyz，νxz，Gxy，Gyz，Gxz。如果只是需要计算结构位移、屈曲载荷、振动频率以及模态，则可以把多层的复合材料结构处理成均匀的正交各向异性材料单层板，这时无法得到板内部各层及层间的应力分布。

(3)细观方法(Ply Level)。细观方法是用单层(single layers)材料来建立层合板(壳)的方法来分析复合材料结构的应变、应力和失效。与微观方法相比，极大减少了节点和单元数量。与宏观方法相比，增加了复合材料每层及层间的应变、应力和失效分析的分析功能。因此，细观方法(Ply Level)是比较先进的数值计算方法，ANSYS Composite Prepost(ACP)正是采用细观方法来建立复合材料铺层模型、分析复合材料结构。细观方法需输入每个单层材料属性、厚度和纤维铺设角等，最大限度地实现复合材料的精细化几何模型建模。

1.2 单层材料参数计算

单层材料(纤维束或复丝)结构如图1所示。其中，1方向平行于纤维方向，2方向垂直于纤维方向，3方向垂直1、2方向的平面[2]。

图1 单层材料示意图

根据单层材料中纤维与基体的不同配比，可得：

(1)纤维体积分数

(2)纤维质量分数

(3)基体质量分数

(4)纤维体积分数与基体质量分数之间的转换

单层材料可认为是横观各向同性材料，独立的弹性常数有5个，即E1、E2、G12、G23、ν12，在ACP中单层材料的输入按照正交各向异性对待，输入材料参数9个，即E1、E2、E3、G12、G13、G23、ν12、ν13、ν23，分别代表单层材料三个方向的弹性模量、剪切模量及泊松比，这9个参数的计算由混合法则决定，混合法则用于根据纤维和基体的性能估算单层材料的性能。文献提供了多种组合法则，最常见的是Jones法则[3]：

E1=φ·EF1+(1-φ)·EM；

E3=E2；

G13=G12；

G23=φ·GF23+(1-φ)·GM；

ν12=φ·νF12+(1-φ)·νM；

ν13=ν12；

ν23=φ·νF23+(1-φ)·νM

1.3 ANSYS ACP模块

ANSYS Composite Prepost(ACP)是集成于ANSYS Workbench环境的全新的复合材料前/后处理模块，可与ANSYS其他块实现数据的无缝连接。ANSYS ACP模块的理论基础是层合薄壳理论，在处理层合复合材料结构方面具有可视化、易操作的优点，ACP具有强大的结果后处理，可获得各种分析结果，如层间应力、应力、应变、最危险的失效区域等[4]。ANSYS ACP操作界面见图2。

对于简单的铺层结构，可直接在“Rosettes”中定义纤维的参考方向(并可参考直角坐标系、柱坐标系、球坐标系等)，但对于固体火箭发动机复合材料壳体来说，在封头上每个缠绕层纤维的缠绕方向(缠绕角)和层厚从赤道到极孔是不断变化，必须根据封头型面(方程)计算出封头纵向缠绕层缠绕角和单层厚度[5]。

1.4 封头纵向缠绕层缠绕角、层厚计算

一般情况下，复合材料壳体封头为椭球型面，如图3 所示。其中，r0为极孔半径，R0为赤道半径，x、y、z为全局坐标系，X、Y、Z为局部坐标系，L和L′为局部坐标系位置的纵向纤维缠绕方向，θ和-θ为局部坐标系位置的纵向纤维缠绕角，见图3[2]。

图2 ANSYS ACP操作界面

图3 椭球模型示意图

假设封头为椭球型面，方程为

当a=b时为旋转椭球面

则在封头上任意点(x0，y0，z0)处法向量Z为

又由(0，0，0)、(0，0，1)及(x0，y0，z0)三点决定的平面方程为

-y0(1-z0)(x-x0)+z0x0(y-y0)+x0y0(z-z0)-

y0z0(x-x0)+x0(1-z0)(y-y0)-x0y0(z-z0)=0

整理后，得

-y0(x-x0)+x0(y-y0)+0=0

对比点法式平面方程：

A(x-x0)+B(y-y0)+C(z-z0)=0

则(x0，y0，z0)点向量X为(-y0，x0，0)

则(x0，y0，z0)点向量Y为

由于X·Y=0，Y·Z=0，Z·X=0，因此可由X，Y，Z构成新的直角坐标系，坐标原点O(x0，y0，z0)，在新坐标系XYZ下，Y向量为(0，1，0)。

设向量Y向右转动θ角度(右手法则)后，形成的新向量为L，向量Y向左转动θ角度(-θ)后，形成的新向量为L′，则在新坐标系XYZ下，L为(-sinθ，cosθ，0)，L′为(sinθ，cosθ，0)。

根据坐标变换公式，按表1给出新坐标轴OX、OY、OZ的方向余弦时：

表1 方向余弦

则有

即有表2和表3结果。

表2 方向向量

表3 单位方向余弦

则L为

L′为

根据上述封头上各处纤维的参考方向以及铺层的厚度计算方法，可编制成计算机程序计算，并将计算结果import用到 Look-up Table中，供ACP建模参考引用。

2 算例

2.1 计算模型

以某型号发动机复合材料壳体为例(图4)，开展内压载荷下的壳体受力状态分析。

如图4所示，壳体半径1000 mm，前开口比为0.201，后开口比为0.433，环向层与纵向层厚度之比为1.085，后封头补强层数为10层，筒段平均缠绕角为18.5°，设计内压载荷为11.5 MPa等。由于壳体具有轴对称性，计算模型简化为1/60，如图5所示。

图4 某发动机壳体结构简图

图5 计算模型

2.2 材料参数

材料模型建立所需的复合材料单向板、补强炭布、树脂、金属及橡胶的材料参数见表4，主要包括弹性模量、泊松比和剪切模量3个参数[6-7]。

表4 复合材料壳体材料参数

2.3 封头缠绕角及厚度计算

(1)缠绕角(矢量)

在ACP模块根据封头母线上各离散点的坐标值计算出各离散点处的缠绕角，并创建出Excel 格式的Lookup Table，如图6所示，用于对封头各点位置处单层缠绕层材料弹性常数的计算及赋值。

(2)厚度

在ACP模块根据封头母线上各离散点的坐标值计算出各离散点处的厚度，并创建出Excel 格式的Lookup Table，如图7所示，这样可根据各缠绕层的缠绕角，计算出层合复合材料的弹性常数。

图6 缠绕角计算结果

图7 缠绕层厚度计算结果

2.4 模型建立

在Excel 格式的Lookup Table创建后，采用ACP建模可准确模拟壳体封头部位纤维缠绕方向、缠绕层退移或切根、碳布补强，在赤道及筒段部位可模拟裙内纵向缠绕及环向缠绕、裙尖填平、裙外环向缠绕及补强等，如图8所示。复合材料结构有限元模型建立完成后，可和各向性材料的其他件的有限元模型合并，如图9所示。最终有限元模型如图10所示。

(a)封头部位纤维缠绕方向

(b)裙内缠绕层纤维方向

(c)裙外缠绕层纤维方向

(a)与裙及橡胶层有限元模型装配

(b)与后接头等有限元模型的装配

图10 有限元模型

2.5 载荷及约束

载荷及约束如图11所示，施加对称边界条件A、B，约束轴向位移C，内表面施加载荷D：11.5 MPa。

图11 载荷及约束

2.6 计算结果及分析

有限元计算得到壳体位移如图12 所示。可见，壳体加大变形发生在封头中部及壳体圆筒段，最大变形量10.6 mm。在变形较大的封头部位沿纤维方向最大应力值为921.85 MPa，见图13。应变值为7065 με，见图14，最大主应力方向沿纤维方向，如图15所示。 在变形较大的圆筒部位沿纤维方向最大应力值为1559.6 MPa，应变值为10 433 με，如图16和图17所示。

图12 壳体变形

图13 封头部位最大应力

图14 封头部位应变

图15 封头部位应力主方向

图16 筒段最大应力

试验测得筒段位移值为 10.5 mm ，封头部位沿纤维方向最大应变值为7171 με，圆筒部位沿纤维方向最大应变值为9599 με。计算结果与水压试验结果表明，在水压11.5 MPa水压条件下，最大应力1559.6 MPa小于纤维方向材料强度3500 MPa，强度满足要求；封头部位应变计算精度为1.5%，圆筒部位应变计算精度为7.9%。分析原因：计算中用到的缠绕角(18.5°)为平均缠绕角，而在赤道部位实际缠绕角约26°，即从筒段中部到后封头赤道缠绕角从18.5°过渡到26°，导致靠近后封头的圆筒段实际环向刚度增大。因此，测得的环向应变值偏小。

3 结论

本文编制计算程序，计算出封头部位每个纵向缠绕层的缠绕角和单层厚度，采用ANSYS软件中的 ACP模块，充分考虑各补强层的铺放情况，建立了固体火箭发动机复合材料壳体精细化几何分析模型，发挥计算编程与商用软件相结合的优点，提供了一套切实可行的固体火箭发动机复合材料壳体结构细观力学分析方法。计算算例表明：

(1)该方法集成于ANSYS Workbench，可直观地进行复合材料铺层精细化设计与几何建模，具有交互性强、易操作、效率高、计算精度高等优点；

(2)仿真计算结果与试验结果误差可控制在8%以内，计算精度满足工程要求。