平面应力脆性断裂相场AT2模型

1 UEL用法

使用UEL子程序进行计算时，首先通过Abaqus建模生成计算所需的inp文件，然后需要对Abaqus的inp文件进行如下几处的修改，以附件中test\single_edge_notched_tension\length0.01文件夹下的SEN_plane_stress_uel.inp文件为例：

（1） 首先添加UEL的定义

值得说明的是，方框中的定义方式能够使得传入UEL的节点自由度列阵为：

下标代表节点编号，即排列顺序为首先排列所有节点的位移自由度，然后是所有节点的相场自由度。

当然，也可以使用如下定义方式：

此时，传入UEL的节点自由度列阵为：

这一定义方式使得在子程序中需要额外对位移场和相场的编号顺序进行处理（当然也并不麻烦），在这里为简便起见，选择第一种定义方式。

（2） 然后紧跟着定义UEL的单元集合

（3）添加叠加的用于可视化的Abaqus标准单元，该单元使用UMAT材料，通过UMAT接口中的状态变量进行可视化

框内的文件使用python脚本python_generate_overlaying_element_2D.py进行生成，可运行批处理文件run_visualization.bat，具体命令为：

python python_generate_overlaying_element_2D.py SEN_plane_stress_uel.inp 4521

需提供两个参数，分别为inp文件名和单元总个数，可以根据具体job进行修改。

（4）添加UEL和可视化UMAT单元的性质

其中UEL的单元性质分别是杨氏模量、泊松比、断裂韧性、相场特征宽度值、保证数值稳定性的小值、平面应力问题中的厚度值

UMAT的材料性质为杨氏模量、泊松比和单元总个数，其中杨氏模量设置为一个极小的值，不同job需要修改单元总个数的值。状态变量的个数设置为8.

（5）修改分析步的设置

具体数值可以酌情修改，每个变量的含义可以查找Abaqus文档。

（6）添加状态变量的场输出，用于可视化

2 理论

将系统的总势能表示为如下两项：

式中第一项能量为：

考虑损伤带来的退化，弹性能的表达式为：

式中

k为一个小值，用于防止数值不稳定现象。另一项断裂能为：

因此代入具体表达式可将系统总势能表达为：

对上述能量进行一阶变分可得：

即可得弱形式方程为：

具体外力虚功为：

式中本构方程为：

该弱形式方程是后续推导有限元方程的基础。同时，通过弱形式方程也可推导得到强形式的控制方程，即位移场和相场的控制方程。对上述弱形式进行分部积分可得：

因次位移场和相场的强形式控制方程为：

以及相应的边界条件为：

3 有限元离散

为推导有限元离散方程，对位移场和相场控制方程的弱形式进行处理：

对位移场和相场进行插值可得：

m指单元节点的个数。因此相应的梯度场可以插值为：

B矩阵的是由形函数对物理坐标的导数组成的。同理有：

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

使用牛顿迭代法求解上述非线性系统。更新格式为:

刚度矩阵为：

为了保证损伤不能愈合，即：

需要做出一些修改，即取历史上最大的弹性应变能，即：

4 代码

代码在附件中的src文件夹下，包含有函数支持文件AT2_plane_stress_uel_pack.f90和主函数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进行结果可视化
    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 !对高斯积分点进行循环
        DO K1=1,NSVINT  !对每个高斯积分点上的状态变量进行循环
            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 单边缺口板拉伸实验

对于单边缺口板的拉伸案例，我们选取四种相场特征裂纹宽度的情况，来展示特征宽度对于结果的影响。

第一个取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.