ANSYS的生死单元模拟焊接过程

ANSYS的生死单元模拟焊接过程

1 概述

焊接模拟计算在CAE仿真是比较大的一块内容，也是比较复杂的一个过程，几个比较关键的问题是热源函数的描述、单元的融覆、热源的移动等等，通过单纯的GUI操作，无论使ANSYS还是Abaqus都不大可能完成这个过程，通常需要借助软件的内置语言。

本次主要介绍单元生死的应用，单元生死主要用于单元缺失的场合，比如凝固溶解过程，断裂过程，焊接过程等等，这些过程都是非线性或者时间历程过程，计算需要很多子步和迭代，为了在此过程中避免一遍一遍修改单元，便引入生死单元的概念，通俗的讲就是通过一些方法让单元失效，具体的改变是单元的弹性模量的改变，当单元死时，修改其弹性模量为非常小的值，让其在求解过程中不起作用。

详细地说，激活单元死这个状态时，ANSYS程序将单元刚度矩阵乘以很小的因子，程序默认值为1E-6，死单元的单元载荷为0，从而不对载荷向量生效，同样的，死单元的质量、阻尼、比热等等参数也设置为0，单元的应力应变也因此为0。

 

2 前处理

前处理包括单元定义、材料定义和建模，单元定义是需要注意单元属性，此次定义13号二维耦合单元，具有温度和位移自由度。

材料属性包括结构参数和热参数，具体包含弹性模量，泊松比，屈服强度，塑性属性，材料密度，热膨胀系数，热传导系数，比热容。焊接时温度较高，定义材料通常需要定义多个温度下的值。例如定义各材料在各温度点下的屈服应力和屈服后的弹性模量：

tb,bkin,1,5  

tbtemp,20,1  

tbdata,1,1200e6,0.193e11  

tbtemp,500,2  

tbdata,1, 933e6,0.150e11  

tbtemp,1000,3  

tbdata,1, 435e6,0.070e11  

tbtemp,1500,4  

tbdata,1, 70e6,0.010e11  

tbtemp,2000,5  

tbdata,1, 7e6,0.001e11

建立的二维模型如图1所示，中间三角区域为焊接区域。

图1 几何模型

划分网格后如图2所示：

图2 Mesh

3 单元排序

焊接区域的单元排序是一个必要的过程，只有将单元按照坐标顺序排列好，激活才好控制。

由于焊接时单元是从底部一层一层铺，因此需要将三角区域的单元按照Y坐标排序。单元的坐标是其形心处的坐标，如果一层的单元Y坐标一致，则将这一层再按照X坐标排序。

排序用ANSYS的APDL实现，通过Do循环一遍一遍的遍历单元。

其主要流程如图3所示：

图3 遍历过程

在进行这个过程之前，需要先选出焊接区域的所有单元，之后：

*get,nse,elem,,count       ！焊接区域单元数目nse

*dim,ne,,nse              ！定义数组，存储单元编号

*dim,nex,,nse             ！定义数组，存储坐标x

*dim,ney,,nse             ！定义数组，存储坐标y  

*dim,neorder,,nse         ！定义数组，存储按照坐标排序的单元标号

mine=0  

!**************************************************************

*do,i1,1,nse  

  esel,u,elem,,mine           ！提出排序了的单元

  *get,nse1,elem,,count        ！得到剩下的单元数目nse1

  ii=0  

  *do,i,1,emax  

    *if,esel(i),eq,1,then       ！表示单元是否被选中

      ii=ii+1  

      ne(ii)=i  

    *endif        

  *enddo                    ！此段循环用于得到剩下的单元的编号

  *do,i,1,nse1  

    *get,ney(i),elem,ne(i),cent,y  

    *get,nex(i),elem,ne(i),cent,x  

  *enddo                   ！此段循环得到剩下单元的坐标

  miny=1e20  

  minx=1e20                ！定义辅助常数，用于比较

  *do,i,1,nse1  

    *if,ney(i),lt,miny,then  

      miny=ney(i)  

      minx=nex(i)  

      mine=ne(i)  

    *else  

      *if,ney(i),eq,miny,then  

        *if,nex(i),lt,minx,then  

          miny=ney(i)  

          minx=nex(i)  

          mine=ne(i)  

        *endif  

      *endif  

    *endif  

  *enddo    ！此段循环用于将剩下的单元经过坐标比较，找到Y最小（X最小的单元编号）

  neorder(i1)=mine   ！得到这个编号的单元并存储

*enddo  

注意循环的时候会有图4的这个警告，这个没有关系，因为第一步剔除单元操作时实际并没有剔除成功，因为第一步循环时mine=0。

图4 Wanrning

 

4 加载求解

温度载荷，结构约束，无需多说。

焊接一个单元的时间假设为5s，起始时间为0，通过得到的单元数可以知道整个历程的时间。

求解之前先将所有焊接区域的单元杀死，通过循环遍历实现，命令是ekill，每次杀死一个单元：

*do,i,1,nse  

  ekill,neorder(i)  

  esel,s,live  

  eplot  

*enddo 

杀死后单元如图5所示：

图5 杀死单元

求解过程的APDL为：

*do,i,1,nse  

  ealive,neorder(i)  

  esel,s,live  

  eplot 

  esel,all  

  !******下面的求解用于建立温度的初始条件******

  t=t+dt1  

  time,t  

  nsubst,1  

  *do,j,1,4  

  d,nelem(neorder(i),j),temp,max_tem  

  *enddo  

  solve  

  !****下面的求解用于保证初始的升温速度为零****

  t=t+dt1  

  time,t  

  solve  

  !*********下面的步骤用于求解温度分布********** 

  *do,j,1,4  

   ddele,nelem(neorder(i),j),temp  

  *enddo  

  t=t+dt-2*dt1  

  time,t  

  nsubst,nsub1  

  solve  

*enddo  

t=t+50000   

求解完之后，继续求解冷却过程，冷却过程只需给一个时间就行。

5 计算结果

最后的应力分布如图6，因为在左端施加位移约束，最后由于焊接变形，右边发生翘曲。

图6 Mises应力

不同时刻的应力云图：

文中未就跟多细节一一展开，而是将主要问题和注意的地方提了出来。