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應力

不同時刻的應力云圖：

文中未就跟多細節一一展開，而是將主要問題和注意的地方提了出來。