基于有限元法的空间探测器复合材料结构安全裕度分析

姚 骏，孔祥宏，蒋 锋，张红英，王建炜，王志瑾

（1.南京航空航天大学 航空学院，江苏 南京 210016；2.上海卫星工程研究所，上海 201109；3.上海航天技术研究院，上海 201109）

0 引言

纤维增强复合材料强度、刚度比较高，广泛应用于航空、航天等领域。在工程应用中，复合材料结构的刚度、强度是设计和分析的重点。对于卫星等航天器结构，通常以刚度设计为主，兼顾结构强度、尺寸精度和稳定性。赵发刚等［1］对卫星复合材料结构的在轨检测进行了研究，为航天器复合材料结构的设计、制造及维护提供了技术参考。马立等［2］对复合材料结构的热变形进行了研究，成功研制了尺寸高稳定性的卫星桁架结构。复合材料结构的强度设计、分析常见于航空、风机等领域［3-5］。

由于复合材料具有可设计性，因此在航天器的刚度设计中，复合材料结构应用比例非常高。由于复合材料结构强度、剩余强度、分层破坏、损伤扩展、断裂失效等分析需要耗费大量的工作［6-8］，因此，在卫星复合材料结构的强度设计中采用首层破坏的强度准则，同时采用安全裕度（Margin-of-Safety，MS）对结构强度、安全可靠性进行评估。

空间探测器结构在研制过程中采用有限元分析方法进行强度校核，代替零件、组件、部件级的大量强度试验，以节省研制周期和成本。为了提高探测器复合材料结构强度分析的效率和准确性，本文拟提出一种基于有限元方法的复合材料结构安全裕度分析方法，通过Python 脚本实现安全裕度云图可视化，改变以往只能查看应力云图的情况，为复合材料结构强度、安全裕度的评估提供直观的云图。

1 复合材料安全裕度

结构安全裕度分析基于强度分析，结构强度分析通常基于应力分析（或应变分析）。在结构强度或安全裕度分析工作实践中涉及载荷、应力（应变）、安全系数、强度、含损伤结构剩余强度、安全裕度等要素。

1.1 各要素关系

1）载荷与安全系数

载荷通常指结构受到的力，通常包括使用（或限制载荷）、设计载荷（或极限载荷）。设计载荷与使用载荷的比值为安全系数。

2）强度、强度比与安全裕度

在强度分析工作中，需要计算载荷作用下的结构应力（或应变），再把所得应力（或应变）与材料强度进行强度计算。不同材料有不同的强度准则，例如材料力学中介绍常用的4 种强度准则、复合材料力学中介绍的各种强度准则和理论。对于强度分析，通常需要计算强度比R，如下所示：

式中：σ 为应力；［σ］为材料强度。若强度比小于1，则结构强度不足。

从安全性考虑，任何结构都要有一定的强度余量，在工程中用安全裕度来衡量，安全裕度MS计算如下：

对于不同材料进行结构设计时，需要根据标准、规范选取不同的安全裕度，并进行安全裕度指标考核。

3）载荷与安全裕度

对于塑性金属材料，其屈服、极限强度不同，极限强度与屈服强度的比值也不同。因此在进行安全裕度分析时，要充分重视使用载荷-屈服强度、设计载荷-极限强度条件下的安全裕度，并取两者中较小值作为结构的安全裕度，即利用载荷（限制载荷）作用下的结构应力与材料屈服强度进行安全裕度计算，得到安全裕度值MS1；利用设计载荷（极限载荷）作用下的结构应力与材料屈服强度进行安全裕度计算，得到安全裕度值MS2；取MS1和MS2中的较小值作为结构的安全裕度值。造成2 种载荷、2 个材料强度下安全裕度不同的原因是材料极限强度与屈服强度的比值与安全系数不同。

4）失效系数、强度比与安全裕度

失效系数If（Failure Index）与强度比R、安全裕度MS的关系为

