基于弧长法的南瓜型超压气球非线性后屈曲分析

2022-09-29卜亚楼杨燕初蔡榕赵荣张向强

科学技术与工程 2022年22期

卜亚楼，杨燕初*，蔡榕，赵荣，张向强

(1.中国科学院大学航空宇航学院，北京 100040； 2.中国科学院空天信息创新研究院，北京 100094)

南瓜型超压气球是目前长航时浮空飞行平台的主流设计形式。其通过在欧拉体球形[1]表面均匀铺设加强筋、使球膜形成局部高曲率鼓包的方式来减小球膜应力，提高承压能力。因形状似南瓜，因此被称作南瓜型超压气球。南瓜球具有驻空高度稳定、飞行时间长、效费比高的特点，是执行临近空间科学试验和高空飞行任务最稳定的平台之一[2]。南瓜球虽然在承压能力和应力分布方面具有明显优势，但在某些情况下却呈现出扭曲、不对称的构型，使得南瓜球的稳定性研究和分析成为其结构设计的重要问题。

针对南瓜球的失稳现象，Calladine[3]最早从几何稳定性的角度开展研究。通过将球膜的几何刚度与欧拉梁的抗弯刚度进行类比，得到球膜幅数n与凸出角α(球膜局部凸出鼓包的弧长所对应的圆心角)之间的稳定性判据。但在随后的研究中发现，Calladine[3]的结论更适用于按照恒定凸出角度(constant bulge angle, CA)设计的气球。Smith等[4]认为按照恒定凸出半径(constant bulge radius, CR)设计的气球展开稳定性更好；而Schur[5]通过地面试验推断，过量的幅宽可能导致南瓜球的展开不稳定性，并指出按照恒定凸出半径设计的气球更有利于展开；Baginski等[6-9]通过能量最小原理计算气球的平衡构型。对成型气球进行稳定性分析发现，如果气球处于不稳定的平衡状态，那么在正常上升过程中可能不会完全展开。

有限元技术的发展使得利用数值方法来求解结构的失稳问题成为可能。文献[10-11]使用基于动态松弛求解方法的内部有限元程序套件inTENS来模拟南瓜球的失稳行为。该方法基于显式求解技术，其时间步动态分析结合了自动动能阻尼控制技术，通过给模型设定内部压力和相应的薄膜应力进行数值预测。文献[10-11]使用该程序套件对直径为10 m的南瓜球进行求解，得到了内部压力为690 Pa时的扭曲形态，数值预测结果与试验相比有良好的相关性；Deng[12]基于ABAQUS/Explicit计算了部分充气气球的平衡构型。其中考虑了球膜材料的非线性行为，对正交各向异性黏弹性材料进行了分析和建模，并基于变泊松比的起皱准则模拟了球膜的褶皱行为，获得的仿真结果与试验相比有较好的一致性。Deng[12]还在仿真的基础上结合量纲分析构建了南瓜球稳定性判据的经验公式，是南瓜球稳定性研究上的一大进展；Pagitz等[13]认为南瓜球扭曲的稳定平衡构型是一个涉及平衡路径分叉屈曲的问题，他们通过求解整体切线刚度矩阵的特征值和特征向量来研究失稳压强和失稳形态；Xu等[14-15]利用ABAQUS软件开展特征值屈曲分析，预测了气球的临界失稳压强和相应的屈曲模态，对可能影响南瓜球稳定性的因素进行灵敏度分析，获得了重要的研究结论；之后在特征值屈曲分析的基础上基于*static, stabilize开展非线性后屈曲分析以模拟气球超过临界压力后的响应。

Koiter[16]提出了考虑初始缺陷的非完善结构稳定性的一般准则和缺陷敏感度的概念，指出实际结构不可避免地存在初始缺陷，这会降低结构的屈曲载荷。具体地，Baginski等[7]指出南瓜球制造的不精确性、球膜和加强筋材料的不均匀性，以及外界环境的变化(如温度)对稳定性有很大影响。这表明，初始缺陷对南瓜球稳定性的影响不可忽略。Xu等[15]中虽构造了初始缺陷开展非线性后屈曲分析，但其侧重点在于研究气球超过临界压力后的响应，无法获得缺陷模式、缺陷幅值等对南瓜球极限承压能力的影响。因此，目前对于初始缺陷如何影响南瓜球稳定性的问题仅停留在定性认识上，尚未开展深入研究。

