有限元理论基础及Abaqus内部实现方式研究系列48：屈曲分析（1）-理论 原创

（原创，转载请注明出处）

1 概述

本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过

（1） 基础理论

（2） 商软操作

（3） 自编程序

三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。

有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。

一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元软件iSolver，通过自研CAE软件和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。iSolver包括完整的前后处理和有限元求解器，功能如下，有兴趣可直接在下面网址下载：

百度网盘链接: https://pan.baidu.com/s/10d6jHdZ01SBY2JxiS6bffw 提取码: 6fdf

2 屈曲分析

结构失效的方式有两种：1. 材料的失效，也就是强度的破坏，2.几何的失效，也就是结构件的失稳，也叫屈曲。在系列文章4：非线性问题的求解中我们提到了弧长法求非线性屈曲，但非线性屈曲只是屈曲分析的一种方法，在工程上更多是一种简化方法。我们将分两个章节介绍结构屈曲分析，本章将介绍屈曲的理论，下一章节介绍屈曲的实验和模拟。

2.1 屈曲的含义

“屈曲”两个字中：

屈：大丈夫能屈能伸。屈和伸相反，表示结构件受压的状态。

曲：九曲十八弯，曲表示结构件弯曲。

屈曲的含义把“屈”和“曲”连起来就是结构件在受压时弯曲的现象，更具体点就是结构件在压力或载荷作用下，因稳定性不足而突然发生大变形甚至失效的现象。用大白话解释“屈曲”就是：结构件在压力下“突然塌掉”或“弯成奇怪形状”的现象。

屈曲一般发生在细长压杆或者薄板等结构件中。生活中有很多这样的例子，譬如帐篷的支架在大力下或者顶端放个重包突然失去支撑能力，导致帐篷坍塌，又譬如空的易拉罐用手指按压时，按压点会瘪下去，力比较小时，易拉罐外壳还能恢复，当指力足够大时，易拉罐外壳就直接现成一个永久的坑了。

工程上也有很多屈曲的现象，譬如火箭发射过程中，由于火箭加速产生的巨大轴向载荷（如发射时的加速度可达数十倍重力加速度），会导致细长结构（如箭体或储箱）因稳定性不足而弯曲，又譬如船舶在波浪中航行时，船体因波浪起伏产生巨大的纵向弯曲应力，导致中拱或中垂状态，引发肋骨或甲板的屈曲。

2.2 屈曲分析的目点

从以上生活常识可知，总有一个突然就垮塌的过程，就像上面的易拉罐例子一样，如果力比较小，放手后易拉罐还能恢复原样，我们认为这个力对这个结构件是安全的，但如果放手后产生了不可恢复的形变，那么这个力对结构件是不安全的，从安全到不安全的临界点的力就是临界载荷。

屈曲这个过程经历了：线弹性变形过程-》达到屈曲临界点-》（后屈曲）屈曲后的结构变形过程。

在工程上屈曲分析的主要目点是计算结构在轴向压力或弯曲荷载作用下发生屈曲失效的临界载荷值，从而判断当前设计是否安全。

2.3 屈曲分析的方法

屈曲分析有多种方法：

2.3.1 非线性屈曲分析

非线性屈曲分析是将力随着位移的关系表达出来，直到能看出哪点是临界载荷，临界载荷时位移增加时，力将不再增加，反而下降，也就是临界载荷就是载荷Vs位移曲线上的马鞍点位置。

如下面例子：

两个杆下方固定位移，中间往下加压力，每次力加的特别缓慢，确保每个时刻力达到平衡状态。

如果以往下的位移为横坐标，加的力为纵坐标，那么画出一条曲线大体如下：

力一开始随着位移增加而增加，知道顶点A，当过了顶点再往下压时，生活常识告诉我们不需要那么大的力也能往下压了，此时力随位移减小直到在水平位置的力变为0。采用有限元的非线性屈曲分析就是要寻找上述过程中的A点时的压力，有两种方式：

（1） 加强制位移约束，然后输出反力，做出反力和位移的曲线图，直接在图上查看马鞍点位置对应的反力大小。

