03 张弦柱稳定算例 (ANSYS)

《有限元法及其应用》，机械工业出版社，2006, 光盘演示算例

03 张弦柱稳定算例 (ANSYS)

本算例将讨论带有预应力的张弦柱的特征值失稳和非线性失稳问题

知识要点：
(a) 预应力
(b) 特征值稳定
(c) 考虑其他内力的特征值稳定
(d) 添加初始缺陷
(e) 弧长法
(f) 非线性分析和收敛
(g) 荷载位移关系

(1) 设定分析参数，在ANSYS顶部菜单Parameters->Scalar Parameters，在弹出的Scalar Parameters窗口中输入，FORCE=100，OFFSET=0.1。
(2) 建立模型，在本算例中，我们将接触两种新的单元，三维梁单元Beam 4和三维索单元Link 10。Beam 4单元在参数定义等方面比前面介绍的Beam 188单元简单，对于常见的矩形弹性截面，也是一个很实用的单元。Link 10单位为ANSYS提供的空间索单元，用户可以控制该单元只能受压或者只能受拉，默认该单元只能受拉。
(3) 在ANSYS主菜单Preprocessor->Element type->Add/Edit/Delete中选择Beam 4单元和Link 10单元
(4) Beam 4单元和Link 10单元都可以通过实参数来设置截面属性。首先设置Beam 4单元的截面。我们要建立的截面是一个0.1m×0.12m的矩形截面。Beam 4单元如果没有指定截面主轴方向，则截面局部坐标系的Y轴方向将和整体坐标系的X-Y平面平行。在ANSYS主菜单Preprocessor->Real Constants->Add/Edit/Delete，选择Add，指定实参数的关联单元类型为Beam 4单元，输入截面参数为AREA: 0.1*0.12, IZZ: 0.12*0.1**3/12, IYY: 0.1*0.12**3/12, TKZ: 0.12, TKY: 0.1, 如下图所示

(5) 继续添加第二个实参数类型，指定第二个实参数与Link 10单元相关。取Link 10单元的截面积为4×10-6m2，初始应变为2×10-3，这里正的表示初始应变为拉应变。

(6) 下面设定材料属性，在本例子中为了简单起见，所有材料都设定为钢材。在ANSYS主菜单Preprocessor->Material Props->Material Models中添加材料属性。在Material窗口选择Structural->Linear->Elastic->Isotropic，输入弹性模量2×109和泊松比0.27，还可以输入材料的密度，在Material窗口选择Structural->Density，输入密度为7800。
(7) 接下来建立几何模型，在本例子中我们还是严格按照ANSYS建模过程要求的点、线、面、体拓扑关系来建立几何模型。
(8) 首先建立关键点，在ANSYS主菜单Preprocessor->Modeling->Create->Keypoints->In Active CS中建立以下关键点

(9) 下面将关键点用线连接起来，在ANSYS主菜单Preprocessor->Modeling->Create->Lines->Lines->Straight Line中按次序将以下关键点连接起来

得到张弦柱模型如图

(10) 下面对几何体划分网格。在ANSYS主菜单Preprocessor->Meshing->Mesh Attributes->Picked Lines，选择1～6号直线，设定其单元类型、材料类型和实参数都是1。选择7～14号直线，设定材料类型为1，单元类型和实参数都是2。
(11) 下面要控制网格划分的尺寸。对于本次分析而言，中间的柱子是研究的重点，所以我们希望把网格划分得密一些，在ANSYS主菜单Preprocessor->Meshing->Size Cntrls->Lines->Picked Lines，选择1～6号直线，设定网格划分的单元最大长度（Element Edge Length）为0.3。对于周边的预应力索，由于Link 10是一个比较复杂的非线性单元，求解时候容易出现比较困难的收敛问题，因此我们对预应力索只划分成一个单元。选择7～14号直线，设定网格划分段数为1。
(12) 在ANSYS主菜单Preprocessor->Meshing->Mesh->Lines，选择所有的直线，划分网格
(13) 在ANSYS顶部菜单PlotCtrls->Style->Size and shape，设定Display of element为On。得到单元形状如图：