较多的薄膜结构分析采用显式准静态算法[17]，而南瓜球由于整体结构尺寸较大、球膜具有较强的几何非线性且球膜和加强筋弹性模量相差悬殊，传统计算方法面临着计算量大、计算不易收敛的问题。更重要的是，显式准静态算法无法获得缺陷模式、缺陷幅值等对南瓜球极限承压能力的影响。而基于弧长法的非线性后屈曲分析是研究结构缺陷敏感性的重要方法，其在追踪失稳路径、捕捉极限承压力方面具有明显优势。文献[18]采用弧长法对膜结构雷达罩进行了稳定性研究，取得了较好的结果，是弧长法在膜结构领域的最新应用。此外，该方法在球壳受外压屈曲[19-20]、网壳结构屈曲[21-22]、复合材料加筋壁板受轴压屈曲[23-25]等领域也有广泛应用。

因此，在前人研究的基础上，现将弧长法应用到南瓜球稳定性研究上。基于该方法既能获得南瓜球的缺陷敏感特性，也能获得南瓜球前后屈曲的全部信息,是南瓜球稳定性研究的一种新的方法。采用弧长法对受内压南瓜球的非线性屈曲过程进行数值模拟，基于特征值屈曲模态位移构建了两种缺陷模式，对模态缺陷下受压南瓜球的稳定性进行研究。通过对比不同缺陷幅值和缺陷模式下南瓜球的临界载荷，研究南瓜球对几何缺陷的敏感性。针对弧长法下降段难以收敛的问题，使用重启动方法结合显式动力学计算南瓜球屈曲后的响应，获得南瓜球的失稳形态和应力分布，为南瓜球在工程领域的设计制造提供参考。

1 稳定性分析的有限元法

1.1 线性特征值屈曲分析

特征值屈曲分析是一种线性摄动分析，它基于小位移线弹性理论，可以计算包含桁架单元、梁单元、板壳单元、实体单元的结构临界载荷系数和相应的屈曲模态。结构的静力平衡方程为

K+λKGu=Q

(1)

式(1)中：K为结构的弹性刚度矩阵；KG为结构的几何刚度矩阵；λ为特征值；u为结构的模态特征向量；Q为作用在结构上的载荷。

平衡方程[式(1)]失稳的条件是方程存在奇异性，即等效刚度矩阵的行列式的值为零。

|K+λKG|=0

(2)

线性屈曲分析即为求解式(2)的特征值，所求得的特征值即为临界载荷系数。将某一阶屈曲模态特征值与施加的计算载荷相乘，就得到该阶模态下结构的特征值屈曲临界载荷。u为结构的模态特征向量，即屈曲模态位移，它预测了结构可能的失效形式。由于线性特征值屈曲以小位移小应变的线弹性理论为基础，没有考虑结构在受载过程中结构构型的变化，在外力施加的各个阶段，平衡方程始终建立在初始构型上，忽略了实际加载过程中结构刚度矩阵的变化。因此，线性特征值屈曲分析计算结果偏保守，对于缺陷敏感结构需要进行非线性后屈曲分析。

1.2 弧长法

弧长法是结构非线性分析中计算稳定高效的一种数值方法，能够有效地分析结构非线性前后屈曲并跟踪后屈曲路径，从而获得结构失稳前后的全部信息，并有效捕捉极限载荷。

在非线性后屈曲分析中，针对Newton-Raphson迭代方法在屈曲临界点附近结构刚度接近奇异时迭代不易收敛的情况，通过结合修正的弧长法可以稳定求解通过极值点的结构分析全过程。假设载荷成比例加载，迭代方程为

(3)

将位移增量ΔUi与载荷比例因子λt+Δt通过弧长半径L联系起来可得到求解λt+Δt的约束方程，即

ω(Δλi)2‖PW‖2+β‖ΔUi‖2=(ΔL)2

(4)

式(4)中：‖PW‖、‖ΔUi‖为外载荷向量和位移增量的二范数；Δλi为第i次迭代的载荷比例因子增量；ΔL为弧长增量的半径；ω、β为尺度因子，在不同的弧长控制方法中取值不同。

2 有限元建模与特征值屈曲分析

2.1 数值模型建模

取某一南瓜球作为算例，其几何参数如表1所示。

表1 南瓜球几何参数

图1为赤道平面处相邻加强筋之间球膜的截面形状。气球的加强筋位于相邻两幅球膜的焊缝处，气球的赤道半径用R表示。对于n幅球膜的情况，赤道对角为θ=2π/n。

s为鼓包的弧长；c为弦长；r为凸出半径；α为凸出角图1 赤道平面处的凸出截面形状Fig.1 Lobe geometry at equator

