微曲率壳板横向失稳临界载荷计算方法

付 晓，梅志远，陈国涛，张 二

（海军工程大学舰船与海洋学院，武汉 430033）

0 引 言

舰船的几何外形特征较为复杂，微曲率壳板是其中一类常见的结构。舰船在服役期内会承受各类载荷作用，横向载荷极为常见（如静水压力、波浪砰击等）。尤其在遭遇各类恶劣天气时，船体壳板出现失稳甚至破坏将威胁舰船的航行安全。近年来，复合材料在船舶结构中得到了广泛使用，由于复合材料与金属材料力学性能存在明显差异，复合材料壳板的横向承载稳定性也逐渐受到关注[1-2]。虽然舰船因船体壳板失稳导致的事故鲜见报道，但是海军工程大学梅志远团队在针对大尺寸微曲率复合材料板架模型开展力学性能试验时却发现了壳板失稳的现象。由此可见，微曲率壳板的横向承载稳定性机理还需深入研究。

由于壳板轴向力的存在，在横向载荷作用下，微曲率壳板普遍会出现某一时刻位移突增但仍可继续承载的现象，有学者将其命名为“跳跃”失稳[3]。多数学者将此现象归结为壳板的屈曲与后屈曲问题[4-10]，并对此进行了理论分析与求解[11-15]。在20世纪70年代，中科院力学所十二室[5-9]就对加筋曲板的侧压稳定性问题开展了较为系统的理论研究，形成了比较完整的分析思路。但由于当时技术条件的限制，相关研究缺乏仿真与试验的支撑。计算机的发展以及有限元技术的普及，为深入研究“跳跃”失稳现象的产生机理与影响因素提供了条件。Budiansky和Roth[16]对浅球壳受突加外压作用的“失稳跳跃”行为进行了研究，将系统动力稳定性准则与浅球壳的动力学控制方程结合起来,描述了结构的“跳跃”特性。陈伟等[17]基于ABAQUS有限元软件研究了双曲率壳板侧压稳定性屈曲与后屈曲载荷曲线，并对影响因素进行了分析。胡文飞等[18]探讨了扁球壳“跳跃”现象的影响因素，发现壳体厚度和矢高之比是主要因素。虽然借助有限元软件可以直观地了解“跳跃”失稳现象的变化过程与影响因素[19-28]，但是由于此问题理论分析过程较为繁琐，复杂的计算也使解析解的求解难度骤增，导致数值分析的准确性难以得到验证。现阶段，针对微曲率结构失稳问题可供借鉴的资料较为匮乏，相关文献一般多讨论壳板或加筋板受轴向压力时的失稳问题，对于横向载荷作用下壳板稳定性问题的研究也多停留在“跳跃”失稳现象仿真分析层面，缺乏理论支撑。因此，建立可适用于工程实际的微曲率壳板理论分析模型，探讨微曲率壳板的“跳跃”失稳机理，有利于掌握其横向承载规律，避免结构安全问题的发生。

本文根据铁木辛柯曲梁理论，推导了单向微曲率壳板的临界失稳载荷公式，并基于Riks弧长法对典型模型进行数值计算，获取其临界失稳载荷，以期为微曲率壳板的横向承载稳定性预报提供参考。

1 微曲率壳板“跳跃”失稳现象

微曲率板架广泛存在于舰船结构中，舰艇航行时，其主要承受静水压力、流击载荷以及波浪砰击载荷作用，需要具备较高的强度与足够的稳定性。本文在对船体微曲率板架模型进行力学性能试验时，发现壳板在承受外压时出现了大面积凹曲的现象。然而，出现大面积凹曲后，结构并未破坏，仍能继续承载，即认为壳板发生了“跳跃”失稳。试验所用玻璃钢壳板模型如图1所示，试验所得“跳跃”失稳现象如图2所示。

2 微曲率壳板横向承载稳定性理论模型

对于薄板而言，当板的挠度与板厚为同一量级时，在推导板的微分方程时就必须考虑附加在中面内的薄膜应力，几何关系不再满足线性假设，板的小挠度弯曲理论以及克希霍夫—勒夫假设均不再成立[29]。微曲率壳板的理论求解方法过程复杂，应用于指导工程设计存在一定困难。参考曲拱结构中的矢跨比概念[30]，本文定义壳板最高点与最低点竖直方向距离与跨距的比值为矢跨比，矢跨比小于0.05 的壳板可视为微曲率壳板。在几何特征方面，虽然该结构拥有一定的曲率，但是曲率较小，与平板较为相似；在变形特征方面，由于微曲率壳板存在的轴向力，其承受横向载荷时，容易出现“跳跃”失稳等问题，又与曲板变形相似。由船舶结构力学可知，求解薄板柱面弯曲问题时，常采用板条梁的方法，将二维板模型退化为一维梁模型进行分析。针对本文研究的单向微曲率壳板结构，选取单位宽度板条梁模型（板厚为t），进而分析结构在承受横向均布载荷下的稳定性问题。

