面向缺陷容忍的加筋筒壳快速优化设计

马祥涛，王法垚，朱英杰，郝鹏，*，王博，刘观日

1．大连理工大学 工程力学系 工业装备结构分析国家重点实验室，大连 116024

2．北京宇航系统工程研究所，北京 100076

加筋筒壳由于具有较高的比刚度和比强度，被广泛应用于运载火箭等航天装备的主承力结构当中［1］。其在发射过程中承受过载导致的巨大轴压载荷，主要破坏模式为屈曲失稳。因此，轴压载荷下的极限承载能力是航天筒壳结构设计过程中的关键考核指标［2］。然而，筒壳结构的屈曲载荷表现为强缺陷敏感性，即使微幅的几何缺陷也会引起其极限承载的大幅折减，而缺陷又是加工、运输、装配等过程无法避免的。因此在航天薄壁筒壳结构的设计过程中必须考虑缺陷对其极限承载的影响。

在传统的筒壳结构设计过程中，是将完美结构的临界失稳载荷乘以一个较为保守的固定折减因子（Knock Down Factor， KDF）作为设计载荷。即其本质仍是以完美结构承载能力作为设计载荷，忽略了筒壳结构设计的改变对缺陷敏感性的影响，用较低折减因子增加安全裕度来保障结构的安全可靠，这无疑会导致结构超重。因此，为了实现航天薄壁筒壳结构的轻量化设计，各国航天中心和学者们纷纷开展了筒壳缺陷敏感性和折减因子预测相关的研究工作，其通过高精度预测筒壳结构的折减因子，减少不必要的安全裕度实现结构的轻量化设计。例如，美国NASA工程安全中心于2007年立项研究面向运载火箭中金属及复合材料薄壁结构的折减因子（Shell Buckling Knockdown Factor，SBKF）设计规范［3］，并提出了筒壳结构的实测几何缺陷法［4］；欧盟于2012年立项研究面向航天复合材料薄壳结构的缺陷敏感性方法（DESICOS）［5］，提出了单点扰动载荷法（Single Perturbation Load Ap⁃proach， SPLA）［6-7］；中 国 也 于2013年 立 项 研 究计及缺陷敏感性的网格加筋筒壳结构轻量化设计理论与方法［8］，王博等［9-10］针对加筋筒壳结构开展了基于非线性显式后屈曲分析的高精度数值模拟、缺陷敏感性分析、高精度折减因子预测等研究工作，提出了多点最不利扰动载荷法（Worst Multiple Perturbation Loads Approach，WMPLA），并基于筒壳轴压屈曲实验进行了方法验证［11］。上述筒壳折减因子方法已经得到各国航天总体部门的认可，并已经应用到实际型号当 中。此 外，Zhang等［12］还 提 出 了 一 种 基 于SPLA的概率随机扰动载荷法，用来预测筒壳的折减因子下限。Kim等［13］给出了一条基于SPLA的筒壳折减因子预测公式。还有很多学者开展了筒壳轴压屈曲实验，以为薄壁筒壳结构的设计提供帮助和指导。例如，Li等［14］开展了折叠点阵夹层筒壳结构的轴压屈曲实验。Jiao等［15］通过实验研究了非均匀轴压筒壳折减因子和载荷分布的关系和规律。

尽管通过折减因子的高精度预测，可以实现航天筒壳结构的轻量化设计，但其只能被动地挖掘设计空间，设计潜力不足。特别是在中国新一代大直径运载火箭对运载能力跨越式提高的要求下，箭体结构的直径和质量将大幅提高，结构轻量化需求导致航天复合薄壁筒壳等效厚度相对变薄，其缺陷敏感性问题将愈发突出，航天筒壳结构地轻量化设计变得更加困难。

