ANSYS经典案例在Workbench中实现 | 薄壁结构的屈曲与后屈曲分析

在前四期的ANSYS经典案例在Workbench中实现的分享中，我们分享了ANSYS经典案例在Workbench中实现之汽车刹车盘制动噪音分析、密封圈仿真分析以及基于网格重划分的金属成型仿真分析和涡轮机叶片冷却过程的热应力分析（可以点击蓝色文字查看前期内容），本期为大家分享薄壁结构的屈曲与后屈曲分析。

薄壁结构的屈曲与后屈曲分析

案例背景

屈曲分析对于一个成功的结构设计，尤其是包含壳和梁的结构，是至关重要的。虽然线性特征值屈曲分析相对直接与简便，但是也有其自身缺点：因为实际屈曲过程是一个非线性（大变形）过程，如果不能考虑结构非线性，分析只能得到近似结果，另外线性屈曲分析对于结构后屈曲分析无能为力。非线性屈曲分析过程较为复杂，同时可能需要多次尝试才能得到较为可信的结果，但是由于其不存在线性屈曲分析的局限性，所以工程上倾向通过非线性屈曲来评价结构的稳定性。

实际中，工程师很难判断结构究竟何时开始发生屈曲。从工程和科研角度看，人们在整个屈曲过程中，最感兴趣的阶段其实是结构将要产生大变形，但是尚未产生较大变形的阶段，有时结构甚至还未产生变形，因为此时对应的载荷是结构的临界屈曲载荷。非线性屈曲分析可以很好得在这方面提供工程意义上的指导。非线性屈曲分析通过使用以下一些方法，控制整个仿真计算的收敛性，达到用户的工程需求：

1 非线性稳定性控制（nonlinearstabilization）

该方法可以应对屈曲分析中的局部和整体不稳定性，并且可以与其它非线性控制技术联合使用进行仿真（弧长法除外）；

2 弧长法

该方法只能处理力载荷下的结构整体失稳。

3 将稳态分析处理成“准静态”的动力学问题

该方法通过使用动力学效应防止计算发散，但是具体操作较为复杂。

本案例通过承受外部静水压力载荷的周向加强筋圆柱薄壁结构，说明如何通过仿真分析，预测结构的屈曲载荷和后屈曲状态，同时介绍控制非线性屈曲分析中，控制计算收敛性的方法。

问题描述

圆柱薄壁的材料为2024-T3铝合金，由五层横截面为Z型的周向加强筋支撑，圆柱薄壁两端由两个厚盖板（厚度为25mm）密封，并分别由一个L型的铆接条加固。圆柱薄壁承受外部均匀压强，从而使圆柱薄壁上两个Z型加强筋之间的局部屈曲，最终导致结构失效。

具体的几何尺寸如下：

	尺寸（mm）
圆柱薄壁截面半径	355.69
圆柱薄壁深度	431.8
圆柱薄壁厚度	1.034
盖板半径	380
盖板厚度	25
Z型加强筋厚度	0.843
L型铆接条厚度	1.64

Z型加强筋横截面尺寸如下图所示：

图1 Z型加强筋横截面形状及尺寸

L型铆接条横截面尺寸为19*19mm，厚度为1.64mm。

最终的几何模型截面如下图所示：

图2 几何模型剖面图

材料参数输入如下表所示：

Young's Modulus：73GPa

Poisson's Ratio：0.33

Yield Stress：268.9MPA

Tangent Modulus：73MPa （该值为估计值，参考文献中没有明确告知具体数值）

边界条件与载荷

该模型可以看作是潜水艇或航天飞行器的舱体结构，所以可以近似认为结构是悬浮的。本案例中，唯一的约束是控制顶盖板三个位置处的六个自由度，以防止整体结构所有的刚体平动和转动。载荷为外表面均布压强载荷0.24MPa。

仿真步骤

1 特征值线性屈曲分析

