基于ANSYS的轴心受压柱屈曲分析

lalalala

为了了解和掌握轴心受压柱特征值屈曲和非线性屈曲差异，以及考虑在屈曲分析中划分不同单元数量对分析结果的影响，选取适当的单元数量，利用有限元软件ANSYS对结构进行分析。初步了解特征值屈曲与非线性屈曲所得结果差异。在此基础上进行了多例轴心受压柱的仿真模拟分析，同时考虑不同长细比对屈曲分析结果的影响，掌握了长细比变化对轴心受压柱特征值屈曲和非线性屈曲的计算结果的影响规律。提出工程中应尽量采取非线性屈曲分析，并在分析中采取正确的分析方法。

引言：

随着计算机的发展人类实现了一个又一个的突破，大大提高了产品开发、设计、分析和制造的效率和产品性能。有限元理论的发展对于建筑专业更是一个飞跃。在结构线弹性计算中，一般都假定在加载过程中用结构变形前的形状来代替结构变形后的形状。然而结构在实际工程结构中，往往存在大位移、大转角或大应变等问题。这时的平衡条件就应如实的建立在变形后的形状上，以考虑变形对平衡的影响，因此要考虑非线性屈曲分析。在进行ANSYS分析时，如果单元数量选取不当，会使结果产生很大的误差，选取正确的单元数量是计算的前提条件。

一．   划分不同单元数量对特征值屈曲和非线性结果影响的分析

本节讨论特征值屈曲和非线性屈曲结果影响分析受单元网格密度的影响，通过分析时通过改变网格密度，所得计算结果提取第一阶特征值屈曲稳定系数和非线性屈曲系数。通过所得数据进行对比，当前后两个结果满足一定误差要求时，即可认为结果正确，否则应继续改变网格密度进行比较。最终找到本单元类型所需划分最佳的单元数量。

1.有限元模型参数

（1）单元类型：BEAM189

（2）截面尺寸：宽度B=0.05，高度H=0.05，长度L=5m

（3）材料属性：Q235钢， EX= 2.06×1011pa，泊松比NUXY=0.3，
      （4）划分单元数： 变量

（5）约束情况：上、下端均端为铰接

（6）分析类型：屈曲分析

（7）受力特征：上端集中力F=-1N。

2. 建立计算分析模型

3.ANSYS分析结果

表（1）分析结果

分析结果绘成曲线如图:

5.结果分析

通过beam189单元类型进行分析所得数据进行对比可知，当单元数量为100时，特征值屈曲和非线性屈曲前后两个结果已满足一定误差要求时，可认为结果正确。因此之后分析时选着单元数量为100。

二．   轴心受压柱特征值屈曲和非线性屈曲的ANSYS分析

在 ANSYS 中，稳定分析分为两类：线性特征值屈曲分析和非线性屈曲分析。

本节通过改变截面尺寸达到改变杆件长细比的目的，选取模型长度为l=5m，由于长细比（ 的变化只与截面刚度有关。本节针对7种不同的大柔度杆 进行分析，截面尺寸及截面惯性矩如表（2），截面1-4取截面惯性矩变化梯度为200cm4，截面4-7取截面惯性矩变化梯度为10cm4。最后把特征值屈曲和非线性屈曲所得结果进行对比，提出什么情况下可以选着用特征值屈曲什么情况下选择用非线性屈曲。

1.有限元模型参数

（1）单元类型：BEAM189

（2）截面尺寸：宽度B、高度H如表（2）所示，长度L=5m

（3）材料属性：Q235钢， EX= 2.06×1011pa，泊松比NUXY=0.3， =100

（4）划分单元数： 100

（5）约束情况：上、下端均端为铰接

（6）分析类型：屈曲分析

（7）受力特征：上端集中力F=-1N。

表（2）截面尺寸和刚度

2. 建立计算分析模型

如上图（1）

3.   理论分析

（1）特征值屈曲分析