因此，为了提高航天筒壳结构的承载效率，除了通过提高折减因子的预测精度，被动地挖掘减重设计空间，还需要进一步主动地挖掘设计空间。例如，Zhou等［16］基于拓扑优化方法开展了薄壁筒壳创新网格构型设计的研究，通过提出新颖的加筋构型，提高筒壳的屈曲载荷和抗缺陷能力。但为了进一步提高航天筒壳结构地轻量化设计水平，还需要开展面向缺陷容忍的筒壳结构设计，在设计过程中充分考虑结构设计与临界失稳载荷和结构缺陷敏感性的耦合关系，通过同步提升结构的屈曲载荷和抗缺陷能力，在保证航天结构安全可靠的前提下，实现筒壳结构的精细化和轻量化设计。其中郝鹏等［17-18］已经开展了面向缺陷容忍的加筋筒壳和加筋锥壳结构优化设计研究，并提出了一种能够显著提升筒壳结构抗缺陷能力的多级加筋增强构型［19］，可在给定的极限轴压载荷约束下，实现筒壳结构的大幅减重。但对于含加筋等丰富结构细节的复杂薄壁筒壳结构来说，考虑缺陷敏感性的单次承载分析极其耗时，其优化所需的计算成本更是难以接受，因此亟需开展面向缺陷容忍的加筋筒壳结构快速优化设计方法研究。

针对上述问题，基于不完全折减刚度法（in⁃complete Reduced Stiffness Method， iRSM），提出了一种面向缺陷容忍的加筋筒壳快速优化设计框架，以解决大直径加筋筒壳优化过程由于非线性显式后屈曲分析耗时导致计算成本高的问题。首先介绍了含缺陷非完善筒壳结构的快速承载分析方法，随后介绍了面向缺陷容忍的加筋筒壳快速优化设计框架，并以一个直径1.6 m的正置正交网格加筋筒壳结构为例，对所提出的优化框架进行详细说明。

1 含缺陷非完善筒壳结构的极限承载分析

1.1 面向非完善筒壳的非线性显式后屈曲分析

非线性显式后屈曲分析可以准确地模拟网格加筋、蒙皮桁条等航天薄壁筒壳结构的后屈曲行为，并可稳定地计算其极限承载能力，已被应用于网格加筋等航天薄壁筒壳结构的后屈曲分析和优化当中［18，20-21］。其可获得每一时刻的结构位移和当前的结构构型状态，并且不存在收敛性问题。

而对于含缺陷非完善的筒壳结构，可通过对模型初始构型的修正，将缺陷引入到数值模型后，再开展非线性显式后屈曲分析。对于有限元模型，可通过修调有限元模型节点坐标的方式，实现上述目的。其中，考虑初始几何缺陷的非完善加筋筒壳结构有限元模型初始构型可表示如下形式：

式中：X为含缺陷非完善筒壳结构有限元模型的初始构型状态；Xp为完善筒壳结构有限元模型的初始构型状态；XI为缺陷引起的有限元节点坐标偏移；δ为缺陷的最大幅值；α为无量纲缺陷幅值；hr为筋条高度；ts为蒙皮厚度；N为归一化后的缺陷基矢量，表示几何缺陷的形式，可以为模态缺陷、正弦波、凹坑缺陷等形式。

尽管非线性显式后屈曲分析可以较为准确地模拟加筋筒壳结构从初始状态到后屈曲直至破坏的全过程，并获得其极限承载能力。但由于筋条等结构细节特征会使得时间迭代步长很小，致使加筋筒壳结构的极限承载分析非常耗时，最终导致面向加筋筒壳结构的优化设计需要极高的计算成本。因此，亟须开展面向含缺陷非完善筒壳结构的承载快速分析方法研究，解决加筋筒壳等航天工程薄壳结构分析及优化效率低的瓶颈问题。

1.2 非完善筒壳的承载快速分析

为了实现含缺陷非完善筒壳结构的承载快速分析，笔者［22］基于折减刚度法（Reduced Stiffness Method， RSM），提出了一种用于替代非线性显式后屈曲分析方法的非完善筒壳承载快速分析方法——不完全折减刚度法（incomplete Re⁃duced Stiffness Method， iRSM），其通过建立筒壳缺陷和膜刚度的折减关系函数，可在线性屈曲的分析框架下实现含缺陷非完善筒壳的承载力快速计算。经算例验证，该方法可在满足计算精度要求的同时，实现计算成本的大幅降低。下面对折减刚度法和不完全折减刚度法分别做简要介绍。

