ABAQUS非线性屈曲分析实战案例 原创

 屈曲分析通常指通过求解特征值问题来确定结构失稳的临界载荷。在已知材料弹性模量和泊松比的线弹性前提下，结合结构的几何形状与边界条件，可计算得到屈曲特征值，进而得出临界载荷，为后续的材料强度评估或非线性屈曲分析提供参考依据。后屈曲分析则进一步研究结构失稳后的力学行为，通常采用非线性分析方法考察结构的屈曲后承载能力，从而更全面地评估结构在实际工况下的安全性与稳定性。

1引言

如果我们在abaqus中建立一个细长杆，赋予基本的材料参数（弹性模量、泊松比），将杆的一段固定，另外一段施加载荷（力/位移），可以发现杆件仅仅会沿着载荷施加方向无限缩小，并不会出现预想的杆件弯曲，如下图所示。这显示是不符合实际情况的，幸运的是abaqus中提供了针对屈曲问题的分析方法，可以实现对受压物体屈曲结果的预测，本文使用的是2022中文版的abaqus。

2 理论部分

abaqus中的屈曲分析又可分成线性分析和非线性分析，线性分析是非线性分析（后屈曲分析的基础）。

2.1 线性屈曲问题

线性屈曲，可以理解成简单的特征值求解问题。在仅仅知道弹性模型和泊松比的情况下，求解出模型的最大承受的载荷。屈曲问题的特征值通用求解方程可以表示为：

其中，表示模型在承受外部载荷下的切线刚度矩阵，这些载荷包括力、位移、热加载等，即包含了模型当前的边界条件和模型本身的刚度矩阵。表示位移解。

在增量求解模式中，特征值求解方程又可以进一步表示为：

 表示模型的基本结构刚度矩阵，包含上一个分析步的边界条件，但不含当前分析步的边界条件，通常在线性分析中。表示由于边界条件引起的初始应力和加载刚度矩阵。 表示第阶屈曲模态下的特征值（这是在云图结果中可以直接看到的），表示第阶屈曲模态下的特征向量。M、N则分别表示模型的自由度，如果是实体模型，M、N取值范围为1-3，如果是平面问题的求解，M、N取值范围为1-2。由于是线性问题，因此所有的刚度矩阵都是正定的，如果是非正定的，则特征值的求解较为困难。

因此，在线性求解中，模型承受的最大载荷，可以表示为：

为由于位移、力、热等引起的边界条件

2.2 非线性问题

非线性屈曲问题又可成为后屈曲问题，是在线性屈曲的基础上，通过在后屈曲分析步中引入微小缺陷，使模型呈现几何非线性，进而求解模型变形和承载能力的问题。在后屈曲分析中，材料可以为线弹性、弹塑性、弹塑性等，材料性质的选择可以更加接近于实际，因此求解出的承受载荷也更加接近于实际。由于非线性求解问题较为复杂，这里不再以公式的形式描述问题的求解过程，大家感兴趣的，可以去abaqus 的官网细看 。

3 案例实战

这里以中科院二区（Thin-Walled Structures）论文 《Experimental investigation on the stability of aluminium alloy 6082 circular tubes in axial compression》针对铝合金圆形管轴向压缩稳定性研究为基础，详细介绍仿真流程，并将仿真结果与计算结果进行对比，验证模拟的合理性。

3.1 O89-L500试件的非线性屈曲分析流程

3.1.1 建模过程

论文中，每种不同尺寸的圆形钢管有三个样品，我们这里以三个样品实际测量尺寸的平均值作为建模尺寸。模型尺寸如下图所示：

​  

abaqus建模尺寸图如下所示，模型为三维实体建模，采用mm单位制。

  为便于后续模型修改，建议对模型尺寸进行标注。

3.1.2 弹塑性本构模型

材料本构模型采用Ramberg–Osgood弹塑性变形模型，表达式为：

