ANSYS梁单元截面残余应力施加方法及行为影响研究-峰设教育

在钢结构中，结构构件常用到轧制或者焊接工艺，这些工艺会使构件温度场产生不均匀变化，从而使构件中产生复杂的残余应力分布。
残余应力是一种自相平衡的力系，当构件承受荷载时，如受拉、受压等，荷载引起的应力将与截面残余应力相叠加，从而使构件某些部位提前达到屈服强度，并发生塑性变形，故会严重降低构件的刚度和稳定性以及结构疲劳强度。
本文以工字钢截面为例，系统阐述了梁单元截面残余应力的施加方法及行为影响。所用工字钢截面尺寸及残余应力分布如下所示（尺寸单位为mm，应力单位为Mpa）：
（为扩大残余应力的影响，这里残余应力取值较大，具体代码中已经进行了参数化，可随意更改残余应力大小从而研究其影响）

一、残余应力的施加
梁单元残余应力的施加目前只有通过命令Inistate施加，无GUI对应的操作，WB界面也只能通过插入命令流来解决。翻阅Inistate命令的帮助，如下所示：

针对梁单元截面而言，命令的参数解释如下：
ELID:单元的编号
EINT:单元积分点编号，对于Beam188而言，这里指I节点和J节点，也即一个单元有两个积分点。
KLAYER: 截面积分点所在的栅格编号。
PARMINT: 截面积分点的编号。
CXX~CXZ:具体的应力值。
从如上参数可见，要想对梁截面施加残余应力，首先应该了解截面积分点的组成。
在下面的这篇文章中，水哥已经详细的介绍了什么是梁单元截面的栅格和积分点，有兴趣的同学可以翻阅下，此处不在详细阐述。
如何获取梁单元截面栅格点和积分点计算结果
1、梁单元积分点位置坐标的获取
通过对命令的理解，残余应力的施加实际上就是对截面不同位置积分点初应力的施加。而这个初应力的大小又和截面的残余应力分布有关，也即可通过积分点的位置坐标来计算具体积分点位残余应力的大小。
所以在施加残余应力之前，第一步便是要获取截面积分点的位置坐标。
目前，水哥暂未发现在前处理中可以直接获取截面积分点坐标的命令，但是可通过列表截面详细来获取，列表命令Slist。
如下命令流所示：
sectype,1,beam,I
secdata,120,120,240,8,8,6
slist,1,,,full
上述命令流根据前面的截面参数，定义了一个工字钢截面，然后用slist命令显示出了截面的详细信息。

弹出来的界面里面，包括了四个方面的内容：
1）、截面属性，包括面积、重心、惯性距等信息；
2）、截面栅格点的位置坐标，也即Beam Section Cells Nodal Coodrinates
3）、截面栅格点位组成，也即Beam Section Cell Connectivity
4）、截面积分点位置坐标，也即Beam Section Cell Integration Point Coordinates
而我们需要的就是提取积分点位置坐标供下一步具体残余应力值的计算。

点击File-Save as，将该文件存为一个后缀名为lis的文件，然后用Matlab或者Python解析该文件，得到积分点坐标文件，本文采用Python来处理。
2、积分点残余应力的计算
积分点位置坐标文件得到后，进入ANSYS，然后读取坐标文件，获取积分点的个数，然后定义一个数组用以存储每个积分点的残余应力。
根据前文残余应力的分布，计算每个点位的残余应力，代码如下所示：

!读入积分点坐标
*create,AAA,mac
*dim,locint,array,32,2
*vread,locint(1,1),D:/ANSYS/locint,txt,,jik,2,32
(2F15.3)
*end
AAA
!根据坐标点计算残余应力
!定义数组存储每个积分点的残余应力
*dim,Rev_int,array,32,1
!首先确定截面原点位置，然后在根据Z坐标来确定是腹板还是翼缘，从而确定采用的公式
RevMaxI=23.4
RevMaxJ=-23.4*2.4
*do,i,1,32
*if,locint(i,2),ge,-112,then
*if,locint(i,2),lt,0,then
Rev_int(i)=RevMaxI
*endif
*endif
*if,locint(i,2),lt,112,then
*if,locint(i,2),ge,0,then
Rev_int(i)=RevMaxI
*endif
*endif
!上翼缘
*if,locint(i,2),gt,112,then
Rev_int(i)=RevMaxI+(RevmaxJ-RevmaxI)/60*abs(locint(i,1))
*endif
!下翼缘
*if,locint(i,2),lt,-112,then
Rev_int(i)=RevMaxI+(RevmaxJ-RevmaxI)/60*abs(locint(i,1))
*endif
*enddo

3、施加残余应力
在施加的过程中，应注意积分点的排序是根据栅格的次序来定的，每个栅格对应四个积分点，具体代码如下所示：

*get,elenum,elem,0,count
inistate,set,-2,stre
*do,i,1,elenum
*do,j,1,8
inistate,define,i,1,j,1,Rev_int(4*(j-1)+1)
inistate,define,i,1,j,2,Rev_int(4*(j-1)+2)
inistate,define,i,1,j,3,Rev_int(4*(j-1)+3)
inistate,define,i,1,j,4,Rev_int(4*(j-1)+4)
inistate,define,i,2,j,1,Rev_int(4*(j-1)+1)
inistate,define,i,2,j,2,Rev_int(4*(j-1)+2)
inistate,define,i,2,j,3,Rev_int(4*(j-1)+3)
inistate,define,i,2,j,4,Rev_int(4*(j-1)+4)
*enddo
*enddo

4、残余应力的复核
施加残余应力后，要对施加情况进行及时的复核，方法为不加任何外部荷载，求解结构，然后打开单元形状，绘制X方向轴向应力查看,且注意，由于查看截面应力的时候是查看的截面栅格点的应力，有的时候可能应力要比施加的应力要大，这主要是积分点应力结果外推的影响，可以用eresx命令关闭外推，复制积分点的应力结果，这样更加准确。
本次计算的结构断臂柱，长度为1.6m，研究其对弱轴的影响，故强轴方向平动均约束。
求解之后的截面应力云图如下所示:

从云图中可见，截面正应力的值和设计值一致，但是截面压应力的最大值与设计值具有一定的差异，造成这个问题的主要原因在于默认的截面网格划分较为粗糙，如下所示的默认栅格布置，积分点位置在箭头所致处，而这离端部还具有一定的距离。

如果为改善残余应力施加效果，可通过自定义截面的方式人工加密截面，如这里手动对网格进行加密，加密后的截面栅格如下所示：

按上述的操作步骤，重新进行残余应力的计算、施加并复核，得到的截面应力云图如下所示：

从图中可见，两者应力已经非常接近。但加密截面后也会对加大模型的计算量，实际运用中如何考虑还需综合评定。

二、行为影响研究

为研究残余应力对构件受力的影响，本文以长度为1.6m的悬臂柱为例，柱材料为Q235,材料模型理想弹塑性本构模型BKIN,分别从静力分析、模态分析、特征值分析、稳定性分析四个方面来进行比较,其中稳定性分析考虑L/500的初始缺陷。

此处分别考虑三个模型。

模型一：未考虑残余应力

模型二：考虑残余应力并采用默认网格划分

模型三：考虑残余应力并手动加密

分析的力学模型如下所示：

静力分析：