折减刚度法是Croll［23］提出的一种用于快速估计筒壳等薄壁结构折减因子下限的方法。在折减刚度法中，薄壳结构的刚度被分为膜刚度和弯曲刚度2部分，其假设薄壳的膜刚度会由于初始缺陷的存在而折减甚至消失。具体而言，如果筒壳存在初始缺陷，在壳体初始后屈曲状态，会因缺陷致使的多种模态相互竞争作用导致膜刚度成分的腐蚀甚至消失，从而引发承载力的大幅折减。基于上述假设的折减刚度法可在线性屈曲的分析框架下，为筒壳折减因子下限的估计提供一个快速简单的计算工具，其相关公式如下所示：

其中，式（4）为完美模型的全刚度线性屈曲分析控制方程。式中：K为完美模型的刚度矩阵；KG为几何刚度矩阵；λ为完美模型的屈曲载荷；Φ为屈曲模态。式（6）为折减刚度模型的线性屈曲方程。式中：Km为结构的膜刚度；Kb为弯曲刚度；λ′为折减刚度模型的线性屈曲载荷；Φ′为对应的屈曲模态；β为膜刚度成分的衰减因子（β≥1）；η为所预测的折减因子。当β趋近于无穷时，筒壳的膜刚度成分即被完全消除，并可获得筒壳结构的折减因子下限。

折减刚度法能稳定地预测筒壳结构的折减因子下限是因为其对筒壳模型全域的膜刚度做了相同程度甚至完全的折减。但实际结构的初始缺陷具有随机性和不确定性，不同区域的缺陷表现各不相同。以几何缺陷为例，不同区域的几何缺陷幅值是完全不同和随机的，有的区域甚至趋近于完美模型，没有缺陷。根据假设，膜刚度的折减是由于初始缺陷诱发的。所以，应对筒壳具有不同缺陷水平区域的膜刚度分别进行不同程度的折减，并且其膜刚度折减程度应与缺陷水平相关。此外，通过把膜刚度成分完全消除来预测折减因子下限也是不合理的，若筒壳膜刚度成分为0，则其不能承受任何大小的轴压载荷，所预测的折减因子应该为0。因此，笔者［22］在折减刚度法的基础上提出了不完全折减刚度法，并且继承了折减刚度法的全部假设。

在不完全折减刚度法中，通过建立缺陷水平和膜刚度成分的折减关系函数f，可依据缺陷具体的水平和分布，对筒壳有限元模型各单元刚度进行准确的量化折减。然后通过组建并求解如式（8）所示的不完全刚度折减模型的线性屈曲方程，即可进行快速的缺陷敏感性分析，并获得该非完善筒壳结构的折减因子。

式中：Kem为单元的膜刚度矩阵；i为单元的编号；λ′′和Φ′′分别为不完全折减刚度模型的屈曲载荷和屈曲模态。另外，针对最为常见的几何缺陷，笔者［22］给出了一个几何缺陷相关的刚度折减关系函数范例，如下所示：

图1 不同参数下的刚度缺陷折减关系函数曲线Fig. 1 Curves of reduction function with different parameters

可以看出，参数A控制着刚度折减关系函数的收敛速度，参数B控制着膜刚度折减程度的极限。当单元几何缺陷水平为零时，不对单元的膜刚度进行折减，膜刚度折减程度随着几何缺陷水平的增加单调递增，并逐渐收敛到一个固定值1/(1+B)。因此即使几何缺陷水平非常大，壳体的膜刚度也不会完全消除。折减关系函数的单调递减趋势表明，筒壳缺陷水平越大，膜刚度折减越多，非完善筒壳折减因子也越小。另外，单元缺陷水平基于单元厚度做了归一化处理，表明相同几何缺陷幅值下，壳体越薄，对几何缺陷越敏感。刚度折减关系函数的上述特征都符合现有实验观测到的薄壁筒壳结构缺陷敏感性规律。因此，所提出的折减关系函数能够描述模型膜刚度和缺陷之间的关系，只需要确定折减关系函数中的待定参数，就能在线性屈曲分析框架下实现非完善筒壳的承载快速分析。

