断裂参数J积分的数值计算方法

线弹性断裂力学（linear elastic fracture mechanics, LFEM）通常假设裂纹尖端处于小范围屈服状态（裂纹尖端的塑性屈服区域远小于应力强度因子K主导的区域），此时将应力强度因子K作为描述裂纹尖端奇异性的参量是有效的。但是在许多情况下，如果裂纹尖端存在大范围屈服，则此时运用线弹性断裂力学是不合适的。

弹塑性断裂力学则可以应用到存在非线性行为（如弹塑性变形）的材料中。弹塑性断裂力学采用J积分来描述弹塑性材料裂纹尖端的状态，并且类比与应力强度因子K，同样可以作为裂纹扩展的计算参量。

下面讨论如何通过数值方法计算J积分，并间接计算得到应力强度因子。

1. J积分概念

Rice提出了一个不依赖于路径的围线积分用于裂纹的分析，随后他将该积分值命名为J积分。J积分本质上是一个非弹性裂纹体在裂纹扩展过程中的能量释放率。考虑任意一个围绕裂纹尖端的逆时针路径Γ，如图1.1所示，则J积分可以表示为：

其中w为应变能密度，T_i为张力矢量的分量，u_i为位移矢量的分量，ds为沿着积分路径的长度增量。应变能密度w可以定义为：

图1.1 围绕裂纹尖端的积分路径

其中σ_ij和ε_ij分别为应力和应变张量，而张力矢量T为积分路径上某一点上的应力矢量。如果我们建立一个被积分路径包围的自由体，则张力矢量T为作用在边界上的应力。张力矢量的分量可以表示为：

Rice证明了J积分的取值与选取的围绕裂纹尖端的积分路径无关，并且等于裂纹扩展中的能量释放率。特别的，对于线弹性材料，有：

其中K_I为I型断裂模式下的应力强度因子。

对于平面应力，有：

对于平面应变，有：

其中E为弹性模量，v为泊松比。

2. J积分的数值计算

上述关于J积分的计算表达式不利于数值计算的实现。Shih和Raju提出了等效区域积分法来计算J积分，上述J积分的表达式可以等效为对裂纹尖端附近的一个有限区域的面积分：

对于平面裂纹问题上式可以展开为：

对于二维问题，首先需要确定积分区域A。积分区域A的内边界Γ₀通常被选定为裂纹尖端，此时积分区域A即为由积分区域外边界Γ₁所包围的区域。积分区域外边界Γ₁通常与单元边重合。

函数q(x,y)只是便于积分表达式可以进行数值计算。但在积分区域的所有单元节点上，函数q(x,y)都必须有确定的值。Shih已经证明J积分的计算值对于假设的q函数的形式并不敏感，q函数的形式可以是任意的，但在积分区域的边界上q必须有正确的取值。例如，对于平面应力或平面应变问题，在Γ₀上有q=1，这通常对应于裂纹尖端，而在积分区域外边界Γ₁上有q=0。

图2.1给出了常用于二维问题的两种q函数的形式。一种为金字塔形，函数q(x,y)在裂纹尖端处取值为1，然后线性变化到外边界上，此时取值为0；另一种为平台形，q(x,y)只在外边界上取值为0，而在其他区域取值均为1。

图2.1 二维问题中两种常见的q函数形式

在某一个单元内任意位置处的q函数的取值可以通过单元形函数插值得到：

其中n为该单元具有的节点数，q_I为节点处的q函数的取值，N_I为单元形函数。

因此q函数的偏导可以表示为：

其中ξ_i为单元的等参坐标。

在求得q函数的偏导后，只需确定位移矢量u和v对坐标x的偏导，以及应力分量和应变能密度w的取值，则可利用高斯积分法来求解J积分。而这些变量在高斯点处的取值同样可以通过单元形函数插值得到。

2.1 平面单元的J积分计算方法

下面以平面应力单元为例来给出J积分的具体计算方法，对于平面应变单元，也可以采用类似的方法来进行计算。

二维平面问题的四节点的等参单元如图2.2所示。

图2.2 4节点平面等参单元

图2.2中x和y为单元的全局坐标，s和t为单元的等参坐标。采用等参单元在进行高斯积分时将非常方便。

单元中的变量均可通过单元节点上变量的取值利用形函数插值得到。考虑单元内的任意一点，其坐标可用节点坐标和形函数来表示：

该点处的位移分量同样采用节点位移和形函数来近似：

该点处对应的q函数的取值同样可表示为：

上面各式中，N为形函数矢量，x和y为单元节点处的坐标矢量，u和v为单元节点处的位移矢量，q为单元节点处的q值构成的矢量，分别可表示为：

上式中的单元节点处的变量均可以通过有限元软件直接计算得到，下面的关键是根据单元类型确定形函数。

从图2.2中4节点平面等参单元的定义可以看出，等参坐标系s-t位于单元的正中心，并且等参坐标系为自然坐标系，坐标轴s和t无需正交，也无需与全局坐标轴x和y平行。通过定义等参坐标取值为+1或者-1，可以完全限制等参单元的边线和角点。假设图2.2(b)中的畸形单元变为了如图2.2(a)中所示的规则单元，此时s-t坐标系可以与全局坐标x和y关联起来：

其中x_c和y_c为单元中心处的全局坐标。

下面假设全局坐标x和y可以通过如下形式用等参坐标s和t表示出来：