这是论文中的公式1的形式，、、、分别表示为弹性模型、应变、应力、塑性硬化指数，表示弹性应变为0.2%时的对应应力。

 abaqus官方文档中给出的公式形式为：

由于论文中测得的应力值为正，已知，并对公式进行移项变换可得

因此，通过对比可以得出塑性缩放因子的值。

我们取论文中测量的平均值作为材料参数，如下表所示    

abaqus设置材料参数的界面和流程如下所示：

3.1.3 分析步设置

我们首先进行线性屈曲分析，选择分析步-创建分析步-线性摄动-屈曲，进行设置，具体设置界面如下所示。

​

注意的是，这里我们只对模型进行了第一阶模态的求解，这里采用子空间求解，最大迭代次数默认为30（这个最大可设置为300）。在求解的特征值数量<10时，采用子空间求解速度较快，在特征值数量>10时，采用Lanczos求解速度较快。

3.1.4 边界条件

首先在圆筒模型的上下两个分别创建两个参考点，并将这两个参考点分别与上下两个底面进行耦合约束，具体如下所示：

        ​

由于实验中的边界条件约束为两端铰接，因此参考点RP2的边界条件设置为U1=U2=U3=UR3=0,参考点RP1的边界条件设置为U1=U2=UR3=0,UR3暂定设置为-50，这样设置是为了后屈曲分析更加便利。实际上，在线性分析阶段，这个值的大小和最终求的特征值的乘积是固定的，即线性最大承载载荷，而对后屈曲分析，UR3对引入缺陷并无影响。

3.1.6 网格划分

采用六面体网格划分，通过预先的网格敏感性分析，结合计算量，设置网格尺寸为3mm，采用中性轴算法对模型进行扫掠划分。网格划分流程和划分结果如下所示。

3.1.5 关键字编辑

在提交分析计算前，需要进行设置参考点模态位移的输出，以便于后续引入缺陷，设置过程如下所示：

​

*NODE FILE表示节点文件输出，U表示输出的变量为节点位移。

3.1.6 线性屈曲结果

所有设置均完成后，进行作业提交，这里作业命名为Q89-500-GG，线性屈曲模态结果如下所示：

所求的特征值的大小为0.29579，这里施加的是位移载荷，因此，线性屈曲分析得到的模型的最大承载能力为0.59579×50=14.7895mm。这里模态显示的位移是归一化之后的结果，最大为1。

3.1.7 后屈曲分析

首先将Q89-500-GG模型进行复制，重命名为Q89-500-F-GG。这里需要注意，不要重新对模型进行网格划分，因为在引入微小缺陷时，实际引入的是节点位移，这要求后屈曲分析的模型的节点编号要和线性屈曲分析阶段的节点编号相同。

3.1.7.1 关键字编辑

在Q89-500-F-GG模型中进行编辑关键字时，首先把上一步添加的关键字删除（*NODE FILE U）,提高计算效率。新的关键字的添加位置为材料参数定义的下一行，具体如下所示。

​

其中*imperfection表示引入缺陷，file=Q89-500-GG表示引入节点缺陷的文件，即线性屈曲模态的节点计算结果，step=1表示从作业Q89-500-GG文件中的第1个分析步读取数据，“1，0.4”表示提取第一节模态节点位移0.4倍作为初始缺陷引入几何模型。

这里需要注意的是在非线性屈曲中，通过将线性屈曲的节点位移作为初始缺陷引入时，如果引入的初始缺陷较大，就相当于真实模型具有较多缺陷，与实际不符，致使得到的最大承载力偏小、结果失真，如果缺陷设置的较小，又很难引起屈服。这里初始缺陷设置为模型长度的0.1%。以O89-500一阶模态分析结果为例，模型的最大位移为1.231，长度为500mm，故引入的初始缺陷应该为0.5，即一阶模态节点位移的0.4左右。

