1 UEL用法

使用UEL子程序進行計算時，首先通過Abaqus建模生成計算所需的inp文件，然后需要對Abaqus的inp文件進行如下幾處的修改，以附件中test\single_edge_notched_tension\length0.01文件夾下的SEN_plane_stress_uel.inp文件為例：

（1） 首先添加UEL的定義

值得說明的是，方框中的定義方式能夠使得傳入UEL的節(jié)點自由度列陣為：

下標代表節(jié)點編號，即排列順序為首先排列所有節(jié)點的位移自由度，然后是所有節(jié)點的相場自由度。

當然，也可以使用如下定義方式：

此時，傳入UEL的節(jié)點自由度列陣為：

這一定義方式使得在子程序中需要額外對位移場和相場的編號順序進行處理（當然也并不麻煩），在這里為簡便起見，選擇第一種定義方式。

（2） 然后緊跟著定義UEL的單元集合

（3）添加疊加的用于可視化的Abaqus標準單元，該單元使用UMAT材料，通過UMAT接口中的狀態(tài)變量進行可視化

框內(nèi)的文件使用python腳本python_generate_overlaying_element_2D.py進行生成，可運行批處理文件run_visualization.bat，具體命令為：

python python_generate_overlaying_element_2D.py SEN_plane_stress_uel.inp 4521

需提供兩個參數(shù)，分別為inp文件名和單元總個數(shù)，可以根據(jù)具體job進行修改。

（4）添加UEL和可視化UMAT單元的性質(zhì)

其中UEL的單元性質(zhì)分別是楊氏模量、泊松比、斷裂韌性、相場特征寬度值、保證數(shù)值穩(wěn)定性的小值、平面應力問題中的厚度值

UMAT的材料性質(zhì)為楊氏模量、泊松比和單元總個數(shù)，其中楊氏模量設置為一個極小的值，不同job需要修改單元總個數(shù)的值。狀態(tài)變量的個數(shù)設置為8.

（5）修改分析步的設置

具體數(shù)值可以酌情修改，每個變量的含義可以查找Abaqus文檔。

（6）添加狀態(tài)變量的場輸出，用于可視化

2 理論

將系統(tǒng)的總勢能表示為如下兩項：

式中第一項能量為：

考慮損傷帶來的退化，彈性能的表達式為：

式中

k為一個小值，用于防止數(shù)值不穩(wěn)定現(xiàn)象。另一項斷裂能為：

因此代入具體表達式可將系統(tǒng)總勢能表達為：

對上述能量進行一階變分可得：

即可得弱形式方程為：

具體外力虛功為：

式中本構(gòu)方程為：

該弱形式方程是后續(xù)推導有限元方程的基礎(chǔ)。同時，通過弱形式方程也可推導得到強形式的控制方程，即位移場和相場的控制方程。對上述弱形式進行分部積分可得：

因次位移場和相場的強形式控制方程為：

以及相應的邊界條件為：

3 有限元離散

為推導有限元離散方程，對位移場和相場控制方程的弱形式進行處理：

對位移場和相場進行插值可得：

m指單元節(jié)點的個數(shù)。因此相應的梯度場可以插值為：

B矩陣的是由形函數(shù)對物理坐標的導數(shù)組成的。同理有：

代入到弱形式方程中可得殘值方程；

使用牛頓迭代法求解上述非線性系統(tǒng)。更新格式為:

剛度矩陣為：

為了保證損傷不能愈合，即：

需要做出一些修改，即取歷史上最大的彈性應變能，即：

4 代碼

代碼在附件中的src文件夾下，包含有函數(shù)支持文件AT2_plane_stress_uel_pack.f90和主函數(shù)AT2_plane_stress_uel_main.f90，可運行run_obj.bat生成相應的obj文件，即AT2_plane_stress_uel_main-std.obj。

UEL需要更新單元剛度矩陣和單元殘值，具體公式在上述理論部分已詳細給出。主程序代碼如下：

include "AT2_plane_stress_uel_pack.f90" 