对于层压复合材料，需要对每层进行安全裕度分析，甚至对层间受力情况进行安全裕度分析。复合材料的失效准则通常有最大应力、最大应变、帕克（Puck）、哈辛（Hashin）、蔡-希尔（Tsai-Hill）准则及蔡-吴（Tsai-Wu）张量理论。而使用Hashin、Tsai-Hill、Tsai-Wu 等失效准则计算的失效系数与应力为二次方关系，以Tsai-Hill 失效准则为例，失效系数计算为

式中：σ11、σ22、σ12分别为单层板纤维方向应力、基体横向应力、面内剪应力；X、Y、S分别为单层板纤维方向强度、基体横向强度和面内剪切强度。

在使用式（4）失效系数计算安全裕度时，需对失效系数开根号，即：

5）失效类型及损伤扩展

复合材料结构的失效类型比较多，主要分为层内失效和层间失效。层内失效包括多种形式，例如：纤维拉/压/屈曲、基体拉/压、纤维与基体层内剪切、法向挤压失效。层间失效主要为分层，包括张开型（拉开）、撕开型（剪开）。面内受压状态下分层后可能进一步导致分层处局部屈曲。根据失效系数可以判断是否失效，但失效类型的确定需要结合应力组合情况进行判断。在损伤演化分析中，需要根据失效类型对失效区域对应的刚度按照一定幅值进行折减，并进一步进行应力、失效分析，从而完成损伤演化或损伤扩展分析。损伤区域刚度折减的幅度通常用损伤变量控制。

6）损伤扩展与断裂失效

复合材料结构的断裂情况与结构形式、受力情况有极大的关系，出现的断裂严重程度也不尽相同。例如，可能发生局部全厚度、局部部分铺成、整体断裂。以复合材料杆件为例，破坏形式可能有接头处部分铺成断裂。由于胶接工艺限制，接头与杆件表层胶接，在承载情况下接头拉脱造成杆件部分铺成断裂。如果接头处采用较强连接方式，则可能造成杆件某处横截面全断裂。从理论及数值分析角度看，复合材料结构损伤有产生、发展的过程，即损伤扩展，当损伤扩展到一定程度时引起结构断裂失效。

损伤扩展并不是引起断裂失效的唯一原因。因为在复合材料结构中，复杂受力状态下，不仅能够造成材料失效从而造成结构失效，而且也极可能造成结构失稳，瞬间破坏。因此，需根据具体情况辩证看待损伤扩展与断裂失效的关系。

1.2 强度分析方法

在复合材料层压板的强度分析中，通常要对每一个铺层进行强度分析。强度分析方法通常分为2 类：1）不考虑损伤扩展、渐进失效的应力分析，只需一次计算，根据得到的各铺层应力和相应强度准则计算结构强度；2）考虑损伤扩展、断裂失效的应力、破坏载荷分析，需要多次迭代计算，在迭代过程中按照一定的强度准则和刚度折减方法对材料刚度进行折减。

对于第2 种方法，孔祥宏等［9］使用ABAQUS 有限元软件和Python 前处理程序、USDFLD 用户子程序（Fortran 子程序）对复合材料结构的渐进损伤强度分析进行了研究。郑双等［10-13］使用ABAQUS和VUSDFLD 子程序对复合材料结构进行了渐进损伤强度分析。此外，还可以通过ABAQUS 和UMAT 子程序进行类似分析。但是，这种分析方法通常关注结构在破坏过程中每个铺层的应力、材料刚度的变化，很难从结构整体评估结构的强度。

而对于第一种方法，经过计算可以得到复合材料层压板每个铺层的应力，将应力值代入相应的强度准则（例如Hashin、Tsai-Hill、Tsai-Wu 等强度理论），可以得到失效系数If，进而可以得到每个铺层的强度比R及安全裕度。

1.3 复合材料安全裕度

对于复合材料结构，通常采用首层失效原则计算结构安全裕度。因此，要对每个铺层进行安全裕度分析。以层压板为例，得到各铺层的安全裕度后，取其中最小值作为层压板的安全裕度：