有限元法对结构静力屈曲失稳问题的分析，对于解决线性屈曲问题，应用特征值计算方法。特征值算法是通过特征值分析计算屈曲载荷，该类屈曲分析主要是针对平衡临界状态的求解，其中包括临界载荷和屈曲模态的求解；按特征值分析屈曲、失稳临界载荷是一种简便的稳定性分析方法，可以获得平衡路径的分叉点。

对于受压结构，随着压应力的增加，结构抵抗横向变形力的能力下降。当载荷大到某一水平，结构总体刚度变为零，丧失稳定性。屈曲分析研究失稳发生时的临界载荷和失稳形态。基于结构失稳前系统刚度阵出现奇异，可将失稳问题转化为特征值问题处理。线性屈曲载荷的计算，属于结构小位移材料线弹性的屈曲范畴。对于总体 Lagrange 式的几何非线性的有限元方程可以写为：

其中 是与应变表达式中非线性应变相关的部分，而 是与应变表达式中线性应变相关的部分，是由于初始应力引起的，通常称为初应力矩阵。 是相关的外力项。另外 ，其中 为初位移刚度矩阵或大位移刚度阵， 为初应力刚度阵或几何刚度阵。

对于特征值稳定问题，载荷可以表示为 。其中 是载荷模式， 是载荷幅值。求解过程应该首先求解对应于载荷 的线性平衡问题

其中 是结构的线弹性刚度矩阵。从上式解得 ，进而可以得到结构内的应力分布 。结构临界载荷 ，可以通过求解关于 的特征值问题得到。

如果认为在结构初始失稳时，初始位移 仍然很小，则可以在有限元方程中忽略其影响，并且可以忽略大位移刚度阵 (1)式变为：

在总体Lagrange式中，将 代入上式，并考虑到结构达到稳定的临界载荷时，可认为 为 0，则得到下列方程：

这就是结构稳定的求解问题。要使（4)有非零解，则需保证

上式为一个广义的特征值方程，求解式(5)解得各阶特征值λ ，从而得到相应的其它物理量。

对于大柔度轴心受压杆件理论上采用欧拉临界力计算公式：

为欧拉临界应力，E为材料弹性模量， 为杆件计算长度系数， 为杆件实际长度。

线性特征值屈曲分析省略了非线性项，作为一种线性屈曲分析方法，是对理想弹性结构的理论屈曲强度的预测，满足于经典的解析理论。忽略了各种非线性因素和初始缺陷对屈曲失稳载荷的影响，使屈曲问题大大简化，从而提高了屈曲失稳分析的计算效率。但是，由于材料的缺陷和非线性，往往导致结构在理论弹性屈曲强度之外的点位发生屈曲。因此，线性特征值屈曲分析经常得到的是非保守结果，得到的失稳载荷可能与实际相差较大。通常情况下不能用于实际的工程分析。

来源：CAE技术联盟

lalalala

（2）非线性屈曲分析

    因几何变形引起结构刚度改变的一类问题都属于非线性问题。非线性通常分为大应变、大位移和应力刚化。以上三种大应变导致结构刚度变化的因素，即单元形状改变、单元方向改变和应力刚化效应。此时应变不再假定是“小应变”而是有限应变或“大应变”。

非线性屈曲分析采用几何非线性的荷载一位移全过程跟踪有限元分析。由能量原理可得到修正的拉格朗日(U．L．)形式的非线性增量有限元基本方程：

(6)式中， 为结构在状态的切线刚度矩阵， ，其中[KE]为结构的线弹性刚度矩阵； 为 次迭代时初应力刚度或称几何刚度矩阵(轴向力规定以拉力为正)，它考虑了单元内力对结构变形的影响； 为结构 次迭代时初位移刚度矩阵或称大位移矩阵，它考虑了结构位置变化对平衡的影响(或结构的变形对刚度的影响)。 为结构在 +1次迭代过程中位移增量列阵， 为 +1次迭代过程的荷载比例系数；{P}为初始选定不变的节点荷载向量； 为 次迭代时各单元内力等效的节点力向量。方程(1)的求解采用把弧长法(Arc．Length Method)和Newton—Raphson法相结合的增量迭代法。荷载增量采用弧长法自动加载。弧长法将荷载比例系数和未知位移同时作为变量，用曲线弧长来控制荷载步长，可使Newton．Raphson法平衡迭代沿一条弧收敛到其平衡路径，以避免矩阵在那些奇异点处变为奇异矩阵，从而避免了结构在加载时某些点可能出现的物理意义上的不稳定(即结构的荷载一位移曲线的斜率为零或负值)，并控制收敛性，帮助稳定数值求解。