subroutine uel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,                  &
                props,nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,dtime,   &
                kstep,kinc,jelem,params,ndload,jdltyp,adlmag,predef,npredf, &
                lflags,mlvarx,ddlmag,mdload,pnewdt,jprops,njprop,period)

    use AT2_plane_stress_uel_pack

    include 'aba_param.inc'

    dimension   rhs(mlvarx,*),amatrx(ndofel,ndofel),props(*),               &
                svars(*),energy(8),coords(mcrd,nnode),u(ndofel),            &
                du(mlvarx,*),v(ndofel),a(ndofel),time(2),params(*),         &
                jdltyp(mdload,*),adlmag(mdload,*),ddlmag(mdload,*),         &
                predef(2,npredf,nnode),lflags(*),jprops(*)
    
    ! **************************************************************
    ! 使用umat進行結(jié)果可視化
    REAL*8 UVAR(70000,NSVINT, NGAUSS)
    INTEGER JELEM,K1,K2,KINTK
    COMMON/KUSER/UVAR
    ! **************************************************************
    
    real(8) :: E
    real(8) :: nu
    real(8) :: Gc
    real(8) :: L0
    real(8) :: k
    real(8) :: thick ! 平面應力單元厚度
    real(8) :: coords_plane(2,4)
    real(8) :: Ke(12,12)
    real(8) :: fint(12)

    E  = props(1)
    nu = props(2)
    Gc = props(3)
    L0 = props(4)
    k = props(5)
    thick = props(6)

    ! call uel
    coords_plane(1:2,1:4) = coords(1:2,1:4)
    call AT2_plane_stress_uel(E,nu,Gc,L0,k,thick, coords_plane, GAUSS_PT, u, du, Ke, fint, svars)
    
    rhs(:,1) = -fint
    amatrx = Ke

    ! **************************************************************
    ! for visualization
    DO KINTK = 1,NGAUSS !對高斯積分點進行循環(huán)
        DO K1=1,NSVINT  !對每個高斯積分點上的狀態(tài)變量進行循環(huán)
            UVAR(JELEM,K1,KINTK) = SVARS(NSVINT*(KINTK-1)+K1)
        END DO
    ENDDO
    ! **************************************************************

    return
end subroutine uel

! ******************************************************************
! umat進行可視化
subroutine umat(stress,statev,ddsdde,sse,spd,scd,                       &
                rpl,ddsddt,drplde,drpldt,                               &
                stran,dstran,time,dtime,temp,dtemp,predef,dpred,cmname, &
                ndi,nshr,ntens,nstatv,props,nprops,coords,drot,pnewdt,  &
                celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)

    use AT2_plane_stress_uel_pack

    include 'aba_param.inc'

    character*80 cmname
    dimension    stress(ntens),statev(nstatv),ddsdde(ntens,ntens),   &
                ddsddt(ntens),drplde(ntens),                        &
                stran(ntens),dstran(ntens),time(2),                 &
                predef(1),dpred(1),props(nprops),coords(3),         &
                drot(3,3),dfgrd0(3,3),dfgrd1(3,3)

    ! **************************************************************
    ! 可視化
    INTEGER NELEMAN, K1
    real(8) :: E,nu
    COMMON/KUSER/UVAR(70000,NSVINT, NGAUSS)

    E = props(1)
    nu = props(2)
    elem_num = nint(props(3))

    NELEMAN=NOEL-elem_num

    DO K1=1,NSVINT
        STATEV(K1)=UVAR(NELEMAN,K1,NPT)
    END DO

    ddsdde = AT2_plane_stress_uel_compute_De(E,nu)
    stress = matmul(ddsdde,stran + dstran)
    
    return
end subroutine umat

5 測試

5.1 單邊缺口板拉伸實驗

對于單邊缺口板的拉伸案例，我們選取四種相場特征裂紋寬度的情況，來展示特征寬度對于結(jié)果的影響。

第一個取l=0.1，相場分布圖如下：

第二個取l=0.05，相場分布圖如下：

第三個取l=0.025，相場分布圖如下：

第四個取l=0.01，相場分布圖如下：

5.2 單邊缺口板剪切實驗

單邊缺口板剪切實驗的相場分布圖如下：

6 參考文獻

[1]. 胡小飛, 張鵬, 姚偉岸. 斷裂相場法. 北京: 科學出版社; 2022.

[2] Kristensen PK, Martínez-Pa?eda E. Phase field fracture modelling using quasi-Newton methods and a new adaptive step scheme. Theoretical and applied fracture mechanics. 2020;107:102446.