(14) 下面对模型添加边界条件和荷载。首先进入ANSYS主菜单Solution->Define Loads->Apply->Structural->Displacement->On Keypoints，选中关键点3，约束三个水平运动自由度UX, UY, UZ和转动自由度ROTZ，选中关键点2，约束两个平动自由度UX和UY。
(15) 进入ANSYS主菜单Solution-> Define Loads-> Apply-> Structural-> Force/Moment-> One Keypoints，选中关键点2，荷载方向为FZ，大小为-FORCE
(16) 先要进行一次静力分析，进入ANSYS主菜单Solution->Analysis Type->New Analysis，设定分析类型为Static，进入Solution->Analysis Type->Sol'n Controls，设置分析为小变形分析并考虑预应力，关闭自动时间步长控制并设定分析子步数(Substeps)为1

(17) 进行一次求解操作，进入ANSYS主菜单Solution->Solve->Current LS。
(18) 下面求解该模型的特征值屈曲。进入ANSYS主菜单Solution->Analysis Type->New Analysis，设定分析类型为Eigen Buckling，进入ANSYS主菜单Solution->Analysis Type->Analysis Options，设定求解一阶稳定荷载。

(19) 再次求解操作，进入ANSYS主菜单Solution->Solve->Current LS。
(20) 这时进入后处理操作，进入ANSYS主菜单General Postproc->Read results->Last Set，读入最后一步结果，选择ANSYS主菜单General Postproc->Plot Results->Deformed Shape，就可以得到结构的失稳形式以及相应的失稳荷载放大倍率，即118.602，如图所示。

(21) 进入ANSYS顶部的Parameters->Get scalar data，在Get Scalar Data中选择Results data，在结果中选择Model results (模态结果)，如图所示。点击OK进入下一个窗口Get Model Results，选择将模态结果存入名称为Freq1的变量。模态为第一阶模态。

(22) 这时从ANSYS顶部窗口Parameters->Scalar Parameters中就可以看到提取出来的一阶频率为118.6021，如图所示

(23) 以上过程为分析特征值失稳的一般过程。但是，ANSYS在分析特征值失稳的过程中有一个缺陷。众所周知，所谓特征值失稳计算就是用结构的材料刚度矩阵减去荷载作用下结构的几何刚度乘以一个系数，当总刚度矩阵奇异时的就是失稳特征值。ANSYS在处理荷载引起的刚度矩阵时不能区分我们需要分析的外力荷载（在本算例中是顶部集中力）和不需要的结构内力（例如本算例中预应力）对几何刚度矩阵的贡献。因此，得到的特征值屈服荷载118602N（等于顶部荷载的初始值100×特征值失稳的一阶频率118.602）是不正确的。因此，必须通过以下的迭代计算来解决该问题。
(24) 迭代计算的基本思想是根据第一次算出来荷载的放大倍数，调整加在结构上的外荷载，再求解新荷载下的特征值失稳放大倍数。重复上述操作直至特征值失稳的放大倍数基本等于1。这时加在结构上的外荷载就是真正的特征值失稳荷载。而且不会和内力结果发生矛盾。迭代计算的命令流如下：
! 设定最多迭代100次
*DO,I,1,100
FINISH
/SOLU
! 给结构施加新的荷载
FK,2,FZ,-FORCE
! 进行静力分析
ANTYPE,0
!设定时间步
TIME,1
AUTOTS,0
NSUBST,1, , ,1
SSTIF,ON
SOLVE
FINISH
! 进行特征值失稳分析
/SOLU
ANTYPE,BUCKLE ! Buckling analysis
BUCOPT,LANB,1 ! Use Block Lanczos solution method, extract 1 mode
MXPAND,1 ! Expand 1 mode shape
PSTRES,ON ! INCLUDE PRESTRESS EFFECTS
SOLVE
FINISH
! 得到当前的特征值失稳一阶频率（放大倍率）
*GET,FREQ1,MODE,1,FREQ
*IF,ABS(FREQ1-1),LT,0.01,THEN !如果频率误差小于1%，则退出循环
*EXIT
*ENDIF
FORCE=FORCE*FREQ1 ! 否则，将荷载乘以新的放大倍率再次计算
*ENDDO

