浅谈屈曲与稳定，理论结合案例来看ANSYS workbench的线性屈曲分析【勘误版】
屈曲（Buckling）是指结构在受到压缩或剪切等载荷时，突然发生几何形态的失稳变形，导致承载能力显著下降的现象。它是结构稳定性失效的一种典型表现，常见于细长杆件、薄壁构件或受压结构。
在生活中，我们也经常看到屈曲现象，如下所示的薄壁罐子受压屈曲(图1)、某刚构件在荷载作用下屈曲(图2)。在工程中，也偶有钢结构的屈曲失稳造成工程事故，如图3。
图2 某钢构件屈曲
图3 某筒仓结构屈曲破坏
在结构工程中，结构或构件的屈曲不稳定是一种非常危险的现象，荷载的微小增加就可能导致灾难性的破坏。因此，在进行结构件设计时，除了满足强度、挠度要求外，结构还必须满足稳定性要求。今天咱们就来聊一聊屈曲与稳定!
--01--
屈曲与稳定，什么是分支点失稳? 什么又是极值点失稳?
我们知道，屈曲是结构因稳定性丧失而发生突然的几何变形（如弯曲、扭转或弯扭组合变形），导致承载能力急剧下降的现象。屈曲的本质是结构的平衡状态在外界扰动下发生了不可逆跳跃，从初始稳定形态转移到新的稳定形态。那么，在这里“稳定”又代表了哪些内涵呢？
在一个结构中，通常稳定性（Stability）是指结构在受到扰动后维持原有平衡状态或恢复平衡的能力。它是衡量结构抗干扰性和安全性的核心指标，通常以不稳定（失稳或屈曲）来衡量。那么，不稳定又是怎么表现的呢？
从过往的工程结构和试验中结构或构件失稳(屈曲)的现象来看，不稳定通常表现为:在载荷没有实质性变化的情况（存在微小的载荷扰动和几何缺陷），结构的位移发生非常大的变化。
一、分支点失稳(第一类稳定问题)
在结构的失稳(或屈曲)分析中，我们常用分支点（bifurcation）来描述这种不稳定的表现，即分支点为载荷历程中的一点，这一点代表着两个平衡路径的交点，表征屈曲失稳的萌生位置，如图4所示。
图4 平衡路径的分支点
从图4中，分支点后的表现可以看出，结构的载荷路径在越过分支点后可能表现出稳定、中性稳定和失稳。整个结构稳定性概念可以用图5来说明:
① 小球在图中1~2弧线内是稳定的，当有微小扰扰动时， 小球依然会返回到初始位置；
② 小球在图中2~3直线内是中性稳定的，扰动时，小球会保持在一个新的直线位置，荷载 F = Pc，此时载荷F即为临界载荷；
③ 小球在图中3点上是失稳的，扰动时圆会向下滚落。小球在图中3~4线内表现为快速通过(大变形)，跳跃到另一平衡位置，称为后屈曲。
从图5中可以看到，当荷载达到临界荷载时，结构的平衡状态发生突变，出现与原平衡形式完全不同的新变形模式，且平衡路径在临界点处分支（如理想轴压杆的欧拉屈曲），我们称为第一类稳定问题。
此时，结构既可在初始位置平衡，也可在偏离后的新位置平衡，具有平衡的二重性。临界荷载称为屈曲荷载，通常通过线弹性特征值法计算，即我们通常所说的线性屈曲。
线性屈曲的典型特征是:失稳前后变形性质改变（如直杆突变为弯曲状态）；且平衡路径发生分叉，通常数学上表现为特征值问题。
分支点失稳是在理想条件下求得的，通常假设:材料线弹性、无初始缺陷，通过欧拉公式或特征值屈曲分析求解临界荷载。例如，两端铰支轴压杆的临界荷载为：
在工程中，其意义在于给我们提供了理论上的稳定性上限，但因忽略实际缺陷，与临界荷载比现实情况偏大。
二、极值点失稳(第二类稳定问题)
实际结构中，材料存在非线性(如钢材的弹塑性)，其构件件也存在初始挠度等缺陷，因为扰动和非线性行为，即使构件受到的载荷低于临界载荷，结构也会变得不稳定，如图6所示。
这种由材料塑性或几何初始缺陷（如杆件初弯曲、偏心荷载）引发的失稳，我们称之为极值点失稳(第二类稳定问题)，其临界荷载称为极限荷载或压溃荷载，它需通过非线性分析方法求解。
图6 结构屈曲类型
从图6可以看到，有缺陷的结构在荷载增加过程中未发生变形模式的突变(未达到分叉点)，但力-位移曲线存在极值点（峰值荷载）。超过该极值后，荷载下降，变形持续增大直至结构压溃。
由图6还可以看到，以理想线弹性结构的理论屈服强度(分叉或分支)为界，将屈服分为前屈曲和后屈曲。前屈曲分析主要表现为线性特征值屈曲（也可以进行非线性屈曲分析），而后屈曲分析表现为考虑非线性因素的非线性屈曲（包括有缺陷结构的理想载荷路径的非线性屈曲和塑性行为、接触、大变形响应的非线性屈曲)。
极值点失稳的典型特征是: 失稳前后变形性质不变（如梁持续弯曲直至压溃），且实际工程中极值点失稳的临界荷载通常低于分支点屈曲荷载。
线性屈曲的分析——线弹性特征值法
线弹性特征值屈曲分析方法通过提取使线性系统刚度矩阵奇异的特征值获得结构的临界失稳载荷及失稳模态。
其推导过程如下:
① 由于线性屈曲分析时，结构状态位于前屈曲，因此满足线弹性的载荷-位移方程:
其中，[Ke] —为刚度矩阵，由此可得载荷{F0} 的位移结果{u0 } ，再由此得到对应的应变和对应的应力{σ} 。
② 假设前屈曲位移较小，可以得到载荷、位移和应力的增量方程：
式中，[Kσ(σ)] 为{σ} 应力状态下的初始应力矩阵。
③ 由于前屈曲状态下，载荷可以认为是一线性函数，即：
可得到：
④ 将上式带入②中的载荷、位移和应力的增量方程，可得：
⑤ 根据失稳的定义，在载荷变化很小时，结构将产生一个大的变形 {∆u} ，即{∆F}=0，上述④中方程变为：
其中：λ—为屈曲载荷因子（特征值）；
{∆u} —为屈曲模态形状（特征向量）。
该方程的意义在于：在n个自由度的有限元模型中，方程求得的λ的n阶多项式，此刻的{ ∆u} 表示屈曲时叠加到系统的变形，再由λ的最小值得到弹性临界载荷Fcr 。
总结: 线性特征值屈曲分析忽略了各种非线性因素和初始缺陷对屈曲失稳载荷的影响，大大简化屈曲分析，提高了屈曲失稳分析的计算效率，而且计算的特征值对结果稳定性评价有一定帮助。例如，当求解出密集排列的数值相差不大的特征值时，就表明该结构对缺陷敏感。
由于线性特征值屈曲分析基于线弹性的假设，得到的失稳载荷可能与实际相差较大，从特征值分析失稳，只能得到描述结构失稳时各处相对位移变化的大小，不是真实变形，或称为失稳模态，无法得到失稳后结构最大位移，但是失稳模态的形状可以作为非线性屈曲分析的初始几何缺陷。
线性屈曲案例分析
从本文前述内容，我们已经了解了屈曲与稳定的概念，也对线性屈曲分析的理论分析过程有了理解，接下来本文以一个实际案例讲解线性屈曲分析在ANSYS workbench中的分析过程。对于非线性屈曲，我们将在后续的文章中做深入的分析。
一、建模及材料参数
① 杆件尺寸及材料参数
本案例以一个支吊架C形槽钢为例，杆件采用欧标EN 1326 S250钢材(密度ρ=7850kg/m3,屈服强度fy=280MPa，弹性模量E=210GPa,剪切变形模量81GPa)，杆长为1m，截面尺寸详见图示。
图7 某装配式支吊架用C型槽钢截面
② 模型创建
利用ANSYS中的spaceclaim进行模型的创建。此次分析中，杆件的截面高度H=52mm，壁厚t=2.5mm,因为(5～8)＜H/t=25＜(80～100)范围为薄板,因此可认为杆件为薄壁杆件，拟采用壳单元进行模拟。在spaceclaim中建立实体模型后，通过抽取中间面形成壳单元所需的面体，如下图所示。
图8 1m杆件抽中面模型
③ 材料参数设置
杆件采用钢材：密度ρ=7850kg/m3,屈服强度fy=280MPa，弹性模量E=210GPa, 泊松比0.3，剪切变形模量81GPa，切线模量取2100MPa。ANSYS中材料参数输入如下图所示。
图9 材料参数设置
二、分析模型设置
① 单元选取和网格划分
在ANSYS中，单元选取和网格划分对求解的收敛和结果的准确性至关重要。ANSYS中，壳单元有Shell181和Shell182单元，其中 Shell181采用线性多项式作为形函数，Shell182采用二次多项式作为形函数。本次模拟中试件在进行单元网格划分时，为减小计算成本，采用Shell181单元。同时，保证模型准确和使模型能够收敛，网格划分为较为细密，采用2mm的网格。构件的网格划分见图10。
图10 网格划分结果及网格的Aspect Ratio检查
Aspect Ratio:长宽比，最佳为1，即正方形和正三角形。1~5 较好，结构分析时，为确保质量必须小于20。本例中，网格长宽比最大1.42左右，长宽比网格质量较好，可满足结构分析要求。
图11 网格的Jacobian Ratio检查
Jacobian Ratio: 此次网格划分雅可比大部分在1左右，单元的最大雅可比2.88，结构分析必须小于40。网格划分雅可比满足结构分析要求。
② 约束和荷载施加
本案例中设构件两端铰接，通过位移约束模拟铰接，约束设置如下：
1）杆两端约束X和Y方向的平动位移；
2）两端施加沿Z向的1000N的压力，为避免Z轴方向出现刚体运动，打开弱弹簧，设置结果如下图所示。
图12 约束施加
图13 荷载施加
需要注意，本文之前采用软件中的Simply Supported作为简支约束，但该约束是存在一定问题的，《使用梁单元和壳单元分析简支梁》一文中对该约束的问题进行了分析，因此本文修正了之前的内容
。
软件中的Simply Supported约束了X,Y,Z三个方向的平动，对于实体单元构件施加此约束相当于固定约束，壳单元也出现了类似情况。本次修改中，Z方向并未限制约束，而是通过施加一对平衡的力，并打开弱弹簧防止Z向刚体运动。这么做的目的在于避免形成Simply Supported同样的约束状态（尝试后，不正确）。
③ 线性静力分析及线性屈曲设置
完成上述设置后，我们可以在静力学模块中完成静力学的分析，结果如下图14。之后我们将“Toolbox(工具箱)”中的“eigenvalue Buckling(特征值屈曲分析)”命令直接拖曳到项目A(静力分析)的A6栏的“Solution”中，如下图15。
图15 拖入屈曲分模块进行数据传递
此时，更新完成后数据已经全部导入项目 B 中，之后双击项目 B中B5栏的“Setup”命令即可直接进入Mechanical界面进行线性屈曲分析。
重新求解分析树下的静力学模块的solution。求解完成后，选择“Outlines”中的“Linear Buckling(B5)”-“AnalysisSettings”命令，在下面出现的“Details of“Analysis Settings’”选项的“Options中设置“Max Modes to Find”栏中输入“10”，表示10 阶模态将被计算。设置完成后，点击求解，进行线性屈曲分析。
三、线性屈曲结果及分析
图16为分析的结果，在上图中，右下方的选项组中可以查到第一阶屈曲荷载因子为 49.085。
图16 分析结果
由于我们施加载荷（静力学荷载=摄动荷载）为1000 N，故可知杆件的屈曲压力为：
屈曲荷载=荷载因子×摄动荷载=1000×47.66N ≈47.66kN。
第一阶临界载荷为47.66kN，由于第一阶为屈曲载荷因子的最小值，因此这意味着在理论上，当加载荷载达到47.66kN时，结构将失稳。
此外，根据仿真计算的变形结果可知，第一模态时，构件存在明显的扭转变形，失稳属于弯扭失稳，与薄壁单轴对称截面的失稳形态相符。
此外，根据某产品手册中提供的受压容许承载力为51.11kN，如图17所示，相差（51.11-47.66）/47.66=7.2%，但仿真结果却比厂商提供的数值要小，这是为何?
图17 某厂商产品受压承载力表
笔者分析，这可能是由于产品在做分析时并未考虑杆件上开洞的影响，导致结果的差异。此外，网格划分差异、约束设置方式均可导致结果的差异。
屈曲失稳的临界荷载很大程度受到杆件约束的影响。该案例中，笔者也尝试过其他铰接约束设置方法，一旦设置了Z向约束，即使是某一个端点被设置了Z向约束，就会造成结果的较大差异。而采用弱弹簧的方式，结果与实际最相符。弱弹簧的影响从force reaction结果来看，影响很小，如下图所示：
当关闭弱弹簧后，会发生了轻微刚体 位移，但屈曲载荷因子没有变化。
关于线性屈曲案例的一些探讨
对于薄壁开口截面，理论上由于其剪切中心和形心并不重合，因此其轴心受压时其屈曲失稳形式一般为弯扭失稳。那么我们能否通过软件得到构件绕主轴的理论上的线性屈曲临界载荷呢？这是可以的。此次分析讨论，结合之前的分析，我们还考虑一下截面削弱的情况。
根据临界荷载欧拉公式计算，杆件两端铰接时，取计算长度系数为1.0，如下图18。不考虑截面削弱时，查询第一主轴毛截面惯性矩为Ixx=139412.0285mm^4（通过spaceclaim获得），如下图19。
由此，可计算出绕第一主轴Ixx的临界失稳荷载理论值：
Fcr=π2EI/(μl)2=3.142×210×109×139.412×10-9÷(1×1)2=288.7kN
图20 约束杆件X方向的位移
由此，重新计算，得到绕第一主轴Ixx的临界失稳荷载的可仿真分析结果：
图21 绕第一主轴Ixx的屈曲模态
仿真计算的绕第一主轴Ixx的屈曲荷载因子为276.14，则临界失稳荷载为：1000×276.14N ≈276.1kN。理论计算与仿真值相差约（288.7-276.1）/288.7=4.4%，误差在5%以内，说明通过仿真估计杆件临界荷载可行。但注意，仿真分析时我们的杆件实际上再腹板出有开洞情况。这有可能是造成其计算数值与理论值小的原因。
图22 杆件净截面参数
考虑截面削弱时（净截面），第一主轴净截面惯性矩为Ix=120835.7972（mm4）。计算其理论的临界失稳荷载为：
Fcr=π2EI/(μl)2=3.142×210×109×120.84×10-9÷(1×1)2=250.2kN
可以看到：
杆件无开洞按毛截面计算的理论临界失稳荷载为：Fcr1=288.7kN
杆件开洞按削弱截面仿真计算临界失稳荷载为：Fcr2=276.1kN
可以看到，杆件开洞按削弱截面仿真计算的结果在毛截面计算和净截面计算结果之间，符合实际情况。
因此，对于这种开洞后造成截面削弱的杆件，其真实的截面参数，如惯性矩我们可以通过反算的形式得到。例如本案例中，我们可以反算得到杆件腹板开洞后的计算截面惯性矩：
Ixx=Fcr(μl)^2/(Eπ^2)=133348.4mm4
仿真分析的主轴失稳临界荷载结果（276.1kN）与考虑净截面的主轴失稳临界荷载结果（250.2kN）相差25.92kN！说明完全按照开孔后的净截面计算，结果是偏安全。但全按照开孔后的净截面计算没有考虑腹板孔与孔之间板件对稳定性起到的作用。相对地，可以想到杆件截面如果开孔过多，杆件的稳定性将会受到削弱。
此外，开口薄壁型钢的临界失稳为弯扭失稳定，不是绕第一主轴的的失稳，弯扭失稳的临界荷载相比绕第一主轴失稳的临界荷载小很多。本例中，两端铰接时，杆件的弯扭失稳临界荷载为49.09kN，仅为绕第一主轴的失稳临界荷载276.1kN的17%左右。因此，在设计中，应尽量避免杆件发生弯扭失稳，对荷载、稳定性要求高的杆件，应当采用闭口对称截面，以提高材料的利用率！