带入节点坐标x₁,x₂,x₃,x₄,y₁,y₂,y₃,y₄以及对应的s和t的取值，可以将a_i用s和t来表示，最终得到：

其中上式中的形函数分别为：

可以得到形函数对等参坐标的导数为：

因此有：

Jacobian矩阵为：

形函数对全局坐标的导数为：

故应变可表示为：

同样可以得到q函数的微分为：

在得到应力后，可以通过线弹性本构关系得到应力，对于平面应力问题有：

在求得应力和应变后，应变能密度w可表示为：

最后，在某个单元内的J积分可以表示为：

上式中：

该面积分可以采用高斯积分法进行求解：

式中的四个积分点坐标分别选取为：

重复上述过程计算积分路径围绕区域内的所有单元上的J积分，将结果累加即可得到裂纹尖端的J积分取值。

3. 中心裂纹板J积分计算

下面采用上述流程计算中心裂纹板裂尖的J积分，中心裂纹板的模型与“基于外推法的应力强度因子的计算”中采用的模型完全相同。为了便于计算，将积分区域的内边界Γ₀选取为裂纹尖端，积分区域内的单元及每个单元节点对应的q值如图3.1所示。

图3.1 积分区域单元及q值

图3.1中具有红色边线的单元即为覆盖在积分区域内的单元，黄色数字为单元节点号，蓝色数字为单元号。注意到由于模型具有对称性，因此只需评估图3.1中的两个单元再利用对称性即可计算J积分。

将上述过程编制成MATLAB程序，便于结果计算。

function Je = J_Integral(E,PR,u,v,x,y,q)%函数用于计算单元的J积分值%E为杨氏模量(MPa)%PR为泊松比%u,v为节点位移(mm)%x,y为节点坐标(mm)%q为节点处q函数取值
coord_s=[-1/sqrt(3) 1/sqrt(3) 1/sqrt(3) -1/sqrt(3)];   %积分点等参坐标scoord_t=[-1/sqrt(3) -1/sqrt(3) 1/sqrt(3) 1/sqrt(3)];   %积分点等参坐标tC_matrix=E/(1-PR^2)*[1 PR 0;PR 1 0;0 0 (1-PR)/2];  %本构矩阵
Je=0;    %单元J积分值初始化为0for ip=1:4      %积分点上累加函数I(s,t)取值    s=coord_s(ip);    t=coord_t(ip);     %读取积分点等参坐标        Ns=[-1/4*(1-t) 1/4*(1-t) 1/4*(1+t) -1/4*(1+t)];  %形函数对等参坐标s的偏导    Nt=[-1/4*(1-s) -1/4*(1+s) 1/4*(1+s) 1/4*(1-s)];  %形函数对等参坐标t的偏导        xs=Ns*x;xt=Nt*x;ys=Ns*y;yt=Nt*y;     %全局坐标对等参坐标的偏导        J_det=xs*yt-xt*ys;    %Jocbian矩阵行列式        Nx=(yt.*Ns-ys.*Nt)./J_det;     %形函数对坐标x的偏导    Ny=(-1*xt.*Ns+xs.*Nt)./J_det;  %形函数对坐标y的偏导        ux=Nx*u;      %u对坐标x的偏导    vx=Nx*v;      %v对坐标x的偏导        exx=Nx*u;            %x方向正应变    eyy=Ny*v;            %y方向正应变    exy=Nx*v+Ny*u;       %xy剪应变        qx=Nx*q;             %q函数对x的偏导    qy=Ny*q;             %q函数对y的偏导        strain=[exx;eyy;exy];       %应变    stress=C_matrix*strain;     %根据本构矩阵计算应力        sxx=stress(1);    %x方向正应力    syy=stress(2);    %y方向正应力    sxy=stress(3);    %xy剪应力        w=0.5*stress'*strain;     %计算应变能密度        Je=Je+((sxx*ux+sxy*vx-w)*qx+(sxy*ux+syy*vx)*qy)*J_det;  %累加函数I(s,t)取值  endend

利用有限元软件ABAQUS可以提取得到该单元四个节点的坐标x和y，以及节点位移u和v的取值，而各节点对应的q函数的取值已经在图3.1中标出，提取得到的单元16214的数据如表3.1所示,单元30000的数据如表3.2所示。

表3.1 单元16214节点输出结果

表3.2 单元30000节点输出结果

将上表中的数据输入到MATLAB中，如对于表3.2中的数据，输入：

E=200e3;PR=0.3;x=[20;20;21;21];y=[1;0;0;0];u=[-2.161e-3;-2.771e-3;-2.509e-3;-2.291e-3];v=[7.731e-4;0;0;3.759e-4];q=[0;1;0;0];J1=J_Integral(E,PR,u,v,x,y,q)

计算得到单元16214上的计算结果为0.144N/mm，单元30000上的计算结果为0.0165N/mm。由于模型具有对称性，故实际整个积分区域上的J积分取值为：

由此计算得到应力强度因子取值为：

由“基于外推法的应力强度因子的计算”中的解析公式可以计算得到应力强度因子的解析值为243.6MPa∙mm^1/2，计算误差为4.01%。相比于通过外推法计算得到的应力强度因子，采用J积分计算得到的结果更接近解析解。

实际上本文为了简便起见，选取的积分区域非常接近裂纹尖端，此时J积分的路径无关性将无法得到保证。如果选取远离裂纹尖端的积分区域，计算得到的结果将更接近于解析解。