(25) 重复上述过程后，得到结构实际的特征值失稳荷载为198174N，可见如果不排除预应力等内力对几何刚度矩阵的影响，计算出来的特征值失稳荷载(118602N)要偏小很多。
(26) 最后我们来进行非线性屈曲分析。特征值失稳计算得到的失稳荷载事实上是对于理想构件理想材料的上限解。实际结构由于种种初始缺陷或者材料非线性的影响，其失稳荷载往往要比特征值失稳荷载要小。因此有必要进行非线性失稳分析。而在非线性失稳分析中，需要引入初始缺陷。初始缺陷的选择方式很多，比较常用的是将结构的特征值失稳形状作为最不利的初始缺陷加在结构上。ANSYS中提供了UPGEOM命令，可以很方便的实现如上过程，下面详细介绍如下
(27) 首先我们要得到特征值失稳模态中结构的最大变形是多少，首先在ANSYS主菜单General Postpro->List Results->Sorted Listing->Sort Nodes，选择对总位移(USUM)的结果，对所有节点的最大位移进行排序。

(28) 接着，在ANSYS顶部菜单中选择Parameters->Get Scalar Data，在Get Scalar Data中选择Results data->Other operations

(29) 在Get Data from Other POST1 Operations窗口中选择从刚才的排序结果中获取数据。需要的数据为刚才排序的最大值，将该结果放到DMAX变量里面。

(30) 这时从Scalar Parameter窗口中就可以看到DMAX=1
(31) 下面用UPGEOM命令来调整结构形状，施加初始缺陷。在ANSYS主菜单中选择Preprocessor->Modeling->Update Geom。UPGEOM将从结果文件中读取结构变形并调整结构几何形状。我们设定放大倍数为OFFSET/DMAX，即我们需要的最大初始缺陷为OFFSET=0.1，而结果文件中的最大变形为DMAX，所以放大倍率为OFFSET/DMAX。最后指定特征值失稳分析的结果文件名称为Case03.rst。

(32) 更新完结构的几何形状以后，就可以进行非线性求解了。我们已知结构出现特征值失稳的荷载为198174N，我们设定结构受到的最大荷载为三倍的特征值失稳荷载，即在ANSYS窗口顶部的输入栏输入FORCE=FORCE*3.。将该荷载加载结构上，同时设置分析类型为静力分析。
(33) 最后我们要在求解器控制里面设定求解方法为弧长法，这样能够自动跟踪失稳路径。进入ANSYS主菜单Solution->Analysis Type->Analysis Options，首先在基本选项中设置分析类型为大位移分析，考虑预应力。设定分析子步数为20步，并每步都输出结果。

(34) 在Solution Controls窗口里面选择Advanced NL页面，点击激活弧长法(Arc-length Method)，本算例中弧长法参数可以采用默认值。

(35) 进入ANSYS主菜单Solution->Solve->Current LS，开始计算当前问题，由于本问题是我们接触到的第一个非线性问题，所以有必要介绍一下ANSYS非线性计算窗口中输出结果的含义。我们知道，非线性计算都有一个收敛的问题，在本算例中，ANSYS在后处理窗口中输出4个计算中间变量：不平衡力的2范数(F L2)，不平衡力的收敛容差(F CRIT)，不平衡弯矩的2范数(M L2)和不平衡弯矩的收敛容差(M CRIT)。窗口中坐标轴水平方向为时间（计算进程），纵坐标为相应的数值。以力为例，如果不平衡力的2范数(F L2)高于不平衡力的收敛容差(F CRIT)，说明没有收敛，要继续计算；如果(F L2) 小于(F CRIT)，在图中表现为(F L2)和(F CRIT)两条曲线相交，则说明该步计算已经收敛，可以计算下一个荷载步。从本算例看收敛的情况还是很好的，一般迭代一两次就收敛了。

(36) 计算完成后，进入ANSYS的时程后处理器，选择TimeHist Postpro，在Time History Variables中点击工具栏中添加变量按钮，选择添加节点Z方向变形，如图所示。选择节点2。

(37) 最后在Time History Variables窗口中设定以UZ_2为X坐标，将Time 点亮作为Y轴，点击工具栏上的绘图按钮，就可以得到如图所示的相对荷载－顶点变形曲线。

(38) 在通用后处理器General Postproc中，选择ANSYS主菜单General Postproc->Read Results->By Time/Freq，设定显示时间为0.2(即相当于1/5的最大荷载)，选择ANSYS主菜单General Postproc->Plot Results->Deformed Shape。绘制出结构当前变形如图。将结果时间逐步增大，显示Time=0.4、0.6步时的变形如图。