根据几何关系，弦长c可用赤道半径和赤道对角表示：

(5)

此外，弦长也可用凸出半径和凸出角表示：

(6)

根据式(5)和式(6)，可得出赤道半径与凸出半径的关系为

(7)

鼓包弧长为

s=αr

(8)

由于鼓包的圆弧与凸出角或凸出半径有关，目前南瓜球的设计主要采用两种方式，一是凸出半径沿子午线保持恒定的凸出半径设计，二是凸出角沿子午线保持恒定的凸出角度设计。祝榕辰[26]证明了恒定凸出半径的设计方法能使球膜表面的应力分布更加均匀；文献[5]从展开稳定性的角度论证了恒定凸出半径设计的气球具有更好的稳定性。因此，恒定凸出半径设计方法得到了较广泛应用。本文的计算模型采用恒定凸出半径设计。根据以上几何参数及几何关系完成建模，三维模型图如图2所示。

图2 南瓜球三维模型示意图Fig.2 Schematic layout of a 3D pumpkin balloon

气球球膜和加强筋对其所用的材料有多种特性要求，随着材料科学的发展，球膜和加强筋的特性在不断改善。线性低密度聚乙烯以其优越的性质取代了以前的低密度聚乙烯作为球膜材料；PBO(p-phenylene benzobisoxazole)以其优越的力学性能、热稳定性和优异的强度重量比被选择作为加强筋的材料。目前，世界上几个开展高空气球项目的国家均在采用这种材料，如ULDB(ultra long duration balloon)项目。具体材料属性如表2所示。

南瓜球的线弹性各向同性球膜采用M3D4R膜单元建模。这是由于膜单元仅有平面内刚度，没有弯曲刚度，且M3D4R是一个4节点四边形单元，为减缩积分单元，具有计算高效的优点；桁架单元通常用于模拟仅在轴向承受拉伸或压缩应力的细长结构，气球中的加强筋虽然无法承受压应力，但在超压状态下所承受的拉应力与桁架单元沿轴向传递的拉伸特性是相同的，因此采用桁架单元T3D2对加强筋进行建模，其中加强筋的刚度定义为弹性模量和横截面的乘积。

南瓜球超压时，载荷会从球膜转移到位于焊缝处的加强筋上，并连接到位于气球顶部和底部的刚性法兰盘上。位于焊缝处的加强筋通过封套约束在球膜表面，加强筋和球膜之间沿轴向为滑动摩擦约束。但实际数值仿真中，由于筋膜之间的滑动摩擦约束关系不易建模，且筋膜之间容易相互穿透导致计算不收敛，因此通常将筋膜建模为共节点约束以代替滑动摩擦约束来传递轴向力。赵荣等[27]已通过数值仿真证明了加强筋和膜单元之间采用共节点约束与采用滑动摩擦约束对球膜承力并无很大差别；文献[28]同样采用共节点约束开展静力学分析，取得了较为合理的结果。因此，为了建模简单起见，本研究同样采用共节点约束。顶部法兰盘和底部法兰盘均为全约束。此外，为了获得较好的网格质量，本分析中没有构建法兰盘部件，而是以桁架单元代替法兰，这样做可以提高计算的收敛性而不影响计算结果。单元数量对于模拟结果的准确度影响很大，需要进行无关性验证，一般认为单元数量的选取应该使前后两次分析结果相对浮动量小于等于1%。经计算比较，最后选择的划分单元数量为T3D2：13 600个，M3D4R：64 600个。取单位载荷1 Pa为计算载荷，以均布力形式施加在气球内表面。

表2 材料属性

2.2 线性特征值屈曲分析

首先对南瓜球进行线性特征值屈曲分析，得到屈曲压强和模态位移。在非线性后屈曲分析中，需要用模态位移构造初始几何缺陷。

计算得到前8阶屈曲压强和屈曲模态，结果如表3所示。从表3中看出，在前8阶屈曲模态中存在着模态形状相似、特征压强相近的4组相邻模态，如1阶和2阶、3阶和4阶，其屈曲模态构型除了在位移幅值上有微小的差别外，相近特征值下的两个模态为绕竖直轴旋转约15°(图3)，说明即使在相等的外载荷作用下，南瓜球仍存在多种可能的屈曲位移趋势，证明南瓜球具有分支屈曲的特性。

