几何非线性-屈曲分析

YORUNOKUNI

已于 2023-08-14 11:04:28 修改

阅读量6.2k

点赞数 17

文章标签：算法

于 2023-04-19 22:32:08 首次发布

本文链接：https://blog.csdn.net/YORU_NO_KUNI/article/details/129946228

版权

1.前言

本文首先对屈曲分析的一些概念和原理进行介绍，之后描述如何使用ANSYS Workbench进行屈曲分析，对软件使用的算法也会有适当介绍，如在上一篇几何非线性文章中，将载荷分步加载，使用牛顿法和一些收敛准则控制求解，但是针对本文的后屈曲问题，牛顿法并不能很好的得到收敛解，为解决该类问题就要使用别的手段，例如非线性稳定性控制或弧长法等。

2.失稳与屈曲

2.1 失稳

2.1.1失稳的定义

结构系统的破坏分为两种，一类是由材料软化而引起的强度破坏，另一类是由几何软化引起的结构较大变形，从而使结构丧失承载能力。本文所要讨论的问题属于第二类，稳定和屈曲具有不同的概念，但结构失稳的原因都可以归结为几何软化。

在材料力学中是这样定义压杆失稳的：中心受压直杆(理想模型)在直线状态下平衡，由稳定平衡转化为不稳定平衡时所受的轴向压力的界限值，称为临界压力。中心受压直杆在临界压力的作用下，其直线形态的平衡开始丧失稳定性，简称为失稳。显然，这个定义只适用于理想的压杆模型，那么关于失稳是否存在一般性定义呢？在钱若军的《结构屈曲分析理论和方法》中，给出了失稳的静力学定义和一般性定义。

基于静力学理论的失稳定义

参考《结构屈曲分析理论和方法》的1.3.2.3

一般的失稳定义

参考《结构屈曲分析理论和方法》的1.3.2.4

2.1.2 失稳的分类

结构的失稳可以分为五类，第一类是结构系统因几何欠约束发生结构的几何可变；第二类是结构系统因外部约束不充分发生刚体位移或瞬变；第三类是内部几何稳定且外部约束充足的结构系统，在外载荷作用下发生弹性或弹塑性屈曲；第四类是结构系统的倾覆；第五类是结构可能发生的松弛。本文主要研究第三类失稳问题。

2.2 屈曲

2.2.1 屈曲的定义

屈曲是平衡失稳的一个具体模式，发生在整体或局部受压的弹性体上，或者说屈曲是结构的失稳形式。经典屈曲的定义可以归结为结构的平衡转移，当稳定的平衡状态受到任意小的外加干扰后失去平衡，结构因而处于不平衡状态，这种从稳定到不稳定的转移，称为结构发生屈曲。

2.2.2 屈曲的分类

我们来讨论一下欧拉的压杆公式，由欧拉解求得临界力是多值的，它与杆件屈曲的波形有关。令n为失稳弯曲的半波数，则临界力为 $n^{2}P_{cr}$ ，当n=1时为最低的临界力， $P_{cr}$ 表示了半波平衡的分枝点。如下图所示，OC为失稳前的平衡路线，CB为失稳后的平衡路线，C点即为分枝点。根据上述性质，可以将屈曲分为分枝型屈曲和极值型屈曲。

分枝型屈曲：如上图所示，C点为失稳点，在此处发生平衡状态的分歧，结构发生失稳。因此，分枝屈曲的特征是在稳定平衡的基本状态附近存在这另一个相邻的平衡状态，而在分枝点处发生从基本平衡分枝到屈曲平衡分枝的转换。该类屈曲分析实质上是解决平衡方程的多值性问题。

极值型屈曲：如下图所示，极值型屈曲没有明显的分枝点，但在变形途径中存在一个最大载荷值，达到最大载荷值后变形迅速增加，载荷反而下降。极值型屈曲中还存在一类情况，如下图中的左侧图所示，OA和CD段是稳定的，AC段不稳定，当载荷增加到A点时，平衡状态发生明显跳跃，突然过渡到另一具有较大位移的平衡位置，这种情况称为跳跃失稳(snap through)。对应工程中的案例有受横向均布压力的球面扁壳，Williams双杆结构等。