式中：MSLam为层压板的安全裕度；MSi为第i层的安全裕度；n为层数。

由此计算安全裕度为一个值，对于复合材料结构，用一个值反映其安全裕度从指标符合性上可以达到目的，但从结构性能评估、改进角度来看，一个安全裕度数值的意义不大。为了识别复合材料结构不同区域强弱的情况，可以对结构每一个细分区域进行安全裕度计算，从而了解整个结构的安全裕度情况。使用安全裕度云图可以非常直观地评估结构力学性能。

1.4 工程应用方法

实际工程结构中的复合材料结构通常比较复杂，经常存在局部应力集中现象。应用式（6）安全裕度对复合材料结构进行安全性评估时通常过于保守。因此，在工程上，可以计算多个位置的安全裕度，也可用有限元的方法计算每个单元的安全裕度，并以云图的形式显示安全裕度分布情况。

以复合材料层压板为例，在层压板的有限元模型中对每个单元赋予n层复合材料铺层的材料参数，即相当于每个单元有n层，单元在铺层法向上的积分点为n的整数倍。用某一层积分点的应力状态表示相对应的某一铺层应力状态。由此计算一个单元安全裕度的方法为

式中：MSE为一个单元的安全裕度；为单元第i层的安全裕度，使用式（2）和式（3）求解；n为赋给单元的复合材料铺层的层数。

将有限元模型中所有单元的安全裕度MSE的值以云图形式显示，可以观察整个模型的安全裕度分布，为安全裕度低的结构补强和安全裕度高的结构减重提供依据。

2 安全裕度可视化方法

复合材料结构静力分析得到应力结果为安全裕度计算的输入。使用ABAQUS 有限元软件和Python 前后处理程序可以将单元的应力经计算后得到安全裕度，并将安全裕度值通过有限元模型进行可视化处理。可视化处理的方法有2 种，分别为分组可视化、云图可视化。

2.1 安全裕度计算

复合材料结构有限元模型中单元安全裕度计算流程如图1 所示。图中，Ei为有限元模型中第i个单元；MSi为第i个单元的安全裕度；MSj为第i个单元第j个铺层的安全裕度；nPly为第i个单元的铺层数；nElem为有限元模型的单元数。单元的安全裕度取各铺层安全裕度的最小值，在计算过程中需要比较MSi与MSj的大小，在计算每个单元的安全裕度前，先给MSi设置一个合理的较大初始值，方便后续计算。

图1 复合材料结构单元安全裕度计算流程Fig.1 Calculation process of MS for the elements of the composite structure

上述计算过程可以通过编写Python 程序从ABAQUS 分析结果文件（ODB 文件）中读取单元应力数据，并进行计算。图1 中仅表示复合材料结构单元安全裕度计算的大致流程。具体编写Python后处理程序时，需要根据ABAQUS 分析结果文件的数据库结构编写相应Python 程序。

2.2 分组可视化

分组可视化是比较简单的可视化操作。根据计算得到的单元安全裕度值，将单元进行分组，为每一组单元定义一个颜色，通过在有限元模型上显示各组单元的颜色来观察安全裕度的分布。该方法首先将单元编号分组，根据单元编号，可以在有限元模型（CAE 文件）中创建单元集合，也可以在结果文件（ODB 文件）中创建单元集合。然后，按照单元集合设置显示颜色。对此，需要编写不同的Python 前处理或后处理程序。

2.3 云图可视化

云图可视化方法是将单元安全裕度的值作为场变量写入ABAQUS 的结果数据库文件（ODB 文件）中，以云图的形式显示模型的安全裕度。由于每个单元只有一个安全裕度值，因此不需要显示各铺层的安全裕度，而是直接显示各铺层中最小的安全裕度值。通过编写Python 后处理程序，可以在ODB 文件中快速创建场变量，从而可以进行可视化操作。