图4显示了南瓜球前4组相近特征值对应的屈曲模态。需要注意的是屈曲模态u为归一化向量，最大位移分量为1.0，这并不代表在临界载荷下变形的实际大小，但能够表明结构的可能失稳模式。屈曲模态构型表现出整体对称的n上n下(n为自然数)形态。高阶模态对应更大的n值；随着n的增大，需要更大的临界压强通过分叉加载点，但临界压强与n值并非简单的线性关系。n不超过4时，相邻特征压强相差较小；当n大于4时，相邻特征压强表现出显著的差距。

表3 南瓜球特征值屈曲分析结果

图3 具有相近特征值的两个屈曲模态对比图Fig.3 Comparison diagrams of two buckling modes with similar eigenvalues

图4 前四组特征模态的总位移等值线图Fig.4 Contourplots of total displacements of the first four eigenmodes

3 非线性后屈曲分析

3.1 初始缺陷的建立

气球受压时，在屈曲前主要受均匀的面内力，保持其球形均匀扩张的平衡状态；到达屈曲临界点时，有多个可能的失稳分支路径；当受到扰动时，将沿某一路径发生变形失稳，进入到新的平衡状态，属于分支点屈曲问题。基于线性屈曲分析的模态位移构造初始缺陷分布用以诱导结构的后屈曲行为，利用弧长法进行非线性后屈曲分析，获得南瓜球的缺陷敏感特性。

其中应用式(9)来表述缺陷几何的扰动网格：

(9)

式(9)中：ψi为第i阶模态形状；θi为相应的缩放因子，即缺陷比例系数，表示对模态位移进行缩放；Δxi为缩放后的模态位移。

由式(9)可知，构建初始缺陷的方式是多样的，既可以对某一阶模态位移进行缩放，也可以对多阶模态位移缩放后进行线性组合。因此在进行非线性后屈曲分析时，采用了两种方法来构造初始缺陷。第一种是临界缺陷模态法，即基于4上4下临界模态位移的缺陷气球。该理论认为，临界屈曲模态对应最小势能状态，当结构缺陷分布形式与临界模态吻合时，结构最容易发生屈曲，即最低阶模态缺陷是最不利的几何缺陷构型。运用临界缺陷模态法对模型进行分析时，先后引入5种缺陷幅值进行对比分析，最小缺陷比例系数为0.01。

图5显示了带有1%初始缺陷的气球的几何形状，可以看出，气球轻微变形，但足以触发后屈曲行为。由于缺陷均是基于临界模态位移，所以5种缺陷幅值所对应的构型一致，只是初始挠度不同。

图5 基于临界屈曲模态的初始缺陷气球Fig.5 Balloon with initial imperfection based on first buckling eigenmode

另一种是组合缺陷模态法。Triantafyllidis等[29]发现，在一些不稳定结构中，由于其特殊的几何形状，在相同的载荷水平下会出现整体屈曲和局部屈曲。对于这些结构，不稳定点对初始缺陷的形状非常敏感。该理论认为，不止一阶屈曲模态用来构造初始缺陷，多阶屈曲模态的组合也可用来构造几何缺陷。因此，本文的第二种初始缺陷的构造方法是使用前8阶模态位移进行线性组合来触发后屈曲行为，其中每种模态的缺陷比例系数均为0.01。

3.2 缺陷敏感度分析

采用弧长法对缺陷气球进行非线性有限元分析。图6为某一缺陷幅值下载荷和应变能随迭代弧长变化的曲线。当结构失稳时，载荷-位移曲线出现负刚度，对应的载荷-弧长曲线出现下降段，结构同时释放应变能维持平衡。由图6可知，载荷和应变能在同一时刻发生转折，说明结构发生失稳。

图6 载荷和应变能随迭代弧长变化曲线Fig.6 Curve of load and strain energy with iteration arc length

图7和图8分别为临界缺陷模态法和组合缺陷模态法在5种缺陷幅值下气球的载荷-弧长曲线。屈曲载荷结果对比如表4所示。

图7 基于临界缺陷模态法不同缺陷幅值下的载荷- 弧长曲线Fig.7 Load-arc length curves under different defect amplitudes based on critical defect mode method

图8 基于组合缺陷模态法不同缺陷幅值下的载荷- 弧长曲线Fig.8 Load-arc length curves under different defect amplitudes based on combined defect mode method

由表4可知，缺陷状态下非线性屈曲临界载荷小于线性特征值屈曲临界载荷，特征值屈曲分析结果偏保守。因此，工程实际中有必要对南瓜球进行非线性后屈曲分析。

表4 屈曲载荷结果对比