分枝型屈曲和极值型屈曲对应的载荷值称为临界载荷，相应的状态成为临界状态。将到达临界状态之前的平衡状态成为前屈曲平衡状态，超过临界状态之后的平衡状态成为后屈曲平衡状态。

一般概念认为，对于弹性体系其屈曲载荷可以作为体系承载能力的依据。但是，对于有些类型的结构，如四边支撑的受压薄板，在发生屈曲后仍可继续加载，体系的承载能力要比屈曲载荷大很多。而对另外一些结构，如轴向受压的圆柱壳，其实际承载能力又远小于理论指出的屈曲载荷。这些现象说明，根据屈曲分析得到的屈曲载荷并不总是与体系的承载能力相联系的。之所以产生这种原因，关键在于体系后屈曲平衡状态并不总是稳定的，所以对后屈曲的研究也显得至关重要。

3.有限元屈曲分析

ANSYS中的屈曲分析可以分两类：线性屈曲分析(特征值屈曲)和非线性屈曲分析。

特征值屈曲可以预测一个理性弹性结构的理论屈曲强度，即图(b)中的分叉点。这种分析得到的结果与材料力学中的欧拉公式计算结果相同，但是由于结构的初始缺陷等一系列因素，理论计算结果通常要大于实际的临界载荷。而且特征值屈曲分析本身是一种线性分析，如果存在非线性的设置，都将会按照线性进行处理

非线性屈曲用于计算结构的最大载荷，同时将结构的初始缺陷，以及塑形行为、接触、大变形等因素都考虑其中，因此非线性屈曲的分析结果更接近真实结构的屈曲极限。

因此，特征值屈曲得到的结果通常作为非线性屈曲分析的前置条件：

特征值屈曲分析得到的理论屈曲强度值可以为非线性屈曲分析要施加的载荷提供参考，非线性屈曲施加的载荷可以略高于(10%至20%)理论屈曲强度值，而且特征值屈曲分析得到的屈曲模态形状可以作为非线性屈曲分析的初始缺陷；
在特征值屈曲分析中，如果有的特征值很接近，说明结构对初始缺陷比较敏感。这时应对结构设置初始缺陷，进行非线性屈曲分析；

3.1 特征值屈曲分析

在ANSYS中的特征值屈曲分析原理如下：

首先求解线弹性屈曲下结构的载荷-位移关系：(在特征值屈曲分析前，一般要做预应力分析)

$\left \{ P_{0} \right \}=\left [ K_{e} \right ]\times \left \{ u_{0} \right \}$ (3.1)

其中， $u_{0}$ 为位移结果， $P_{0}$ 为载荷， $K_{e}$ 为弹性刚度矩阵；

假设特征值屈曲的位移很小，则在任意状态下载荷，位移和应力的增量平衡方程如下：

$\left \{ \Delta P \right \}=\left ( \left [ K_e \right ] +\left [ K_\sigma (\sigma ) \right ]\right )\times \left \{ \Delta u \right \}$ (3.2)

$K_{\sigma }\left ( \sigma \right )$ 为某应力状态下对应的初始应力矩阵；

假设特征值屈曲行为是一个线性函数，则：

$K_{\sigma }\left ( \sigma \right )=\lambda K_{\sigma }\left ( \sigma_{0} \right )$ (3.3)

$\sigma _{0}$ 为位移结果对应的应力，将式(3.3)带入式(3.2)后可得：

$\left \{ \Delta P \right \}=\left ( \left [ K_e \right ] +\lambda \left [ K_\sigma (\sigma_0 ) \right ]\right )\times \left \{ \Delta u \right \}$ (3.4)

当结构承载达到临界载荷时，一个很小的载荷增量 $\Delta P$ 就可以产生较大的变形 $\Delta u$ ，因此当 $\Delta P\approx 0$ 时，上式变为经典的特征值问题：

$\left ( \left [ K_e \right ] +\lambda \left [ K_\sigma (\sigma_0 ) \right ]\right )\times \left \{ \Delta u \right \}=0$ (3.5)