另外，笔者［22］已经针对含实测几何缺陷的筒壳结构，基于不完全折减刚度法开展了非完善筒壳的承载快速分析。结果显示，相比非线性显式后屈曲分析，不完全折减刚度法可在满足分析精度的同时，实现计算成本降低90%左右。但需要说明的是，目前不完全折减刚度法所采用的刚度折减关系函数具有一定的局限性，因为其仅对单一类型的几何缺陷有效，无法同时兼顾多种形式的几何缺陷并保持满意的分析精度。仍需继续开展相关研究，提出能够覆盖多种形式几何缺陷的刚度折减关系函数，并进一步揭示薄壁结构几何缺陷敏感性的内在机理。

2 面向缺陷容忍的加筋筒壳快速优化设计框架

为了解决计及缺陷敏感性加筋筒壳分析优化计算成本高的问题，本文基于不完全折减刚度法，提出了一种面向缺陷容忍的加筋筒壳快速优化设计框架。主要包含以下步骤：首先，基于加筋筒壳的初始设计，根据给定的几何缺陷形式确定刚度折减关系函数，建立含缺陷非完善筒壳结构的极限承载低保真度分析模型；然后使用不完全折减刚度法替代非线性显式后屈曲分析，作为非完善筒壳结构的承载分析方法，开展面向缺陷容忍的加筋筒壳快速优化设计；最后，针对优化结果进行不完全折减刚度法的精度验证，若不满足分析精度，则基于优化结果对刚度折减关系函数进行更新，若满足精度则可获得加筋筒壳的最终设计。其优化流程如图2所示。

图2 面向缺陷容忍的加筋筒壳快速优化设计框架Fig. 2 Flowchart of fast optimization design of stiffened cylindrical shell oriented to imperfection tolerance

3 算例说明

3.1 模型描述

本节以正置正交网格加筋筒壳结构为例，对所提出的面向缺陷容忍的加筋筒壳快速优化设计框架进行详细说明。其中，加筋筒壳的半径为800 mm，高度为1 000 mm。初始设计的具体尺寸参数包括：蒙皮厚度ts=1.6 mm，筋条厚度tr=2.5 mm，筋条高度hr=11.9 mm，环筋数量NC=17，纵筋数量NA=138，如图3所示。其材料为2A14铝合金，弹性模量为76 169 MPa，密度为2 700 kg/m3。有限元模型分析中的载荷工况和边界条件如下：筒壳下端为固支边界条件，在筒壳上端面圆心处建立一个参考点和上端面完全耦合，固定耦合点除轴向位移外的其他所有自由度，并在参考点处施加轴压载荷。

图3 正置正交网格加筋筒壳初始设计Fig. 3 Initial design of orthogonal grid stiffened cylinder

3.2 基于iRSM的非完善加筋筒壳承载分析

为了建立基于iRSM的非完善加筋筒壳承载快速分析模型，本文以基于假设凹坑几何缺陷的单点扰动载荷法（Single Perturbation Load Ap⁃proach， SPLA）为基准，采用优化的方式获得式（10）中的2个待定参数，从而确定刚度折减关系函数。

对于单点扰动载荷法，首先在筒壳外壁中部施加一个扰动载荷，以在筒壳上产生单个凹坑形式的几何缺陷。然后，将几何缺陷引入到筒壳完美的有限元模型中，并基于非线性显式后屈曲分析获得非完善筒壳结构的极限承载能力。随后，逐渐增加扰动载荷的大小，就可以获得如图4所示的几何缺陷敏感性分析曲线。可以看出，对于加筋筒壳初始设计，其极限承载能力随着扰动载荷的增加逐渐降低并趋于收敛。并且以扰动载荷为8 kN时，结构的极限承载作为承载下限和设计载荷，可得初始设计的设计载荷为3 209.3 kN。

图4 SPLA缺陷敏感性分析曲线Fig. 4 Imperfection sensitivity analysis curves by SPLA

然后，以上述SPLA 的分析结果为基准，采用如式（11）所示的优化列式，进行折减刚度函数参数的求解。需要说明的是，对于加筋筒壳结构，筋条的作用就是增加整体筒壳结构的抗弯刚度，因此本文只对蒙皮的膜刚度按照如式（10）所示的函数关系进行折减，对筋条刚度不作任何处理。