对比相同缺陷模式下不同缺陷幅值所对应的屈曲载荷发现，临界载荷随缺陷幅值的增加而减小。假设实际制造的南瓜球其初始缺陷符合临界屈曲模态构型1%的情况，其屈曲临界载荷为224.5 Pa。当缺陷比例系数为0.02时，其屈曲临界载荷为215 Pa，承载力相对下降约4.2%。当缺陷比例系数为0.05时，承载力相对下降13%，说明几何缺陷幅值对屈曲临界载荷有显著影响。

图9 充气过程中位移等高线图Fig.9 Contours of displacement during inflation

此外，缺陷幅值相同时，不同缺陷模式对应的非线性临界屈曲载荷各不相同。基于组合特征缺陷模态法构造的几何缺陷所对应的屈曲载荷要低于同幅值下的临界缺陷模态法。同样地，假设实际制造的南瓜球其初始缺陷符合组合屈曲模态构型1%的情况，其屈曲临界载荷为212 Pa。当缺陷比例系数为0.02时，其屈曲临界载荷为193.5 Pa，承载力相对下降约8.7%，当缺陷比例系数为0.05时，承载力相对下降53.5%。下降幅值高于相同缺陷比例系数下的临界缺陷模态情况，且随着比例系数的增大，结构承压力大幅度下降，这表明：①缺陷幅值和缺陷模式对失稳载荷有显著影响，南瓜球是一种缺陷敏感性结构；②第一阶模态缺陷不一定是最不利的几何缺陷构型，求解的临界载荷有可能高估结构的失稳载荷；组合模态缺陷也有出现的可能，其对应的屈曲临界载荷有可能更低，对南瓜球屈曲特性的研究应考虑多种模态缺陷情况，以便找出最低承载力对应的缺陷模式。

3.3 失稳形态和应力分布分析

针对弧长法下降段难以收敛的情况，利用重启动方法结合显示动力学计算了两种缺陷模式下屈曲之后的失稳形态和应力分布随载荷的变化情况。需要注意的是，由于球膜的大变形，球膜之间会发生接触，需对接触进行设置。

对于单个临界缺陷模态的初始构型，发现气球最终稳定为4上4下的整体屈曲构型。图9显示了气球从初始缺陷构型充气到500 Pa的过程中位移变化情况。在充气刚开始时(100 Pa之前)，气球仅发生轻微变形，位移变化较小，位移云图开始出现不均匀性；继续充气，压力从100 Pa变化到230 Pa的过程中气球仍未出现明显变形；从230 Pa开始，变形速率显著增大；在300 Pa之后再次减小。气球变形速率变化的转折点位于230 Pa附近，大约为非线性后屈曲的临界失稳压强。在整个充气过程中，气球的缺陷模式被逐渐放大。

球膜的Mises应力云图如图10所示。高应力区位于气球凸起的波峰，这可能是试验气球在高压下爆裂的区域。低应力区位于波峰之间的接触区域。观察应力图例条可知，应力值出现两极分化，应力分布极不均匀，且最大应力值(近似取22.1 MPa)远高于未失稳状态时的最大应力值(约 4 MPa，可参考文献[27])。南瓜球在承压能力、应力分布方面的优势，在面对失稳问题时消失。

图10 Mises应力云图Fig.10 Mises stress contour

对于组合缺陷模态的构型，气球在500 Pa压强下呈现局部变形，如图11所示。与整体变形构型(4上4下)不同，图11显示了一种非均匀变形模式，其中几幅球膜形成局部塌陷，而球体表面的其他部分比较均匀。与整体变形模式中应力分布不同，高应力区位于局部塌陷的区域，并且应力水平高于整体变形构型中的应力水平。

中外南瓜球的地面试验曾多次出现以上两种失稳构型。一种是整体变形模式，气球具有循环对称的变形形状，如4上4下；另一种是局部变形模式，美国国家航空航天局(National Aeronautics and Space Administration，NASA)在气球试验中多次发生S形裂缝的失稳形态，中国科学院在地面试验中曾出现过局部塌陷(图12)的失稳形态。两种不同的失稳形态可能是由于初始缺陷模式不同造成的。

图11 Mises应力云图Fig.11 Mises stress contour

图12 地面试验中南瓜球出现局部失稳Fig.12 Local instability of pumpkin balloon in ground test

在整体屈曲和局部屈曲构型中，加强筋出现严重变形。图13显示了裂缝区域周围严重变形的加强筋形状。其应力云图显示，无论是整体屈曲还是局部屈曲，在球膜接触区域内，加强筋应力要相对偏低，并且其他区域内的加强筋具有几乎均匀的应力值。这可以解释为，在球膜的接触区域，加强筋均位于较低位置，处于相对“松弛”状态，这将导致应力降低。