保证上式有解，则矩阵的行列式为0：

$\left | \left [ K_e \right ] +\lambda \left [ K_\sigma (\sigma_0 ) \right ] \right |=0$ (3.6)

由计算出的 $\lambda$ 最小值给定临界载荷：

临界载荷=初始载荷×λ(临界载荷系数)

实例1：

以一端固定一段自由的压杆模型为例，杆件长500mm，截面尺寸5mm×5mm，弹性模量为2e5MPa，密度为7.85e-6 kg/mm3，初始载荷设置为1N，Y轴负方向。

理论计算结果： $F_{cr}=\frac{\pi ^2EI}{4l^2}=102.8N$

(a) 一般情况：建立静态结构和特征值屈曲的分析流程，如下图所示，具体流程可下载分析文件后查看。算例中提取了10阶模态，我们只需要第一阶的负载乘数即可，图中可以看到第一阶负载乘数为102.8。已知初始载荷为1N，因此临界载荷为102.8N，与理论计算结果相同。

但是，如果模型中存在恒定载荷(如重力)，或者需要考虑非线性问题时，临界载荷的计算便不能沿用上述公式了。

(b) 存在恒定载荷时：对上述模型施加重力载荷，其余参数不变进行特征值屈曲分析，得到第一阶负载乘数为79.697。可以计算出模型的质量为0.098kg，因此重力为0.96N。如果仍按照上面的公式计算临界载荷，那么载荷值应为79.697+0.96=80.657N，显然结果是不正确的。

因此，对于存在恒定载荷的模型，必须控制那些可变载荷使得第一阶负载乘数的计算结果接近1，此时在分析中所加的总载荷就是临界载荷。将外载荷由1N变为102N时，第一阶负载乘数为1.005，此时的总载荷为102+0.96=102.96N，即为临界载荷。

在求解完成后再次出现了一条提示信息，告诉我们此时临界载荷的计算方式改变了，变为

临界载荷=初始载荷×(负载乘数+1)

此时第一阶的负载乘数为101.8，则临界载荷=1×(101.8+1)=102.8N

3.2 非线性屈曲分析

非线性屈曲可以分为前屈曲和后屈曲，上一节特征值屈曲计算的是蓝色虚线上面的分枝点，在非线性屈曲中，以极限载荷点为界，前面的部分，即只计算极限载荷而不关注整体载荷-位移曲线的分析，称为前屈曲分析。追踪整体载荷-位移路径的分析称为后屈曲分析。

在上一篇几何非线性的文章中，使用牛顿法并分步加载求解一系列的平衡方程，即在每个增量步内按给定的载荷增量(载荷控制)或给定的位移增量(位移控制)，迭代求解找到系统的平衡位置。但是该方法在非线性屈曲分析中并不完全适用。如下图所示，在b-d段随着载荷的增加位移量出现了突变，这时系统的刚度是负值，牛顿法无法求解，因此载荷控制不适用，若使用位移控制可求解。但是有些情况位移控制也无法求解，如下图中的f-h段。为解决这些问题，在分析结构的后屈曲时，需要借助一些别的手段，配合着载荷控制或位移控制进行求解计算，这部分在后屈曲中进行详细描述。

3.2.1 前屈曲分析

前屈曲分析的基本方法就是，逐渐增加载荷直到结果开始发散，但是要注意，当施加的载荷快要接近临界载荷时，要保证载荷子步足够小，载荷子步太大的话会影响求解精度。

同时也要注意，一个未收敛的解可能并不意味着找到了最大载荷，也可能是由于数值不稳定导致的不收敛，此时可根据载荷-位移图表进行判断，当出现物理不稳定时，不稳定点切线即刚度是趋向于0的，此时结果可信。

实例2：考虑缺陷的压杆模型