通过优化可得折减刚度关系函数的2个参数分别为A=51.39，B=1.30，而基于此刚度折减关系函数可得到缺陷敏感性分析曲线如图 5所示。可见2条曲线在收敛段基本重合，并且iRSM得到的设计载荷为3 205.6 kN，与非线性显式后屈曲分析结果相比，误差仅为−0.12%。而非线性显式后屈曲分析大约需要13.5 min，iRSM仅需要1 min左右，计算成本仅为非线性显式后屈曲分析的7.41%，这说明iRSM可实现非完善加筋筒壳结构的承载力高精度快速分析。另外需要说明的是，本文所有的分析计算都是在配置为Intel（R） Xeon（R） CPU E5-2670 v3 @ 2.30 GHz and 128 GB RAM的工作站上进行的。

3.3 iRSM误差评估

为了充分验证iRSM的分析精度，本节在如表1所示的加筋筒壳结构设计空间中，对iRSM进行整体的误差评估，并验证基于iRSM和非线性显式后屈曲分析的非完善筒壳结构承载在设计空间内整体趋势是否一致。

表1 加筋筒壳设计变量上下限Table1 Upper and lower limits of design variables for stiffened cylindrical shell

首先，在初始设计基础上分别单独改变不同设计变量的参数，进行误差评估。①蒙皮厚度ts分别取0.8、1.2、1.6、2.0、2.4 mm，而其他变量参数与初始设计完全一样，并分别使用iRSM和基于非线性显式后屈曲分析的SPLA进行缺陷敏感性分析，可得结果如图 5所示，其中当蒙皮厚度为2.4 mm时，iRSM的误差最大为4.25%；②筋 条 厚 度tr分 别 取1.25、1.88、2.5、3.13、3.75 mm，并使其他参数不变，可得分析结果如图 6所示，其中当筋条厚度为3.75 mm时，误差最大为2.58%；③筋条高度hr分别取8.93 、11.9、14.88 mm，并使其他参数不变，可得分析结果如图 7所示，其中当筋条厚度为14.88 mm时，误差最大为5.7%。

图5 不同蒙皮厚度参数下iRSM与SPLA结果对比Fig. 5 Comparison of iRSM and SPLA with different skin thickness

图6 不同筋条厚度参数下iRSM与SPLA结果对比Fig. 6 Comparison of iRSM and SPLA with different stiffener thickness

图7 不同筋条高度参数下iRSM与SPLA结果对比Fig. 7 Comparison of iRSM and SPLA with different stiffener heights

从分析和误差评估的结果可以看出，当在一定范围内改变单一设计变量的参数时，采用iRSM获得的非完善加筋筒壳结构承载能力具有较高的精度。并且尺寸变量中，筋条高度的变化对iRSM的误差影响最大。

然后，使用最优拉丁超立方在设计空间中进行采样，并使用不同方法对采样点分别进行分析，完成设计空间的整体误差评估。本文共采取20个样本点，每个样本点不同分析方法的承载结果如图8所示。其中，大部分样本点iRSM分析结果的误差在10%以内，尽管个别样本点的误差较大（−13.28%），但所有样本点在整体设计空间中仍具有相同的趋势。这说明基于iRSM的非完善筒壳承载力快速分析，可被用于面向缺陷容忍的加筋筒壳结构优化当中，替代非线性显式后屈曲分析。

图8 样本点非完善筒壳的屈曲载荷分布Fig. 8 Buckling load distribution of imperfect stiffened cylinders for sample points

3.4 面向缺陷容忍的加筋筒壳快速优化设计

接下来，以不超过加筋筒壳初始设计质量为约束条件，提高缺陷作用下的轴压极限承载能力为目标，开展面向缺陷容忍的轴压加筋筒壳优化设计，优化列式如下所示：