（2） 上面方法缺点是实际工程我们不知道应该加哪个强制约束点或者如壳的屈曲，强制约束点太多无妨加和实际一样的强制位移，所以实际工程更多的是加力，做出力和位移的关系，此时由于要越过马鞍点，一般在有限元中采用弧长法（Risk）替代Newton迭代来计算，具体可看系列文章4：非线性问题的求解。

2.3.2 基于特征值的线性屈曲

非线性屈曲分析可捕捉整个屈曲从发生到后屈曲整个过程，同时如果变形量比较大或者材料达到塑性段时，那必须考虑非线性问题，只能采用非线性屈曲分析。但只要是非线性分析都有共同缺点，需要更大的计算量，而且非线性屈曲存在马鞍点时收敛更困难。

在实际工程上，很多问题不关心后屈曲的情况，只关心屈曲刚发生时刻点的临界载荷，而且有相当一大部分屈曲发生时材料还没达到塑性，变形也没超过5%，因此此时可以采用线性屈曲分析更快的得到临界载荷。

线性屈曲最通用的数值计算方法是基于特征值来求解。核心思想是如果我们已知时刻0和1的两个位移和刚度K0、K1，那么能不能得到t时刻的刚度呢？

结构有限元的刚度阵按照虚功原理得到。

虚功原理可以理解为外力在虚拟位移下做的虚功=内部应变能的一段小时间内对应变能的积分：

S和E分别表示应力和应变。

求导后得到任意时刻的刚度阵K如下：

也就是刚度矩阵将分为两块，上式的前面一部分依然是以前的BDB形式，只不过B换成了当前时刻的应变位移矩阵，而后面新增项将将转换为S*G刚度阵，称为几何刚度阵，也称为初始应力矩阵（initial stress stiffness）。

因为是小变形线性材料，所以B、D、G为常量，BDB我们可以认为是一个常量，也就是K0，则

• 因为压缩导致刚度减小，我们在上式中假定S和G是正数。

我们任意时刻外力为F，由于线性屈曲假定整个屈曲前过程都是线性问题，那么得到的应力S正比于位移U，也正比于F，设为

带入上式得到

(1)

上式就是一个一次方程。将K和F的关系是一条直线如下

其中，这条直线与y轴相交于已知点（0，K0），与x轴相交于未知点（Fe，0）。Fe即我们要求的临界载荷。

下面仅是怎么求出斜率AG，一种简单的方式就是在直线上找两个已知点就能求出斜率了。既然已经有一个已知点（0，K0），那么取时刻1作为另一个已知点（F1，K1）

2.3.3 基于欧拉应力理论修正的线性屈曲

非线性屈曲分析和基于特征值的线性屈曲看起来已经把有限元屈曲分析的所有情况覆盖了，但实际工程上很多行业还是采用基于欧拉应力理论的线性屈曲。

针对一根压杆的受压，如下图，左端简支，右端约束y、z位移且加压力F。设此压杆是完全弹性的，且应力不超过比例极限，若轴向外载荷F小于它的临界值Fe，此杆将保持直的状态而只承受轴向压缩。如果一个扰动(如—横向力)作用于杆，使其有一小的挠曲，在这一扰动除去后。挠度就消失，杆又恢复到平横状态，此时杆的直的形式的弹性平衡是稳定的。若轴向外载荷F大于它的临界值Fe，柱的直的平衡状态变为不稳定，即任意扰动产生的挠曲在扰动除去后不仅不消失，而且还将继续扩大，直至达到远离直立状态的新的平衡位置为止，或者弯折。此时，此压杆失稳或屈曲。

对此简单问题，我们猜测一下这个Fe大概为多少？

当到达Fe时，压杆开始便变形，根据生活常识，应该大体变形为如下形状：

显然当L足够小时，一定会超过材料屈服强度也会到时结构件失效。

实际工程材料因此如果将结构件失效应力和长度做一条曲线将会是如下形式

这条曲线在L>Ly时是双曲线，在L<Ly时是直线，且失效应力恒定为材料屈服强度。基于特征值的线性屈曲分析结合屈服分析就会得到这条曲线。