为什么要考虑缺陷呢，如果结构是理想的，如实例(1)中的压杆稳定模型，理想压杆在承受轴向压力的情况下是永远不会发生弯曲的，不管压力多大所产生的变形只有轴向压缩。与此同时，在实际生活中几乎所有的结构或多或少都会存在初始缺陷，但要注意，只有在不设置缺陷就无法进行非线性屈曲分析的情况下，如理想压杆等分枝屈曲问题，比例因子才是必须的。如果一个结构在发生屈曲前发生了材料非线性变形或者大变形，如受压的壳体结构，此时就无需再引入初始缺陷，因为这时结构已经很容易失稳了，对初始缺陷并不敏感。

在实例1(a)后面再接一个静态结构，将特征值屈曲的“求解”数据与静态结构的“模型”相连，之后右击“求解”展开属性界面，便可以在“比例因子”处设置缺陷系数。图中的比例因子20，模式(模态)1，可以理解成将第一阶模态的变形放大20倍作为非线性屈曲分析的初始缺陷。

从下图可以看，设置比例因子后模型不再是完全竖直的，而是存在了一定的初始弯曲，这个初始弯曲可以理解为杆件在加工制造时存在误差(或公差)导致的。本文中将比例因子设为20是为了方便演示，当然在实际生活中零件加工时不会存在如此之大的误差。

比例因子的计算公式如下：

比例因子=零件的公差÷特征值屈曲的第一阶模态值

特征值屈曲的第一阶模态值一般都在1左右，因此比例因子可以由零件的制造公差确定。

定义初始缺陷后，将压力设为100N，打开大变形开关进行求解。绘制支反力和位移的图表，在79.125N处结构的变形开始逐渐增大，此时的力为临界载荷。该算例比较简单，也没有考虑材料的非线性，因此计算结构并未出现发散。

3.2.2 后屈曲分析

后屈曲分析是非线性屈曲分析的延续，由于达到极限载荷后，进入刚度下降段，这一段给数值计算带来了困难。在上一节简单提到过，载荷控制只能用于前屈曲问题，位移控制可以用于前屈曲和没有snap back的后屈曲问题。在某种情况下，如分析薄壳受压问题，单纯的载荷控制和位移控制是无法解决的，这时需要借助一些其他方法。

ANSYS有一些方法是专门用于解决后屈曲问题的，如非线性稳定性控制和弧长法。

非线性稳定性控制

载荷控制+重启动+非线性稳定性控制可以用于解决后屈曲问题，稳定性控制一般用于重启动之后，例如：对一个压力容器进行失稳分析，可以先设置大的时间步，打开大变形和重启动，求解计算。当施加的载荷大于失稳载荷时，不出意外计算的结果是不收敛的，之后缩小时间步，打开稳定性控制，选择发散点之前的稳定载荷步，从这个时间点开始重新计算。

稳定性控制用于处理全局屈曲和局部屈曲问题，可以配合线搜索和自动时间步使用；

稳定性可以理解成：在每个单元的节点上设置人工阻尼，该阻尼有两个节点，一个是单元的节点，另一个节点与地面相连，当结构发生较大变形或出现位移突变时，两个节点之间的假定速度(并非真的速度)也会发生突变，此时产生的阻尼力可以使结构保持稳定。阻尼系数可以通过设置能量耗散率获得(设置能量耗散率后软件自动计算)，也可以通过人工手动设置。

能量耗散率=稳定力做所的功÷单元势能，数值介于0到1之间。该数值越大，则稳定力越大，系统越容易稳定，但并不是说稳定力越大越好！

激活稳定性之后会出现“常数”和“减少”两个选项，“常数”表示阻尼系数在整个求解过程中保持不变，“减小”表示阻尼系数在整个过程中线性减少至0。使用场景如下，当存在多载荷步时，第一个载荷步若使用常数，则在第一个载荷步结束时，会有一些稳定力保留下来，如果稳定力没有变成一个很小的值，那么当第二个载荷步不使用稳定性控制时，它的第一个子步会出现收敛困难，原因是上一个载荷步遗留的稳定力突然变为0。因此，在这种情况下，第一个载荷步应该使用“减少”，保证载荷步计算完成后稳定了衰减至0。