3 安全裕度分析

火星环绕器在结构设计过程中进行了多次模态、静力、强度和安全裕度分析。通过多轮迭代分析，改进设计，局部补强，使火星环绕器的结构满足刚度和强度设计要求。根据火星环绕器的静力、强度分析结果，使用云图可视化的方法对火星环绕器复合材料结构系统进行安全裕度分析。根据图1 安全裕度计算流程，结合ABAQUS 结果文件数据库结构，针对多部件有限元模型在多个工况下的分析结果，编写可以对多部件、多工况结果进行数据后处理的Python 程序。

3.1 结构改进前分析结果

火星环绕器结构分系统在某工况下静力分析结果的位移云图，如图2 所示。复合材料层压板的应力每次只能显示一个铺层的一个应力分量云图，此处不再给出应力云图。火星环绕器结构最大位移为3.99 mm。

图2 火星环绕器结构改进前位移云图Fig.2 Displacement contours of China’s Mars orbiter before structure modification

火星环绕器主承力结构安全裕度分析结果如图3 所示，图例中MS 为安全裕度，使用Python 程序在ODB 文件中创建场变量时定义。火星环绕器内部十字形结构安全裕度较低，有大面积区域安全裕度低于0，最低为-0.966，严重不符合设计要求，需要对结构进行改进。

图3 火星环绕器结构改进前安全裕度云图Fig.3 MS contours of China’s Mars orbiter before structure modification

3.2 结构改进后分析结果

由于受其他结构影响，火星环绕器内部十字板的某个最佳传力路径被打断，因此提高其安全裕度的方法为进行局部补强。结构改进后火星环绕器结构分系统在某工况下静力分析的位移云图，如图4 所示。结构改进后，火星环绕器结构最大位移为2.48 mm。

图4 火星环绕器结构改进后位移云图Fig.4 Displacement contours of China’s Mars orbiter after structure modification

结构改进后，火星环绕器主承力结构安全裕度分析结果如图5、图6 所示。由图5 可知，主承力结构大部分区域安全裕度大于0，仅有很小的边角区域安全裕度小于0。在图6 中，大面积的灰色区域安全裕度大于1，剩余彩色部分安全裕度小于1，只有十字板上边拐角区域安全裕度小于0。

图5 火星环绕器结构改进后安全裕度云图Fig.5 MS contours of China’s Mars orbiter after structure modification

图6 火星环绕器结构局部安全裕度云图Fig.6 MS contours of partial China’s Mars orbiter

由图2 和图4 对比可知，结构改进后，火星环绕器结构分系统的刚度提高。由图3 和图5、图6 对比可知，结构改进后，火星环绕器主承力结构复合材料部件的安全裕度大幅提高。

4 结束语

复合材料结构安全裕度分析可以从整体了解结构的强度、安全可靠性。复合材料结构的强度分析通常可以得到各铺层的强度情况，将各铺层中强度最弱的区域投影到一个面上，可以得到结构较弱区域的最大范围。本文的安全裕度分析利用该思想实现了对复合材料结构整体安全裕度的计算，以云图可视化的形式显示了结构安全裕度，便于观察和评估结构性能。

1）通过对火星环绕器结构分系统进行静力分析和安全裕度分析，得到主承力结构的安全裕度云图。

2）通过与静力试验结果（试验图片略）对比，发现试验中结构发生问题的区域均在图3 中安全裕度小于0 的区域内。

3）通过安全裕度分析，可以对火星环绕器复合材料结构的改进和补强提供依据。根据分析结果和试验结果对火星环绕器主承力结构进行局部补强后，结构满足强度要求。结构补强后的安全裕度分析结果也验证了该结论。

4）在ABAQUS 中利用Python 程序进行前后处理工作效率高，可以将任何已有的包含应力数据的分析结果进行安全裕度分析。避免编写Fortran 用户子程序（如UMAT、USDFLD 等）后，再对结构进行分析，从而节省50%以上的工作量。