这个简单例子用实验很容易得到，下章我们也将做一个简单的实验（估计应该是全网第一个采用简单的家用工具做出来的屈曲试验）。做实验发现是如下修正曲线：

实际没有一个明显的转折点，很多行业规范会人为分为两段，在L>Lp时还是双曲线，而在L<Lp时采用抛物线，两者的相交点为(Lp, )。在船舶行业，取为/2。

这是一根压杆得到的曲线，模拟的最终目点还是和实验尽量接近，既然它比基于特征值的线性屈曲分析更接近试验，那么在实际工程中也更受欢迎。船舶行业的线性屈曲就采用基于欧拉应力理论修正的线性屈曲。长方形壳单元可以看成是压杆截面的一个维度取为实际平面尺寸的一个应用。同时，为了适用一般的壳形状，船舶行业的规范规定了三步的模拟：

（1） 先确定板格的位置，周围由桁材、纵骨或者不在一个平面的面板围出来的图形就是板格，如果是有限元模型，板格一般由多个板单元组成。如下便是一个板格的示例，每个不同颜色表示不同的板格：

（2） 板格等效为长方形的规范，上面示例的板格都是长方形，但实际问题很多板格都是不规则形状，如下分别就是一个带弧度的壳单元和三角形等效为长方形的规范示例。

（3） 套用基于欧拉应力理论修正的线性屈曲公式分析是否屈曲。

2.4 屈曲的本质

以下面简支梁为例：

模型尺寸：

长L=240，截面为10X5。E=1.5e6，打开几何非线性模拟实际情况。

（1）当不加任何力时，打印刚度阵如下，y方向K22=230：

（2）由上面欧拉公式得到欧拉屈曲临界力Fe=26773，一端施加压力Fe=26773，

考虑梁的几何非线性，得到的刚度阵如下，可以看到K11增加，但K22，K33等减少。K22=118.

（3）压力增大为Te=26773*2，得到刚度阵如下，K22继续减少为7.3，已经非常接近0了。

由这个例子，可以知道，杆件受加压载荷时，横向方向刚度随压力反比减小，此时如果有一个小的横向力，将产生大的位移，导致横向抗外力能力的下降，导致失稳弯折。这就是屈曲的本质。

• 几何非线性时，加了线性屈曲欧拉理论的压缩载荷，Y方向刚度减小了，甚至2倍时也不为0，所以线性屈曲的理论和实际还是有一些差异的，有兴趣可以再细致研究一下差异的来源（个人猜测需要研究考虑几何非线性的B导致的几何刚度阵的差异）。

2.5 总结

本文总结了屈曲分析的三种方法：

（1） 非线性屈曲分析

（2） 基于特征值的线性屈曲

（3） 基于欧拉应力理论修正的线性屈曲

在实际工程中，线性屈曲更适合于工程问题，同时，由于线性屈曲理论和试验的差异，部分行业实际上更多采用基于欧拉应力理论修正的线性屈曲。同时说明屈曲的本质还是纵向加压后横向刚度变小，导致横向抗力能力的下降，导致失稳弯折。

2.6 视频讲解和操作验证演示

如果觉得上面的文字太复杂，也可以看一下视频的简要讲解和软件演示，地址如下：

https://www.jishulink.com/college/video/c12884 20理论系列文章48-屈曲分析(1)-理论

3 以往的系列文章

3.1 ========第一阶段========

第一篇：S4壳单元刚度矩阵研究。

http://www.jishulink.com/content/post/338859

第二篇：S4壳单元质量矩阵研究。

http://www.jishulink.com/content/post/343905

第三篇：S4壳单元的剪切自锁和沙漏控制。

http://www.jishulink.com/content/post/350865

第四篇：非线性问题的求解。

http://www.jishulink.com/content/post/360565

第五篇：单元正确性验证。

https://www.jishulink.com/content/post/373743

第六篇：General梁单元的刚度矩阵。

https://www.jishulink.com/content/post/403932

第七篇：C3D8六面体单元的刚度矩阵。

https://www.jishulink.com/content/post/430177

第八篇：UMAT用户子程序开发步骤。

https://www.jishulink.com/content/post/432848

