运载系统运用工程专栏
加筋圆柱壳广泛应用于航天运载器结构当中[1-2]。当加筋圆柱壳的轴压载荷过大时,结构会突然塌陷,局部出现较大的径向变形,瞬间丧失轴压承载能力,即发生所谓的屈曲。因此需在设计时准确预估结构的屈曲临界载荷[3-4]。然而,此类结构对于几何、边界的微小缺陷非常敏感,其实际承载能力往往显著低于理论值[5]。NASA的设计指南SP-8007[6]基于大量的试验数据,建议采用线性解析解乘以1个0~1之间的折减因子(knockdown factor,KDF)的方式计算考虑缺陷的结构屈曲载荷,其中KDF为结构实际屈曲载荷与线性解析解的比值。这种方法十分便捷,非常适合在初始设计阶段各种参数频繁迭代时使用。
随着工艺、材料以及制造技术的进步,早期SP-8007给出的KDF较小而显得过于保守,导致了结构设计重量的冗余[7-8]。人们寻求新的、更合适的KDF选用方法,期望由此能大幅度降低结构重量。代表性的有:NASA的SBKF计划[9]和欧航局的DESICOS计划[10]等,研究人员尝试用少量的试验验证非线性有限元模型的准确性。早期的KDF是通过总结试验结果获得的,因此可称为“基于试验的KDF”(test based KDF)。近年来可以通过高仿真度有限元来获得KDF,可称为“基于分析的KDF”(analysis based KDF)[11]。这种方法成本低、效率高、适用范围广,是目前KDF研究的主流方法。其核心思想是在有限元模型中引入几何或边界缺陷,得到可替代试验的数值仿真结果,而数值解与解析解的比值即为“基于分析的KDF”。典型的有:实际缺陷法[12-13]、模态缺陷法[14]、单点扰动载荷法[15]、多点扰动载荷法 [16]以及边界扰动载荷法[17-18] 等。其中,边界扰动载荷法(single boundary perturbation approach,SBPA)是在壳体模型加载端添加1个凸台产生边界扰动模拟加载缺陷,在多次有限元分析中不断增加凸台的高度,屈曲载荷逐渐下降并最终收敛于一个“下限值”,此时定义该“下限值”为SBPA的极限屈曲载荷。相较于试验值,SBPA的计算结果略微保守[19-20],适合于工程应用。
在航天运载器的初始设计阶段,结构参数会频繁地迭代、优化,需要快速、准确地计算出结构的承载能力。如果每次结构参数发生变化都进行一次非线性有限元分析,则响应速度慢、计算成本高。虽然采用解析解乘以固定值的KDF方法便捷,但却忽略了KDF随设计参数发生的变化,可能导致计算结果失真。本文基于理论解析法和SBPA法计算了正置正交网格加筋圆柱壳的KDF,旨在总结出蒙皮厚度、壳高以及筋条宽度等参数变化时其相应的KDF变化规律。从而,可以在初始设计阶段基于若干次非线性有限元的分析结果而迅速地线性插值出针对某一组新设计参数的KDF。这种动态的KDF显然会比固定值的KDF更准确。
基于分析的KDF可定义为
(1)
式(1)中,Nlb和Nsb分别是基于解析法和SBPA法得到的结构屈曲载荷。解析法是针对完美结构的线性解,而SBPA法是针对包含严重初始缺陷的结构的非线性解。因此,KDF也可以看作是对结构缺陷敏感性的定量描述。
加筋圆柱壳轴压屈曲载荷的解析解可依据NASA SP-8007[6]手册中的公式计算,具体为
(2)
其中,
(3)
(4)
(5)
(6)
(7)
(8)
式(2)—式(8)中:Nx为沿厚度方向的线屈曲载荷;分别为加筋圆柱壳沿轴向和环向的弹性模量;为剪切模量;分别为加筋圆柱壳沿轴向和环向的弯曲刚度;为沿厚度方向的弯曲刚度;为加筋圆柱壳的耦合常数,具体内容详见文献[6],此处不再赘述。
图1给出了在有限元软件中实现SBPA的方法。首先,在结构加载端边缘建立1个高度为h、圆周角为α的凸台,产生边界扰动缺陷,然后,在结构顶端附近建立1个刚性平面并使其沿轴向向下移动一小段距离,然后进行非线性分析。在计算结果中读取载荷位移曲线的峰值作为结构的屈曲载荷。如图1所示,逐渐增大边界扰动高度并进行多次非线性分析,所得到的结构屈曲载荷会随着边界扰动高度的增加而逐渐降低,且当边界扰动高度超过某个阈值h1时,结构屈曲载荷不再发生显著变化,边界扰动高度h1所对应的屈曲载荷N1即为SBPA法得到含缺陷的结构产品的屈曲载荷Nsb。
图1 SBPA原理示意图
Fig.1 Schematic mechanism of SBPA
本节以1个典型网格加筋壳(编号TA01[21])为例,采用SBPA法计算了其非线性、含缺陷状态的屈曲载荷,介绍了主要参数的选取原则。TA01的结构参数及解析解如表1所示。
表1 网格加筋壳TA01参数
Table 1 The parameters of TA01
参数数值壳高度L/mm1 828.8壳半径R/mm1 219.2蒙皮厚度t/mm2.54纵筋数目75环筋数目17纵环筋横截面尺寸b×h/mm22.54×7.62密度/(t·mm-3)2.7×10-9弹性模量/MPa70 000泊松比0.3解析解/MN3.186
TA01的有限元模型如图2所示,在加载端建立1个凸台模拟边界缺陷。凸台厚度与蒙皮一致,宽度设为圆周角2°(R/t>200宽2°;R/t<200宽4°),高度由后面的收敛性考察确定。
图2 TA01的有限元模型
Fig.2 Finite element model of TA01
边界条件与网格划分如图3所示,在结构顶端附近,建立1个解析刚体平面并将其与结构顶端的蒙皮边缘和凸台边缘之间设置为无摩擦硬接触。刚体平面以及壳体顶端边缘均约束了除轴向平动之外的其他全部自由度,并给刚体平面沿轴向施加向下的强制位移10 mm。圆柱壳底端边界条件为固支,约束了全部自由度。蒙皮、筋条、以及凸台均采用四节点减缩积分壳单元(S4R)模拟。通过网格收敛性考察,确定有限元仿真的网格单元尺寸为15 mm,结构共包含 84 525个网格,83 850个节点。整个计算过程在ABAQUS中完成,采用了含人工阻尼的非线性静力分析步进行计算,有限元分析参数如表2所示。
表2 有限元分析参数
Table 2 Finite element analysis parameters
参数数值参数数值阻尼因子1×10-5初始增量1×10-2最小增量1×10-10最大增量1×10-2最大增量数5×102
图3 边界条件与网格划分
Fig.3 Boundary conditions and meshing
加筋圆柱壳发生屈曲时,其轴压承载力会骤然降低,结构坍塌,产生大量变形,图4给出了凸台高度为0 mm和1 mm时的载荷位移曲线以及屈曲模态。其中完美结构的载荷位移曲线呈明显的线性特征直至峰值,峰值屈曲载荷为2.807×106 N,该阶段,完美结构为轴对称的环状模态。峰值后,曲线陡然下降,模态由峰值处的环状轴对称模态转变为沿环向分布的波形模态并逐渐扩展,最终在后屈曲阶段形成吉村模态[21]。
图4 TA01载荷位移曲线以及屈曲模态
Fig.4 Load-displacement curve for TA01 with deformation plots
对于含边界缺陷的结构,加载面先接触凸台的边缘,此时其载荷位移曲线为一小段斜率较低的直线。当加载面接触到上端面之后,载荷位移曲线变为一段斜率较大的直线,对应的结构变形为在凸台下方出现菱形凹陷。载荷在达到峰值1.654×106 N后迅速下降,凹陷开始向两侧传播,最终形成吉村模态。虽然其最终模态和完美结构相同,但屈曲过程和屈曲载荷有显著的区别。这是因为,边界缺陷诱导结构在较低载荷时即发生局部凹陷,而当局部凹陷开始扩散时,结构即损失了大部分的承载能力。
在SBPA分析中, 阻尼因子和边界上凸台高度的取值对计算结果影响显著。结构发生屈曲的瞬间是一个动态过程, 用准静态分析去模拟这个动态过程时需要结构阻尼的存在, 以降低结构的加速度, 从而保证分析能够收敛。阻尼因子取值过小,分析不收敛, 取值过高又会导致分析结果失真。建议在SBPA分析中, 考察屈曲载荷关于阻尼系数的收敛性。二者的关系如图5所示,由图5可以看出,随着阻尼系数的降低,结构屈曲载荷逐渐减小,而当阻尼系数达到1×10-5时,屈曲载荷收敛,因此本例中1×10-5即为合适的阻尼因子。
图5 阻尼因子收敛性研究
Fig.5 The convergence study of damping factor
图6给出了结构的屈曲载荷随边界凸台高度的变化曲线。
图6 屈曲载荷与凸台高度的关系曲线
Fig.6 The curve of buckling load vs.shim height
从图6可以看出,随着扰动高度的不断增加,结构屈曲载荷逐渐减小并趋于收敛,当扰动高度增加至1 mm时,结构屈曲载荷不再发生显著变化且收敛于一个“下限值”N1,此时N1=1.654×106 N,即为SBPA解。需要注意的是,这个解所对应的是含有非常恶劣缺陷的结构的承载能力,从工程实际的角度看是偏于保守的。结合解析解以及SBPA解,可得到基于分析的KDF为0.52。
初始设计阶段结构参数需要频繁的迭代调整,但一般不会出现很大的改动,因此本文仅考虑各个参数在正负20%内变化时对KDF的影响。本节总结出了KDF随纵筋截面宽度、蒙皮厚度以及结构高度等主要几何参数的变化规律。
以TA01算例为基础,其他条件保持不变,令纵筋横截面面积在20%范围内变化,计算其不同尺寸下的KDF,结果列于表3中。
表3 改变纵筋截面面积的计算结果
Table 3 The results of various stringer section widths
筋条横截面积尺寸b×h/mm解析解/MNSBPA/MNKDFA-20%2.032×6.0962.5231.3890.55A-10%2.286×6.8582.8321.5200.54TA012.54×7.623.1861.6540.52A+10%2.794×8.3823.5751.8280.51A+20%3.048×9.1444.0001.9870.50
由表3可知,随着纵筋截面面积的增加——结构的材料用量在增加——屈曲载荷的解析解和SBPA解都随之增加,这是符合工程直觉的。图7给出了纵筋截面宽度对KDF的影响曲线。
图7 纵筋截面面积对KDF的影响
Fig.7 The stringer area of section effect on KDF
由图7可知,随着纵筋截面面积变大,结构的KDF呈现出约10%的线性下降,说明纵筋截面面积对KDF值具有显著的影响。另外,选用具有较小截面的筋,结构会具有更高的KDF,即更低的缺陷敏感性。因此,在注意回避筋条局部失稳的情况下应尽可能地选用低而薄的纵筋。
本节在其他条件保持不变的情况下,令蒙皮厚度在20%范围内变化,计算其不同尺寸下的KDF,结果列于表4中。
表4 改变蒙皮厚度的计算结果
Table 4 The results of various skin thicknesses
蒙皮厚度尺寸/mm解析解/MNSBPA/MNKDFt-20%2.032.5811.2890.50t-10%2.2862.8731.4560.51t2.5403.1861.6540.52t+10%2.7943.5241.8840.53t+20%3.0503.8932.1230.54
蒙皮厚度对KDF的影响如图8所示。由于蒙皮的厚度增加意味着材料用量的增加,因此无论是线性解析解还是SBPA解都随之增加。由图8中KDF随蒙皮厚度的变化趋势可以看出,随着蒙皮厚度的增加,结构的KDF呈现出约7%左右的线性递增,说明蒙皮厚度的变化对KDF值有显著影响,且蒙皮厚度越大结构的边界缺陷敏感性越低。
图8 蒙皮厚度对KDF的影响
Fig.8 The skin thickness effect on KDF
在保持其他条件保持不变的情况下,令壳体高度在20%范围内变化,计算其不同尺寸下的KDF,结果列于表5中。高度对KDF的影响如图9所示。
表5 改变高度的计算结果
Table 5 The results of various heights
壳高尺寸/mm解析解/MNSBPA/MNKDFL-20%1 463.03.3241.7720.53L-10%1 645.93.2501.7080.53L1 828.83.1861.6540.52L+10%2 011.73.1241.6170.52L+20%2 194.63.0681.5840.52
图9 高度对KDF的影响
Fig.9 The height effect on KDF
由表5可知,随着高度的增加,结构屈曲载荷的解析解和SBPA解小幅降低。由图9中的KDF变化趋势可以看出,随着高度的增加,结构KDF仅呈现出2%的变化,说明高度的改变并不会影响结构的缺陷敏感性。
本文中基于SBPA法,对网格加筋圆柱壳的KDF随几何参数变化的规律展开研究,研究结果表明:
1) 结构的KDF与壳体高度几乎无关,但与蒙皮厚度线性正相关、与纵筋截面面积线性负相关。
2) 在初始设计阶段,可以尽量选用较厚的蒙皮和低而薄的纵筋以提高结构的抗缺陷性能,并基于若干次非线性有限元的分析结果而迅速插值出针对某一组新设计参数的KDF,从而大幅度提高计算效率。
[1]WAGNER H,HÜHNE C,ELISHAKOFF I.Probabilistic and deterministic lower-bound design benchmarks for cylindrical shells under axial compression[J].Thin-Walled Structures,2020,146:106451.
[2]WAGNER H,HÜHNE C,NIEMANN S.Buckling of launch-vehicle cylinders under axial compression:a comparison of experimental and numerical knockdown factors[J].Thin-Walled Structures,2020,155:106931.
[3]MA X,HAO P,WANG F,et al.Incomplete reduced stiffness method for imperfection sensitivity of cylindrical shells[J].Thin-Walled Structures,2020,157:107148.
[4]JIAO P,CHHEN Z,MA H,et al.Buckling behaviors of thin-walled cylindrical shells under localized axial compression loads,Part 2:Numerical study[J].Thin-Walled Structures,2021,169:108330.
[5]JIAO P,CHEN Z,MA H,et al.Buckling behaviors of thin-walled cylindrical shells under localized axial compression loads,Part 1:Experimental study[J].Thin-Walled Structures,2021,166:108118.
[6]PETERSON J P,SEIDE P,WEINGARTEN V I.NASA SP-8007,buckling of thin-walled circular cylinders[S].NASA Space Vehicle Design Criteria,1965.
[7]王博,郝鹏,田阔.加筋薄壳结构分析与优化设计研究进展[J].计算力学学报.2019,36(1):1-12.WANG Bo,HAO Peng,TIAN Kuo.Recent advances in structural analysis and optimization of stiffened shells[J].Chinese Journal of Computational Mechanics,2019,36(1):1-12.
[8]焦鹏.局部轴压下薄壁圆柱壳结构的屈曲行为及设计方法研究[D].杭州:浙江大学,2021.JIAO Peng.Research on the buckling behaviors and design method of thin walled cylindrical shell structures under localized axial compression loads[D].Hangzhou:Zhejiang University,2021.
[9]HILBURGER M.Developing the next generation shell buckling design factors and technologies[C]//53rd AIAA/ASME/ASCE/AHS/ASC structures,structural dynamics and materials conference 20th AIAA/ASME/AHS adaptive structures conference 14th AIAA.2012:1686.2012.
[10]DEGENHARDT R.New design scenario for future composite launcher structures[C]//Proceedings of the 30th Congress of the International council of the Aeronautical Sciences (ICAS).2016.
[11]PAN L C,WANG D F.Buckling of thin-walled cylindrical shells of desulphurizing tower under wind loading[J].Applied Mechanics and Materials,2014,662:147-152.
[12]HILBURGER M W,LINDEL M C,WATERS W A,et al.Test and analysis of buckling-critical stiffened metallic launch vehicle cylinders[C]//2018 AIAA/ASCE/AHS/ASC structures,structural dynamics,and materials conference.2018:1697.
[13]SCHULTZ M R,SLEIGHT D W,GARDNER N W,et al.Test and analysis of a buckling-critical large-scale sandwich composite cylinder[C]//2018 AIAA/ASCE/AHS/ASC Structures,Structural Dynamics,and Materials Conference.2018:1693.
[14]European Committee for Standardization.Strength and stability of shell structures:EN 1993-1-6[S].EN Special Publication,2007.
[15]HüHNE C,ROLFES R,BREITBACH E,et al.Robust design of composite cylindrical shells under axial compression—Simulation and validation[J].Thin-Walled Structures,2008,46(7/9):947-962.
[16]HAO P,WANG B,LI G,et al.Worst multiple perturbation load approach of stiffened shells with and without cutouts for improved knockdown factors[J].Thin-Walled Structures,2014,82(9):321-330.
[17]WAGNER H,PETERSEN E,KHAKIMOVA R,et al.Buckling analysis of an imperfection-insensitive hybrid composite cylinder under axial compression-numerical simulation,destructive and non-destructive experimental testing[J].Composite Structures,2019,225:111152.
[18]WAGNER H,HÜHNE C,ROHWER K,et al.Stimulating the realistic worst case buckling scenario of axially compressed unstiffened cylindrical composite shells[J].Composite Structures,2017,160:1095-1104.
[19]WAGNER H,HÜHNE C,JANSSEN M.Buckling of cylindrical shells under axial compression with loading imperfections:An experimental and numerical campaign on low knockdown factors[J].Thin-Walled Structures,2020,151:106764.
[20]WAGNER H,HÜHNE C,NIEMANN S.Buckling of launch-vehicle cylinders under axial compression:A comparison of experimental and numerical knockdown factors[J].Thin-Walled Structures,2020,155:106931.
[21]HILBURGER M,LOVEJOY A,THORNBURGH R,et al.Design and analysis of subscale and full-scale buckling-critical cylinders for launch vehicle technology development[C]//53rd AIAA/ASME/ASCE/AHS/ASC Structures,Structural Dynamics and Materials Conference 20th AIAA/ASME/AHS Adaptive Structures Conference 14th AIAA.2012:1865.