线性屈曲分析可以预测理想线弹性结构的理论屈曲载荷，同时还可以提供用于后续非线性屈曲分析中初始缺陷定义的屈曲模态振型。因此，线性屈曲分析对于整个分析很有必要。在ANSYS中，通过特征值屈曲分析进行结构的线性屈曲评估。具体步骤是：

（1）在静力学中定义所需的边界条件和载荷（通常为单位载荷，有时视具体情况而定）；

（2）进行特征值线性屈曲分析，输出理论屈曲载荷和屈曲振型。

图3 特征值线性屈曲分析流程

图4 约束顶盖板三个位置处的平动自由度

图5 施加压强载荷0.24MPA

定义输出前十阶屈曲模态，点击solve进行计算。

2 非线性屈曲分析

在进行非线性屈曲分析之前，有一点必须明确：如果结构是完全理想的对称结构（事实上没有结构是完全对称的），理论上在仿真中是不会出现非线性屈曲失效，正是由于各种公差和工艺上不可避免的缺陷，才产生屈曲失效的可能。所以，在非线性屈曲分析中，对于本案例这个完全理想的对称结构模型，需要引入结构的初始缺陷。

通常有两种方法对结构定义初始缺陷：

（1）在结构上定义一个很小的扰动载荷。在此并不推荐使用这种方法，因为很难把握应该施加多大的扰动载荷，同时也无法判断在结构的什么位置施加这个扰动载荷。如果扰动载荷过大，那么将会得到与实际完全不同的结果。

（2）通过屈曲模态振型配合缩放系数定义初始缺陷。通过线性特征值屈曲分析，可以得到结构的屈曲模态振型，这些振型可以通过命令流引入，直接对单元节点的坐标位置进行更新。建议引入多阶屈曲振型模态，同时应保证初始缺陷的量级与制造公差的量级接近（通常小于1%）。

对于本案例而言，将线性特征值屈曲分析的前十阶屈曲振型的叠加作为结构的初始缺陷，同时考虑到圆柱薄壁的半径为355.69mm，制造公差的量级约为0.1mm，而屈曲振型的量级在1mm，所以选择缩放系数为0.1。

在Mechanical中插入一下命令：

*do,i,1,10

upgeom,0.1,1,i,buckling,rst

*enddo

在上述命令中，buckling为线性特征值屈曲分析的结果文件，注意该名称应包含文件的路径，此处已省略。相关UPGEOM的命令介绍，可参见ANSYS的Help文件。

打开结构大变形，定义均布压强载荷为0.24MPa，同时使用尽可能小的时间步长，以相对准确地估算屈曲临界载荷（这里最大子步数取500）。在计算过程中，尤其是在结构达到屈曲临界点或后屈曲阶段，计算往往很难收敛，所以推荐使用重启动计算。一旦出现不收敛的情况，可以及时调整设置，在指定时间点重新启动计算。

3 后屈曲分析

非线性静力学计算的发散点通常意味着结构的屈曲失效即将或已经出现，此时结构的变形并不明显，所以屈曲的初始点很难被肉眼所观察到。但是，通过仔细评估仿真的载荷-位移曲线，可以判断屈曲何时发生。

通常采用三种方法判断结构是否发生屈曲：

（1）迭代发生发散导致计算失败；

（2）最大位移值在一个很小的时间增量内发生突变；

（3）最大位移值在符号上发生变化。

图6 非线性屈曲计算的力收敛曲线

在本例中，非线性屈曲计算在第十四个子步，约0.8s左右（不同的工作站以及不同的软件版本可能会导致上述时间点发生偏差）发生不收敛，计算失败，这是由于结构屈曲导致的。接着进行重启动计算，以进行结构后屈曲分析。由于结构的后屈曲过程极其不稳定，需要进行特殊的设置以控制计算的收敛。在静力学分析中，尤其是在结构局部屈曲分析和时间相关非线性材料分析中，非线性稳定性控制（nonlinear stabilization）是一个非常好的选择。