以铁木辛柯为代表的部分学者完成了曲线梁变形的理论推导，认为可以使用承受轴向力以及横向载荷联合作用时的直梁模型求解微曲率曲梁（杆）横向承载问题[31]，其示意图如图3所示。

对于上述简支梁模型，根据文献[31]求解，首先应判断其是否失稳，即比较轴向力T与结构欧拉临界失稳载荷Se的大小关系。临界失稳载荷TE可表示为式中：E表示材料弹性模量；I表示结构抗弯惯性矩；n表示结构的失稳半波数，结合工程实际与试验数据，发现在本文研究背景下的微曲率壳板失稳变形特征可近似为3个半波，因此n暂时取3。

根据文献[31]可计算图3所示梁的最大挠度为

由式（1）～（3）可知：当梁所承受的轴向力T与欧拉临界失稳载荷Se相比很小时，γ的值很小，式（2）中第二个因子数值接近于1，说明在此时轴向力对挠度的影响可以忽略；当轴向力T与欧拉临界载荷接近时，γ值接近π/2，则式（2）中第二个因子无限增大，导致梁的挠度急剧增加。

轴向力对结构失稳至关重要，在研究此问题过程中，在铁木辛柯微曲梁推导的基础上，假设结构的轴向力达到欧拉临界载荷时结构发生失稳。因此，以承受轴向压力的薄板板条梁模型为对象，通过求解该板条梁的大挠度复杂弯曲问题，反推当梁轴向力达到临界失稳载荷时结构所承受的均布载荷q，该载荷即为结构临界失稳载荷qcr，其推导过程如下：

针对微曲梁，首先求解其挠曲线方程，微曲梁的挠曲线方程通过级数法表示，微曲梁初始长度如图4所示。

本文仅取式（4）中第1项表示微曲梁的挠曲线，已知微曲梁的曲率半径R、跨距l，即可由下式求得系数a1。

微曲梁初始状态下梁的长度为

则其初始长度较直梁的伸长量为

然后以直梁模型代替微曲梁模型，进行理论推导。取图5所示板条梁微段dx，其变形后长度为ds，则有

根据结构力学板的柱面大挠度弯曲理论，将式（8）中w'展开为幂级数后，可求得整个板条梁变形后的伸长量为

则微曲梁的整体伸长量为

继而由应变的定义可得

同时，板的柱面弯曲求解过程中参数u与轴向力T之间存在如下关系式：

式中，D表示板的弯曲刚度，根据结构力学中复杂弯曲梁的微分方程，可推导出轴向压力作用下的简支直梁复杂弯曲挠曲线方程为

考虑到微曲梁初始曲率的影响，则简支状态下微曲梁复杂弯曲的挠曲线方程可表示为

则有

又令w0(x)=q·g(x)，则可得

根据假设，已知失稳时，轴向力达到失稳临界载荷，因此，T=TE，将式（11）～（16）联立，即可得到失稳时均布载荷q的表达式，再由式（12）解得u，最后临界失稳载荷qcr可由式（17）计算。

对于ICP-AES分析，可调节的仪器参数主要有射频发生功率、工作气体流量(包括冷却气、辅助气、雾化气)、蠕动泵转速、观测方式等。其中射频发生功率、雾化气流量和辅助气流量是影响分析线信号的关键因素[19-20]。

式（17）用于计算单向微曲率壳板的临界失稳载荷，可对微曲率壳板的极限承载能力进行预测，进而对结构安全进行评估。

3 基于Riks弧长法的微曲率壳板横向承载稳定性数值分析

3.1 改进弧长法

曲壳在凸壳面承受分布压载时，初始阶段容易在曲壳顶部出现反向凹曲的变形特征，此后，随着横向压载荷的增加，凹曲边界将持续稳定地扩展[32]。研究此问题必须考虑几何大变形的影响，平衡方程和几何关系均存在非线性关系，应变的表达式中包括位移的二次项[33]。因此，对于均布载荷下微曲率壳板的变形问题，本文采用非线性方法对其求解。

目前，一般使用弧长法[34-42]对曲壳结构进行非线性分 析，该 方 法 由Wempner 和Risk 提 出，后 经Roma 和Crisfield 等人改进，形成了改进的弧长法。改进弧长法是一种稳定高效的结构非线性分析方法，对于结构的非线性前屈曲以及后屈曲的路径跟踪较为有效，改进弧长法是通过约束方程进行迭代求解收敛点来进行非线性静力学求解。图6为改进弧长法的迭代过程[43]。

3.2 微曲率壳板数值分析模型