一般情况下，结构初始状态是稳定的，因此当载荷子步设置合理时，第一个子步也应该是收敛的。激活第一个子步一般为“否”就好。

稳定力或者能量过大也会影响计算精度，虽然会出现提示信息，但是也需要在后处理对稳定力或稳定能进行检查：

如果稳定能远小于单元势能，则结果可以接受，无需再检查稳定力；

若稳定能较大，可检查所有子步下每个自由度的稳定力，如果稳定力远小于加载的载荷或支反力，则结果可以接受；

或者修改能量耗散率或阻尼系数，不要使稳定力过大。

注意事项：

Tips1：稳定性控制可用于解决不稳定的非线性问题，如后屈曲，跳跃点分析和不稳定材料。稳定性控制多用于重启动之后，也可以在开始计算时就打开；

Tips2：使用稳定性控制时，不要将最后一个收敛的子步作为重启动点，应使用前面的收敛子步，因为程序需要用一个子步去准备稳定性计算的数据；

Tips3：使用稳定性控制时，不同的子步数量可以产生不同的计算结果，原因是子步数量不同，节点之间的假定速度就不同，进一步造成稳定力的不同。所以结构越稳定，子步数量产生的影响越小。

弧长法

弧长法有较好的下降段计算能力，更容易解决后屈曲问题，但弧长法仅限于解决具有斜坡加载的全局屈曲问题，且不能配合线搜索，自动时间步和位移收敛准则使用。在Workbench中，如果要使用弧长法还需插入command，关于弧长法的详细介绍可参考ANSYS的help文档，或网上其他资料。如果要借助弧长法求解后屈曲问题，建议使用ANSYS APDL或ABAQUS软件。

实例3：

后屈曲分析直接使用静态结构，模型如下，圆弧半径25mm，Z向拉伸长度50mm，厚度1mm，材料为非线性结构钢，顶部线施加Y轴负方向的集中力2000N，底部两条边线简单支撑。

(a) 下图为使用非线性稳定性控制得到的载荷-位移曲线，先对模型进行求解，报错后选择发散点之前的子部作为重启动点，打开稳定性控制，再次求解。建立支反力和位移的图标，可以看出，该结构的临界载荷为1219.9N，后屈曲路径如图中所示。

(b) 使用弧长法进行求解，关闭自动时间步，位移收敛，线搜索和稳定性，插入command，输入：

nsubst,100            
arclen,on,5,0.005     
arctrm,u,52,ntip,uy   

##设置100个子步
##打开弧长法，最大弧长系数5，最小弧长系数0.005
##使用uy方向的位移作为收敛判据，大于50mm时计算结束

##最大(小)弧长=最大(小)弧长系数×施加的集中力÷子步个数

求解后可以得到载荷-位移曲线，图中可以看出，结构的临界载荷为1217.5N。

使用位移控制可以得到同样的结果。

附录中给出了Williams双杆模型的算例，分别使用载荷控制和弧长法进行求解。从结果图中可以看出，使用载荷控制得到的载荷-位移曲线无法显示出后屈曲路径中的snap through部分，该部分用直线进行了替代。相反的，使用弧长法得到的载荷-位移曲线中，可以清晰的看出snap through部分。

附录：文中所用案例

实例1：点击下载

实例2：点击下载

实例3：(a)点击下载 (b)点击下载

双杆算例：载荷控制下载弧长法下载

YORUNOKUNI

关注

17
点赞
踩
51

收藏

觉得还不错? 一键收藏
3
评论
几何非线性-屈曲分析

结构系统的破坏分为两种，一类是因材料软化而引起的强度破坏，另一类是几何软化使得结构的几何形状发生极大变化，从而使结构丧失承载能力。在材料力学中是这样定义压杆失稳的：中心受压直杆(理想模型)在直线状态下平衡，由稳定平衡转化为不稳定平衡时所受的轴向压力的界限值，称为临界压力。中心受压直杆在临界压力的作用下，其直线形态的平衡开始丧失稳定性，简称为失稳。显然，这个定义只使用于理想的压杆模型，那么关于失稳和屈曲是否存在一般性定义呢？
复制链接

扫一扫