找回密码
 加入慢享
猜你喜欢
旅行常客论坛

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

[复制链接]
发表于 2022-6-26 22:33:26 | 显示全部楼层 |阅读模式
线弹性断裂力学(linear elastic fracture mechanics, LFEM)通常假设裂纹尖端处于小范围屈服状态(裂纹尖端的塑性屈服区域远小于应力强度因子K主导的区域),此时将应力强度因子K作为描述裂纹尖端奇异性的参量是有效的。但是在许多情况下,如果裂纹尖端存在大范围屈服,则此时运用线弹性断裂力学是不合适的。

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

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

1. J积分概念

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

其中w为应变能密度,Ti为张力矢量的分量,ui为位移矢量的分量,ds为沿着积分路径的长度增量。应变能密度w可以定义为:

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

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

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

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

对于平面应力,有:

对于平面应变,有:

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

2. J积分的数值计算

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

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

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

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

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

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

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

其中n为该单元具有的节点数,qI为节点处的q函数的取值,NI为单元形函数。

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

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

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

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

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

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

2.2 4节点平面等参单元

2.2xy为单元的全局坐标,st为单元的等参坐标。采用等参单元在进行高斯积分时将非常方便。

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

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

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

上面各式中,N为形函数矢量,xy为单元节点处的坐标矢量,uv为单元节点处的位移矢量,q为单元节点处的q值构成的矢量,分别可表示为:

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

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

其中xcyc为单元中心处的全局坐标。

下面假设全局坐标xy可以通过如下形式用等参坐标st表示出来:

带入节点坐标x1,x2,x3,x4,y1,y2,y3,y4以及对应的st的取值,可以将aist来表示,最终得到:

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

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

因此有:

Jacobian矩阵为:

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

故应变可表示为:

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

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

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

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

上式中:

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

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

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

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

下面采用上述流程计算中心裂纹板裂尖的J积分,中心裂纹板的模型与“基于外推法的应力强度因子的计算”中采用的模型完全相同。为了便于计算,将积分区域的内边界Γ0选取为裂纹尖端,积分区域内的单元及每个单元节点对应的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可以提取得到该单元四个节点的坐标xy,以及节点位移uv的取值,而各节点对应的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∙mm1/2,计算误差为4.01%。相比于通过外推法计算得到的应力强度因子,采用J积分计算得到的结果更接近解析解。

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

作者:邓一超

精选:王华军

编辑:刘义美

爱仿真 关注仿真科技控视频号!

回复

使用道具 举报

快速回复 返回顶部 返回列表