本节基于Riks 弧长法，选取典型微曲率壳板对其开展数值仿真，以验证理论计算模型的准确性。为不失一般性，结合船舶结构中常见微曲率壳板实际情况，建立单曲率矩形壳板模型。对各类壳板开展数值计算，其中曲率半径分别为4000 mm、6000 mm、8000 mm、10 000 mm，跨距为1000 mm，长宽比分别为1、1.5、2、3，探讨不同曲率与长宽比矩形曲壳板的横向承载稳定性规律，图7所示为长宽比1.5，曲率半径为6000 mm时的壳板承载示意图。

数值模型采用实体单元建模，材料为碳纤维（T700/350），参数为：E1=58.7 GPa，G12=3.32 GPa，ρ=1.46 g·cm-3，υ=0.045；铺层方式为0°/90°正交铺层，壳板的厚度为16 mm，表面施加200 kPa 均布载荷，边界条件为简支约束。

模型的网格划分如图8 所示，采用四边形SC8R 单元，网格密度为20 mm，壳板厚度方向采用数值积分法来求解刚度矩阵，板壳的面内采用减缩积分。模型采用Riks弧长法进行求解，可防止在载荷施加时壳板屈曲过程中出现短时间内变形过大的现象，导致计算不收敛。

在非线性分析过程中，可通过绘制壳板的载荷比例因子（LPF）曲线，并通过第一次出现的峰值点位置计算壳板临界失稳载荷，不同曲率半径壳板的LPF曲线如图9所示。

图10 是跨距为1000 mm，长宽比分别为1、1.5、2、3的平板LPF 曲线。由图可知，长宽比越大，单位时间内承受载荷越大，曲线斜率越大。与图9所示不同曲率壳板的LPF 曲线对比可知，曲率对“跳跃”失稳现象影响较明显，长宽比一定时，曲率半径越小，“跳跃”失稳现象越容易发生，对应的临界失稳载荷越大。然而，在曲率半径一定时，长宽比对“跳跃”失稳现象的影响较小。但是，矩形壳板与正方形壳板的LPF 曲线仍存在一定的差异。此外，结合图9 和图10 可知，壳板长宽比越小，“跳跃”失稳现象越不容易出现，同一曲率半径不同长宽比壳板对应的失稳载荷差别不大。

3.3 数值与理论计算结果对比分析

根据数值计算结果可知，微曲率壳板的“跳跃”失稳现象与壳板曲率密切相关，在本文研究背景下，在壳板曲率半径超过6000 mm后该现象较难出现。为验证理论计算模型，本文选取曲率半径分别为4000 mm、4500 mm、5000 mm、6000 mm 的不同长宽比微曲率壳板，将3.2 节数值计算中模型的相关要素代入理论计算模型中，通过式（17）计算微曲率壳板的临界失稳载荷，与表2中数值计算所得临界失稳载荷进行比较，并对理论解与数值解存在的差异进行误差分析，数值计算结果与理论计算结果如表2所示。

由表2可知，数值解与理论解存在误差，且随着曲率半径的增加，误差会小幅上升，这是由于曲率半径增加后，壳板的变形特征向平板趋近，因此导致误差出现一定程度累积。随着矩形微曲率壳板长宽比的增加，临界失稳载荷值会小幅降低，正方形壳板较难出现“跳跃”失稳现象，半径越大，壳板越难出现该现象。综上所述，可以认为本文提出的壳板临界失稳载荷计算方法针对曲率半径在4000～6000 mm范围内的单向微曲率壳板具有较高的准确性。

4 结 论

本文基于铁木辛柯曲梁理论与薄板大挠度弯曲理论，提出了单向微曲率壳板的横向承载稳定性理论计算模型，对临界失稳载荷进行了预报，通过理论与仿真对比分析，得到以下结论：

（1）微曲率壳板在横向承载时会出现稳定性问题，本文通过理论分析与数值仿真计算，发现单向微曲率壳板的“跳跃”失稳现象主要与曲率半径相关。半径越小，“跳跃”失稳现象的临界失稳载荷越高；随着壳板曲率半径的增加，壳板承载特征越接近平板，“跳跃”失稳现象越不容易出现；在壳板曲率半径相同的情况下，矩形壳板的长宽比对临界失稳载荷影响不明显。

（2）通过单曲率柱面大挠度弯曲薄板的板条梁模型可以求解微曲率壳板横向承载问题，假设以轴向力接近欧拉失稳临界载荷作为判断“跳跃”失稳的发生条件，求解过程中需考虑轴向力作用以及微曲率壳板初始曲率的影响。

（3）由壳板临界失稳载荷的理论计算结果与数值计算结果对比可知，本文提出的临界失稳载荷计算方法具有较高的准确性，能够较好地预测微曲率壳板的横向承载稳定性。