基于用户子程序的自定义应力场的施加

在材料的成型加工中或者表面处理工艺中，由于材料局部产生了严重的塑性变形，在加工结束后通常会在工件内部留下残余应力，如通过喷丸工艺可以在材料内部引入残余压应力场，提高材料的疲劳寿命。在通过ABAQUS进行有限元分析时，通常应力场可以通过初始条件的方式施加到有限元模型中，但在ABAQUS CAE的界面中仅提供了一些简单的应力场的施加方式，如通过直接指定模型中某一区域的应力分量的取值或者指定通过其他分析步计算得到的应力场结果作为残余应力场。

那么，在ABAQUS中是否可以在模型中引入自定义的应力场呢，答案是肯定的。对于用户自定义应力场的施加，可以通过ABAQUS用户子程序SIGINI来实现。

1. 用户子程序SIGINI简介

用户子程序SIGINI可以用来定义初始应力场。该子程序将在分析开始时每一个可被调用的特定材料点上被调用。该用户子程序可以将材料点上所有被激活的应力分量定义为关于坐标，单元号，积分点号等变量的函数。

1.1 应力分量

在该子程序中被定义的应力分量的数量取决于当前被调用单元的单元类型，相关单元类型对应的应力分量可参考帮助文献。该子程序中应力分量定义的顺序必须与参考文献中给定的定义顺序保持一致。例如对于三维的连续单元，六个应力分量的定义顺序必须为：σ11，σ22，σ33，σ12，σ13，σ23。

1.2 初始应力场平衡

在施加初始应力场用于后续计算时，你必须保证通过用户子程序定义的初始应力场与施加的载荷或分布载荷是平衡的，该平衡可以通过使用静力分析步来检查。

1.3 用户子程序界面

     SUBROUTINE SIGINI(SIGMA,COORDS,NTENS,NCRDS,NOEL,NPT,LAYER,     1 KSPT,LREBAR,NAMES)C      INCLUDE 'ABA_PARAM.INC'C      DIMENSION SIGMA(NTENS),COORDS(NCRDS)      CHARACTER NAMES(2)*80

      用户定义变量SIGMA(NTENS)的程序

      RETURN  END

1.4 需定义的变量

SIGMA(1)

第一个应力分量

SIGMA(2)

第二个应力分量

Etc.

可定义的应力分量数量共NTENS个，NTENS取决于当前被调用材料点所属单元的单元类型。

1.5 常见用于传递信息的变量

COORDS

包含当前被调用材料点初始坐标的数组

NTENS

       需要定义的应力分量的数量，取决于单元类型

NCRDS

       坐标的个数

NOEL

       单元号

NPT

       当前被调用单元的积分点号

更多该用户子程序用于传递信息的变量可参考帮助文献。

2. 焊接钢管的残余应力场定义

为了说明该用户子程序的使用方法，下面构建一个简单的计算实例。现假设有一焊接钢管，其焊接区域由于热应力的作用产生了一残余应力场。为了计算简便，现假设该残余应力场仅在钢管的轴向方向上存在正应力分量，并且在钢管内径表现为残余压应力，在钢管外径表现为残余拉应力，残余应力沿着壁厚方向线性分布，如图2.1所示。

图2.1 模型几何及初始应力场定义

图2.1中模型的全局坐标取在钢管的几何中心处，这样的坐标建立方式将有利于在用户子程序中计算相应的坐标参数。圆形虚线框对应模型中残余应力场的作用区域，残余应力在该区域沿壁厚方向线性分布，其分布规律如图2.1中右侧所示。

现假设残余应力场在管壁轴向方向的作用范围H为18mm，管道外径上作用的最大残余应力σ_max为100MPa，管道内径上作用的最小残余应力σ_min为-100MPa，管道内径R_min为30mm，管道外径（在残余应力场区域包含了坡口的厚度）R_max为37.8mm，则残余应力场内任意区域的残余应力可表示为：

其中R为当前位置距管道轴线的径向距离，可表示为：

在ABAQUS用户子程序中，当前被调用的材料点的坐标被储存在变量COORDS中，COORDS(1)代表坐标x，COORDS(2)代表坐标y，COORDS(3)代表坐标z。

焊接钢管的有限元模型如图2.2所示，其中焊缝区域采用细化网格进行过渡。

图2.2 有限元模型

该有限元模型的分析过程在ABAQUS中完成，分析步采用通用静力分析步，其它步骤不再赘述。

3. 用户子程序实现

完整的用户子程序如下。

SUBROUTINE SIGINI(SIGMA,COORDS,NTENS,NCRDS,NOEL,NPT,LAYER,     1 KSPT,LREBAR,NAMES)C      INCLUDE 'ABA_PARAM.INC'C      DIMENSION SIGMA(NTENS),COORDS(NCRDS)      CHARACTER NAMES(2)*80C     程序用于施加初始管道焊接残余应力C     采用单位制(N,mm,MPa)      REAL,PARAMETER::STRESS_AXIAL_DIST=18         !残余应力在管道轴向分布的长度      REAL,PARAMETER::TUBE_INTERNAL_RADIUS=30      !管道内径      REAL,PARAMETER::TUBE_EXTERNAL_RADIUS=37.38   !管道外径(包含了坡口的径向长度)            REAL,PARAMETER::STRESS_MAX=100               !假设最大残余应力出现在管道外径      REAL,PARAMETER::STRESS_MIN=-100              !假设最小残余应力出现在管道内径            REAL::RADIUS                                 !当前积分点到管道中心的径向距离C     假设残余应力将沿着管道厚度方向均匀分布,从最小残余应力线性过渡到最大残余应力      C     计算管道残余应力分布      IF(ABS(COORDS(3)).LE.STRESS_AXIAL_DIST/2)THEN   !当前材料积分点位于残余应力作用的轴向区域          RADIUS=SQRT(COORDS(1)**2+COORDS(2)**2)      !计算当前材料积分点距离管道中心的径向距离C         计算残余应力分量的取值,假设管道只在轴向方向(Z向)存在残余应力          SIGMA(3)=(STRESS_MAX-STRESS_MIN)/(TUBE_EXTERNAL_RADIUS-     &    TUBE_INTERNAL_RADIUS)*(RADIUS-TUBE_INTERNAL_RADIUS)+STRESS_MIN      END IF      
      RETURN      END

为了ABAQUS在计算过程中能够调用用户子程序SIGINI，需要对inp文件进行修改。点击Model→Edit Keywords修改当前模型的关键字，在分析步前面添加关键字*Initial Conditions，TYPE=STRESS，USER，如图3.1所示。

图3.1 修改关键字

4. 计算结果

最终的计算结果如图4.1所示。

图4.1 有限元计算结果

从图4.1中可以看出由于我们人为定义的残余应力场没有实现应力的平衡，因此当我们使用通用静力分析步进行计算之后将得到平衡的残余应力场。我们提取了轴向应力S33沿壁厚方向的分布，如图4.2所示。

图4.2 残余应力沿壁厚方向分布

从图4.2可以看出由于在初始时刻施加的应力场没有实现应力平衡，因此最终得到的残余应力场与初始施加的应力场存在很大的区别。但由于该方法具有通用性，我们可以很方便地利用该方法将其他软件计算得到的应力场施加到ABAQUS的有限元模型中，因此仍不失为一种处理自定义应力场的最佳方法。