式中：PSPLA为加筋筒壳在几何凹坑缺陷作用下的极限承载能力；M为结构的质量；M0为初始设计的质量约束；Xil和Xiu分别为设计变量的下限值和上限值，如表1所示。并为了形成对比，基于非线性显式后屈曲分析和iRSM分别开展优化。此外，为了尽可能找到全局最优解，本文使用混杂优化策略，即先使用多岛遗传优化算法，再使用序列二次规划等梯度类算法开展优化。

其中，考虑到非线性显式后屈曲单次分析计算成本过高，多岛遗传算法的主要参数设置如下：种群数为30，岛数为4，最大遗传代数为10。从如图 9所示的多岛遗传优化迭代曲线可知，2个优化在全局优化过程中均已收敛。然后以遗传算法的最优解作为初始点，使用序列二次规划方法进行梯度优化，可得到最终优化结果。其中基于iRSM的优化结果为3 768.3 kN，使用非线性显式后屈曲分析进行验证，可得其SPLA设计载荷为3 630.5 kN，iRSM优化结果的误差为3.8%；而基于非线性显式后屈曲分析优化结果的SPLA设计载荷为3 672.2 kN，优化结果详情如表 2所示。

图9 多岛遗传优化迭代曲线Fig. 9 Optimization iteration curves of multi-island genetic

表2 面向缺陷容忍加筋筒壳结构的优化结果Table 2 Optimization results of stiffened cylindrical shells oriented to imperfection tolerance

首先可以看出，相比初始设计，2个面向缺陷容忍加筋筒壳优化设计结果的SPLA设计载荷分别提升13.12%和14.55%，其通过同时提高屈曲载荷和抗缺陷能力，实现了加筋筒壳设计载荷的提升，在保证轻量化水平的同时，提高了筒壳结构的安全裕度，他们的位移承载曲线如图 10所示。另外，2个优化结果的设计载荷近似但变量参数各不相同。这一方面是因为筒壳屈曲优化问题的设计空间本身就具有强非线性、多峰值等特点，其可能存在多个性能表现相近却存在明显差异的设计结果。另一方面是因为基于iRSM的非完善加筋筒壳快速承载分析仍存在微小误差导致的。其中可以发现，基于iRSM的优化结果倾向于增加筋条厚度，而基于非线性显式后屈曲的优化结果倾向于增加蒙皮厚度，这是因为我们在构建折减刚度模型的过程中只对蒙皮的膜刚度进行了折减，对筋条刚度不作任何处理。

图10 优化结果的位移载荷曲线Fig. 10 Load vs end-shortening curves of optimization results for stiffened cylinders

此外，可以发现基于iRSM优化结果的设计载荷仅比基于非线性显式后屈曲分析优化结果的设计载荷低1.24%，但其花费的计算成本仅约为后者的16.67%（其中单次分析均包括1.5 min左右的建模和划分网格的时间成本）。这说明仅通过将优化过程中的非线性显式后屈曲分析替换为不完全折减刚度法，就能在稳定找到优化结果的同时，减少80%以上的计算成本，实现面向缺陷容忍加筋筒壳结构的快速优化设计。另外，还可以通过刚度等效、模型降阶、连续分步类优化策略、创新优化算法等技术途径在此基础上进一步减少优化过程所需的计算成本。

4 结 论

针对航天大直径加筋筒壳结构分析效率低和优化设计成本高的问题，本文扩展了不完全折减刚度法（iRSM）在加筋筒壳上的应用，并基于凹坑形式的人工几何缺陷，以显式非线性分析方法为基准，验证了不完全折减刚度法的精度和加筋筒壳结构在设计空间中极限承载下限分布的一致性。同时建立了一种面向缺陷容忍的加筋筒壳快速优化设计框架，其使用iRSM替代非线性显式后屈曲分析，进行非完善加筋筒壳结构的极限承载能力计算和优化设计。算例结果显示，相比初始设计，优化结果可在质量不变的前提下，实现设计载荷提升10%以上，有效提高航天加筋筒壳结构的承载效率；并且所提出方法能在保证稳定找到优化结果的同时，降低计算成本80%以上，实现面向缺陷容忍加筋筒壳结构的高效优化设计，为中国航天主承力薄壁加筋筒壳结构的快速优化设计提供技术支撑。