4.   基于ANSYS分析

本节选取编号1截面型式进行详细分析，其他编号截面计算步骤相同只给出计算结果。

（1）特征值屈曲分析

选取编号1截面，分析方式为静力分析，并且打开预应力选项，求解。在列杆件屈曲方程时，都假定构件有了一定的侧向变形，预应力效应与此相似。打开预应力效应是把静力分析的结果产生的几何刚度加进去。选择求解方式为Block Lanczos，并且选择提取5阶屈曲模态，并且在载荷步选项卡中设定对 5阶屈曲模态进行扩展，求解。下面几幅图显示了不同模态的结果。

线性特征值屈曲分析所产生的多阶模态结果，直观来看，是对于线性屈曲计算产生的不同特征值所绘制的变形图，然而工程实际是不会对同一个结构产生多种屈曲的，当承载使其达到第一阶屈曲的载荷时，就会发生屈曲，因此分析时只提取一阶屈曲系数，作为实际工程中应用。

提取第一阶屈曲系数FQRT1= 677030

（2）非线性屈曲分析

屈曲问题主要分为两类：分叉点屈曲和极值点屈曲。前面提到的特征值屈曲问题，属于分叉点屈曲。ANSYS 模拟特征值屈曲问题时，对于理想压杆的线性特征值屈曲问题，可以很好的模拟；但是，对于非线性特征值问题，ANSYS 并不能给出让人满意的解答。但是，可以用解决极值点屈曲问题的方式，也就是压溃理论，去求解非线性特征值问题。

由于线性特征值屈曲分析仅限于线性问题，忽略了工程实际中确实存在的非线性项，所得的结果不够准确，所以在实际工程分析中，更多的是采取非线性屈曲分析的方式，解决结构的稳定性分析问题。下面，将在 ANSYS 中实现对理想轴压杆的非线性屈曲分析，这里将采用静力学结构分析的方式，使用 Newton－Raphson算法和弧长法(Arc．Length Method)对杆结构进行非线性屈曲分析。由于本论文采用Q235钢作为材料，极限应力

进入求解器，打开大变形选项卡，并且勾选预应力选项，设定迭代子步数为200，求解。得到非线性屈曲系数为552720。最大位移发生在跨中截面处，绘出跨中截面的荷载位移曲线如图。

（3）ANSYS屈曲分析结果

（4）结果分析

<1>特征值屈曲计算结果与非线性屈曲计算结果相差超出误差允许范围，所以在工程中应采取非线性屈曲分析的方式。

<2>随着截面惯性矩的均匀减小，而特征值屈曲系数与非线性屈曲系数差值越来越小，但差值率基本相差不大。

<3>线性特征值屈曲分析所产生的多阶模态结果，是对于线性屈曲计算产生的不同特征值所得的变形。然而工程实际是不会对同一个结构产生多种屈曲的，当承载使其达到第一阶屈曲的载荷时，就会发生屈曲，因此第一节屈曲是我们设计时参照的重点。

<3>当作用在梁顶端的载荷逐渐增大时，伴随着位移增大的同时，结构发生屈曲。对于理想弹性材料，结构的支反力随着跨中截面UX 方向位移逐渐增大，并且在加载达到一定数值时发生屈曲，曲线在此位置发生了明显的斜率变化在屈曲位置之后，曲线变得平缓，并且随着加载的增加，缓慢的增大。