此外，在进行重启动计算时，尽量避免在计算发散的时间点重启动计算。比如本例中，计算在第十四个子步发散，在重启动计算时，选择在第十个子步设置nonlinear stabilization并启动计算。如果计算依旧不收敛，则应选择更靠前的时间点。

在后屈曲分析中，需要设置非线性稳定性控制（nonlinear stabilization）以保证计算收敛。通常有两种方式实现：定义阻尼系数或定义能量比。用户需要选择通过哪种方式进行控制，然后输入适当的值。理论上，用户应该根据结构失稳类型、网格尺寸与类型，以及时间步长的综合考虑进行设置。用户可能需要经过多次尝试才能寻找到相对合理的输入。总之，nonlinear stabilization的目的是以最小的附加阻尼力获得整个计算的收敛性。

（1）阻尼方法：如果用户使用阻尼方法，意味着阻尼将定义在所有单元上。一旦阻尼值过大，附加在结构上的控制也会很大，因此结构刚度相应偏大，即结构很难产生变形；如果阻尼值过小，那么该设置不会起到控制计算收敛的作用。综上所述，阻尼控制方法不适合用于局部非线性屈曲失效的结构。

（2）能量方法：由于本案例中的薄壁结构的失效是由于局部屈曲失效引起的，所以能量方法更加适合本案例。能量方法的原理是在不同的单元处施加不同的阻尼，所以程序会给容易失稳的单元施加较大的阻尼，而给不易失稳的单元施加较小的阻尼，这样既可以控制计算的稳定性，同时也可以将附加控制力降低到最小。本案例中取能量比为0.000143进行计算。

计算结果

下图所示为线性特征值屈曲分析的屈曲阵形图：

第一阶：0.660

 第二阶：0.660

第三阶：0.664  

 第四阶：0.664

第五阶：0.666

第六阶：0.666

第七阶：0.679

第八阶：0.679

第九阶：0.689

第十阶：0.689

图7 线性特征值屈曲分析振型图

所有的特征值屈曲分析的临界载荷因子（即特征值）大约为0.66，线性特征值屈曲分析的静力学部分所施加的载荷为0.24MPa（见图5），所以理论临界屈曲载荷约为0.66*0.24=0.158MPa。

下图为底盖板中间某节点的轴向位移随外载荷变化的曲线结果。不难发现，在外载荷达到0.15MPa左右，曲线的曲率发生很大的突变，这表明结构发生了屈曲失效，对应的载荷即为结构的屈曲载荷。

图8 底盖板中间节点轴向位移随外载荷变化曲线

下图为结构整体最大变形位置处的径向位移随外载荷变化曲线，通过曲线可以看出，在结构发生屈曲失效时，曲线的斜率发生了方向性的改变，这也很好地体现了结构屈曲失效正在发生。当外载荷达到约0.2MPa时，位移发生瞬间突变，从该点往后，外载荷虽然增量很小，但是位移增量非常大，最后直至结构完全失效。

图9 最大整体变形位置处节点径向位移随外载荷变化曲线

下图所示为外载荷达到屈曲载荷（约为0.154MPa）时，以及最终状态下，结构的整体变形云图。在结构屈曲初始状态下，肉眼很难观测到结构屈曲的发生，但是，实际上结构屈曲已经产生。

结构初始屈曲变形

 结构最终状态下变形云图

图10 结构整体变形云图

注：用户需要校核收敛控制产生的额外能量（stabilizationenergy），并与结构应变能进行对比，前者应远小于后者，计算才相对可靠。

另外，本案例仅对结构非线性屈曲的仿真方法进行描述。虽然本案例中的模型取材于某文献中的参考模型，但是，由于并非所有输入都与文献中输入一致，有的甚至缺少相应输入，所以结果与文献中有所偏差。

延伸阅读：

1、ANSYS经典案例在Workbench中实现 | 汽车刹车盘制动噪音分析