第九篇：编写线性UMAT Step By Step。

http://www.jishulink.com/content/post/440874

第十篇：耦合约束（Coupling constraints）的研究。

https://www.jishulink.com/content/post/531029

3.2 ========第二阶段========

第十一篇：自主CAE开发实战经验第一阶段总结。

http://www.jishulink.com/content/post/532475

第十二篇：几何梁单元的刚度矩阵。

http://www.jishulink.com/content/post/534362

第十三篇：显式和隐式的区别。

http://www.jishulink.com/content/post/537154

第十四篇：壳的应力方向。

https://www.jishulink.com/content/post/1189260

第十五篇：壳的剪切应力。

https://www.jishulink.com/content/post/1191641

第十六篇：Part、Instance与Assembly。

https://www.jishulink.com/content/post/1195061

第十七篇：几何非线性的物理含义。

https://www.jishulink.com/content/post/1198459

第十八篇：几何非线性的应变。

https://www.jishulink.com/content/post/1201375

第十九篇：Abaqus几何非线性的设置和后台。

http://www.jishulink.com/content/post/1203064

第二十篇：UEL用户子程序开发步骤。

https://www.jishulink.com/content/post/1204261

3.3 ========第三阶段========

第二十一篇：自主CAE开发实战经验第二阶段总结。

https://www.jishulink.com/content/post/1204970

第二十二篇：几何非线性的刚度矩阵求解。

http://www.jishulink.com/content/post/1254435

第二十三篇：编写简单面内拉伸问题UEL Step By Step

http://www.jishulink.com/content/post/1256835

第二十四篇：显式求解Step By Step。

https://www.jishulink.com/content/post/1261165

第二十五篇：显式分析的稳定时间增量。

http://www.jishulink.com/content/post/1263601

第二十六篇：编写线性VUMAT Step By Step。

https://www.jishulink.com/content/post/1266640

第二十七篇：Abaqus内部计算和显示的应变。

https://www.jishulink.com/content/post/1273788

第二十八篇：几何非线性的T.L.和U.L.描述方法

https://www.jishulink.com/content/post/1282956

第二十九篇：几何非线性的T.L.和U.L.转换关系

https://www.jishulink.com/content/post/1286065

第三十篇：谐响应分析原理

https://www.jishulink.com/content/post/1290151

3.4 ========第四阶段========

第三十一篇：自主CAE开发实战经验第三阶段总结

https://www.jishulink.com/content/post/1294824

第三十二篇：谐响应分析算法

https://www.jishulink.com/content/post/1299983

第三十三篇：线性瞬态动力学

https://www.jishulink.com/content/post/1302074

第三十四篇：非线性瞬态分析

https://www.jishulink.com/content/post/1787283

第三十五篇：接触求解算法

https://www.jishulink.com/content/post/1792869

 第三十六篇：DLOAD用户子程序开发步骤

https://www.jishulink.com/content/post/1826803

 第三十七篇：梁单元差异（1）-理论基础

https://jishulink.com/content/post/1872208

 第三十八篇：梁单元差异（2）-梁截面方向

https://www.jishulink.com/content/post/1874628

 第三十九篇：梁单元差异（3）-剪力和弯矩

https://www.jishulink.com/content/post/1876013

 第四十篇：梁单元差异（4）-形心、剪心和偏置

https://www.jishulink.com/post/1888000

3.5 ========第五阶段========

 第四十一篇：自主CAE开发实战经验第四阶段总结

https://www.jishulink.com/post/1905963

 第四十二篇：声学分析(1)-有限元

https://www.jishulink.com/post/1912491

 第四十三篇：声学分析(2)-边界元

https://www.jishulink.com/post/1923936

 第四十四篇：声学分析(3)-湿模态

https://www.jishulink.com/post/1928692

 第四十五篇：约束关系（1）-统一形式

https://www.jishulink.com/post/1933077

 第四十六篇：约束关系（2）-Lagrange因子法求解

https://www.jishulink.com/post/1945262

 第四十七篇：约束关系（3）-船舶规范约束导致的Max Ratio问题

https://www.jishulink.com/post/1950396