<4>在求解非线性问题时解法的选择。ANSYS 对屈曲问题不同处理方式有：线性特征值屈曲法，Newton－Raphson 法，弧长法。三种方法各有优劣。线性特征值屈曲法前面已经提到，不适合解决工程实际问题。Newton－Raphson 法可以控制子步数，最为常用。用弧长法作屈曲分析时，一定使用应力刚化，这在非线性屈曲分析中，由于打开了大变形选项而自动满足。弧长法由于程序可以自动控制步长，所以计算速度快，并且可以更为精确的捕捉到极值点。弧长法可以控制最大步长，如果最大步长太大，可能使求解跨过临界载荷点；得到一个错误的结果。 也可以设定最下迭代步长，使得程序可以不断缩小步长至最小步长，来达成收敛。

三、结论

1 .基于ANSYS进行屈曲分析中，单元数量的选择对分析结果影响很大，通过所得数据进行对比，当前后两个结果满足一定误差要求时，即可认为结果正确，否则应继续改变网格密度得到结果进行比较。最终找到本单元类型所需划分最佳的单元数量，本论文进行屈曲分析时选取最佳单元数量为100。

2. 线性特征值屈曲分析所产生的多阶模态结果，是对于线性屈曲计算产生的不同特征值所得的变形。然而工程实际是不会对同一个结构产生多种屈曲，当承载使其达到第一阶屈曲的载荷时，就会发生屈曲，因此第一节屈曲是我们设计时参照的重点。

3. 随着截面惯性矩的均匀减小，而特征值屈曲系数与非线性屈曲系数差值越来越小，但差值率基本相差不大。

4.特征值屈曲计算结果与非线性屈曲计算结果相差超出误差允许范围，所以在工程中应采取非线性屈曲分析。

5.当作用在梁顶端的载荷逐渐增大时，伴随着位移增大的同时，结构发生屈曲。结构的支反力随着跨中截面UX 方向位移逐渐增大，并且在加载达到一定数值时发生屈曲，曲线在此位置发生了明显的斜率变化在屈曲位置之后，曲线变得平缓，并且随着加载的增加，缓慢的增大。

6. ANSYS 对屈曲问题不同处理方式有：特征值屈曲法，Newton－Raphson 法，弧长法。在选择求解非线性问题的解法时，我们要根据问题的需要采取最为正确的问题解法。

附：命令流：

/PREP7

B=0.1

H=0.1

L=5

E=2.06E11

p=-1

N=100

ET,1,BEAM189

MP,EX,1,E

MP,PRXY,1,0.3

SECTYPE,1,BEAM,RECT

SECDATA,B,H

K,1

K,2,,L/2

K,3,,L

K,10,0,L/2,L/2

L,1,2

L,2,3

LATT,1,,1,,10,,1

LESIZE,ALL,,,N

LMESH,ALL

FINISH

/SOLU

/VIEW,1,1,1,1  

/ESHAPE,1.0                              

DK,1,UX,,,,UY,UZ,ROTY

DK,3,UX,,,,UZ,ROTY

FK,3,FY,p

PSTRES,ON

SOLVE

FINISH

!特征值屈曲分析

/SOLU

ANTYPE,BUCKLE

BUCOPT,LANB,1

MXPAND,5

OUTRES,ALL,ALL

SOLVE

FINISH

/POST1

SET,LIST                                 

*GET,FREQ1,MODE,1,FREQ

FINISH

!非线性屈曲分析

/CONFIG,NRES,500

/PREP7

TB,BISO,1,1,2

TBTEMP,0

TBDADA,,2.0E8,0

UPGEOM,0.01,1,1,'gangwenyong1','rst',' '

FINISH                                    

/solu                                             

antype,0                                          

nlgeom,1                                            

outres,all,all                                      

arclen,1,0                                         

arctrm,l                                             

nsubst,500,,,1                                      

fk,3,fY,p*freq1                                    

solve                                               

finish               

/post26                                    

nsol,2,2,u,X,deflection                             

rforce,3,1,f,Y,reactionf                           

/axlab,x,deflection

/axlab,y,reactionf                                 

xvar,2                                                

plvar,3                                            

finish      

来源：CAE技术联盟