需要说明一点，abaqus的官方文档中阐述模态分析的位移是归一化后的位移结果，由于论文中圆柱管的最大长度为1650mm，因此我在对缺陷节点位移的缩放比例分别为0.01，0.1，1进行来了后屈曲分析，结果显示，当缺陷比例在1时，模型产生弯曲屈曲，在0.01和0.1时，模型并没有产生屈曲，这说明，引入的缺陷过小，模型过于完美，大家也可以试着调整引入的缺陷大小，总的来说，只要引入的缺陷大小在1左右，结果不会产生过大的偏差。

 3.1.7.2 修改分析步

 在后屈曲分析时，需要重新定义分析步的类型，这里我尝试分别采用静态通用和静态-Riks两种不同的分析步类型进行分析。通过提取参考点RP1的力和位移，得到分析结果如下图所示：

​

从图中可以看出两种分析步得出的力-位移曲线几乎重合，因此为了方便得到更加光滑的曲线，这里采用静态通用进行分析，分析步设置界面如下。

几何非线性设置为打开，为了得到更加光滑的曲线，初始增量步和最大增量步设置为0.002，这里大家如果为了提高计算效率，可以适当设的大一些，但是得到的力-位移曲线就会粗糙一些，尽量做好平衡。为提高计算的收敛性，最大增量步数设置为1000，最小增量步数设置为1E-20。

3.1.7.3 提交计算

到此，所有的后屈曲分析的模型设置都已经完成，直接提交分析计算即可。

3.2 O89-L500非线性屈曲分析结果

提取参考点RP1的力-位移曲线，得出最大载荷为444.76KN，与式样实验得到的平均值482.71KN之间的误差为8%，小于10%以内，力-位移曲线如下图所示：

​

最终的屈曲失效的应力云图如下所示，发现O89-L500的屈曲失效形状皆为中间弯曲的形态，弯曲部分的塑性行为最先出现且最明显，这符合实际，进一步说明了仿真的可靠性。

3.3 其它试件分析结果

对论文中的其它试样进行仿真时，在建模时的尺寸均选取为试样的平均测量尺寸，具体尺寸参数如下所示。

​

所有试样均采用Ramberg–Osgood弹塑性变形模型，材料参数均相同。所有试样的分析步骤均和O89-L500的分析步骤相同。对其它试样进行分析，最终得到参考点RP1的力-位移曲线如下所示。

将仿真得到的峰值力与文献进行对比，并计算误差，如下所示：

预测值与实验值的平均误差在10.94%左右，多数试样在10%以内，说明了仿真的可靠性。为进一步分析屈曲形态，我们将不同试样的应力和等效塑性应变云图进行输出，如下所示。

​

从图中可以看出不同试样的塑性应变总是最先出现在模型的中部，且随着试样长度的增加，两端的应力水平越小。通过对比文献中仅有的三种试样的最终屈曲形态，可以看出，仿真结果与实验结果吻合较好，都表现出中部弯曲的状态。

4 总结

总的来说，通过abaqus进行后屈曲分析得出的预测结果和实验结果较为吻合，说明该仿真方法的可行性。为避免出现问题，在仿真中，大家尤其需要注意一下几点：

1. 后屈曲分析的模型直接来源于线性屈曲的复制文件，不要另外重新划网格，保持节点标签的一致性

2. 在进行参考点位移约束时，一定要根据实际情况进行约束，如两端完全固定约束和两端完全铰接约束的计算结果相差很大，因为边界条件会直接影响屈曲模态的结果。

3. 在引入微小缺陷时，可以通过逐渐尝试改变引入缺陷缩放因子的大小，来确定最优的值。可以将刚刚能够引起模型屈曲的效果作为评判标准，但缩放因子的值不宜太大，建议不要超过10。

4. 在进行后屈曲分析步的设置时，如果使用静态-Riks分析类型，则要控制弧长增量的大小，以保证计算结果的准确性。在我对该模型进行分析时，静态通用和静态-Riks两种分析结果几乎完全重合，但静态通用的分析时长更短一些。

​ ​