基于权函数法的应力强度因子计算

10月前浏览432

1. 简介

在采用断裂力学进行疲劳评估时，应力强度因子是一个非常重要的计算参量。对于一些简单结构的应力强度因子，通常可以通过查阅应力强度因子手册中给出的经验公式计算得到，但计算手册给出的应力强度因子的计算公式非常有限，对于复杂应力场下的应力强度因子的计算，一般会采用有限元法来完成。虽然有限元可以用于计算任意应力场下的应力强度因子，但采用有限元计算应力强度因子时通常需要在采用奇异单元来构建裂纹网格，在进行前处理时非常繁琐且计算量相对较大。

Bueckner¹提出的权函数法可以用于计算任意应力分布下的应力强度因子，对于I型裂纹，其应力强度因子K可以表示为：

式（1）

其中a为裂纹长度，σ(x)为假想结构无裂纹体时沿裂纹面的应力分布，m(x,a)为权函数，权函数的选取只与裂纹体的几何有关，与载荷形式和取值无关。

从式（1）可以看出，由于权函数在计算时采用的是结构不存在裂纹体时的应力分布，因此该应力分布可以很容易地从有限元计算的结果中提取得到，并且在有限元分析时由于无需构造裂纹，因此极大的简化了前处理过程。

Glinka和Shen²通过推导得到，一维和二维裂纹的权函数可有如下统一形式：

式（2）

参数M₁,M₂和M₃只取决于裂纹体的几何形式，部分裂纹体的参数在已文献[3]给出。

2. 数值积分方法

从式（1）可以看出，在利用权函数求解应力强度因子时需要进行积分运算。由于应力分布σ(x)通常是通过有限元法计算得到的，由一系列离散的数据点构成的，因此只能通过数值方法来进行求解。下面给出求解权函数的数值积分方法。

代表应力分布σ(x)的一系列离散点可以通过多段线性的形式表示：

式（3）

式（4）

式（5）

其中σ_i(x)为位于区段(x_i,x_i₊₁)的代表真实应力分布的线性区段，k_i和b_i分别为该区段对应的斜率和常数，i的取值范围为1,…, N-1，N为离散数据点的个数。采用分段线性函数表示的应力分布如图1所示。

图1 分段线性函数表示的应力分布

由式（1）和式（2）可以得到，I型裂纹的应力强度因子K可表示为：

式（6）

将式（3）中采用分段线性表达的应力分布σ(x)带入式（6），并重新写为如下形式：

式（7）

其中K_IM₀, K_IM₁,K_IM₂, K_IM₃可表示为：

式（8）

式（9）

式（10）

式（11）

式（8）~（11）中的定积分可以直接通过解析方法计算得到：

则K_IM₀, K_IM₁,K_IM₂, K_IM₃的完整表达式为：

式（16）

式（17）

式（18）

式（19）

在已知M₁, M₂和M₃的情况下，首先通过式（4）和式（5）计算代表离散应力分布的多段线性函数的系数，随后通过式（7）和式（16）~（19）即可计算得到当前应力分布对应的应力强度因子。

3. 半无限宽平板中角裂纹的应力强度因子计算

下面采用权函数法对半无限宽平板中角裂纹的应力强度因子进行计算,假设裂纹面上作用的非线性应力分布σ(x)可表示为：

计算中取σ₀为1MPa，裂纹长度a为1mm，则裂纹面上的应力分布如图2所示。

图2 半无限宽平板上的应力分布

对于半无限宽平板中的角裂纹，基于Shen和Glinka²的推导，其权函数的系数M₁, M₂和M₃分别为：

下面采用第2节中给出的数值积分方法来计算半无限宽平板中角裂纹的应力强度因子，对应的Matlab程序如下所示。

  %程序用于计算半无限大平板中角裂纹的应力强度因子clear,clca=1;  %裂纹长度stress0=1;   %远场应力discrete_num=500;   %离散积分点数x=linspace(0,a,discrete_num);stress=stress0*(1+x./a+2*(x./a).^2+4*(x./a).^3+16*(x./a).^4+...    32*(x./a).^5);%采用多段线性函数离散化应力分布stress_discrete=zeros(discrete_num,1);   %存放离散化后的应力line_func_k=zeros(discrete_num-1,1);       %多段线性函数系数kline_func_b=zeros(discrete_num-1,1);       %多段线性函数系数bfor i=1:discrete_num-1     %计算多段线性函数系数    line_func_k(i)=(stress(i+1)-stress(i))/(x(i+1)-x(i));    line_func_b(i)=stress(i)-x(i)*line_func_k(i);end%验证多段线性函数正确性stress_discrete(1)=stress(1);for i=2:discrete_num    stress_discrete(i)=line_func_k(i-1)*x(i)+line_func_b(i-1);end%绘制应力分布及多段线性函数表示的应力分布figurescatter(x,stress)hold onplot(x,stress_discrete)title('应力分布及多段线性函数表示的应力分布');xlabel('沿裂纹长度方向距离/mm');ylabel('应力/MPa');legend('应力分布离散点','多段线性函数');%计算应力强度因子weight_func_M1=0.0719768;      %权函数系数weight_func_M2=0.246984;weight_func_M3=0.514465;KIM0=zeros(discrete_num-1,1);   KIM1=zeros(discrete_num-1,1);KIM2=zeros(discrete_num-1,1);KIM3=zeros(discrete_num-1,1);for i=1:discrete_num-1    KIM0(i)=2*sqrt(2)/(3*sqrt(pi))*((line_func_k(i)*x(i)+2*line_func_k(i)*a+...        3*line_func_b(i))*sqrt(a-x(i))-(line_func_k(i)*x(i+1)+...        2*line_func_k(i)*a+3*line_func_b(i))*sqrt(a-x(i+1)));         KIM1(i)=sqrt(2/(a*pi))*(line_func_k(i)/2*(x(i+1)^2-x(i)^2)+...        line_func_b(i)*(x(i+1)-x(i)));    KIM2(i)=2/(15*a)*sqrt(2/pi)*((3*line_func_k(i)*x(i)+2*line_func_k(i)*a+...        5*line_func_b(i))*(a-x(i))^(3/2)-...        (3*line_func_k(i)*x(i+1)+2*line_func_k(i)*a+...        5*line_func_b(i))*(a-x(i+1))^(3/2));    KIM3(i)=1/(a*sqrt(a))*sqrt(2/pi)*(line_func_k(i)/3*...        (x(i)^3-x(i+1)^3)+1/2*(line_func_k(i)*a-line_func_b(i))*...        x(i+1)^2-1/2*(line_func_k(i)*a-line_func_b(i))*x(i)^2+...        line_func_b(i)*a*(x(i+1)-x(i)));endKIM0_sum=sum(KIM0);KIM1_sum=sum(KIM1);KIM2_sum=sum(KIM2);KIM3_sum=sum(KIM3);K=KIM0_sum+KIM1_sum*weight_func_M1+KIM2_sum*weight_func_M2+...KIM3_sum*weight_func_M3;

由于应力分布是通过多段线性函数来表示的，因此在数值积分中计算得到的应力强度因子的精度取决于表示应力分布的离散点的个数，代表应力分布的离散点以及对应的多段线性函数如图3所示。

图3 代表应力分布的离散点及对应的多段线性函数

在程序中代表应力分布的离散点的个数由变量discrete_num控制，增大该变量的取值可以采用区段数更多的多段线性函数来逼近原始应力分布，从而提高应力强度因子的计算精度，但也会相应地增加计算时间。分别选取不同的离散点个数进行计算，得到的应力强度因子如图4所示。

图4 应力分布离散点与应力强度因子的关系

从图4中可以看出，随着离散点个数的增加，应力强度因子值迅速下降并趋于收敛，从计算时间和计算精度的角度来考虑，在程序计算时可以将离散点个数设置为100个，计算得到的应力强度因子K_I取值为：

而根据Stallybrass⁴解析计算得到的结果，应力强度因子K_I与裂纹长度a和分布压力σ₀之间的关系为：

式（23）

通过式（23）计算得到的应力强度因子K_I的取值为：

本文计算得到的应力强度因子与解析方法计算得到的应力强度因子误差仅为0.41%，因此通过本文给出的数值计算方法可以采用权函数法精确地计算应力强度因子。

4. 总结

为了计算如应力强度因子等断裂力学参数，通常会采用有限元法来进行应力分析。有限元法的优点为可以计算任意几何形状和任意载荷下裂纹体的应力强度因子，但其缺点为为了精确计算应力强度因子，需要在裂纹尖端采用非常精细的单元或者奇异单元来进行网格划分，增加了前处理的难度。在进行疲劳裂纹扩展分析时，裂纹长度不断变化，这需要不断对裂纹体的网格进行重划分，计算量极大。该过程的实现也需要非常高的建模技术和经验。

而权函数法克服了有限元法的上述缺点，只要裂纹体的形式是确定的，便可以计算任意载荷形式和任意裂纹长度下的应力强度因子，而任意裂纹体对应的权函数公式也可以借助于有限元法拟合得到。同时，权函数法的计算量很小，并且具有非常高的计算精度，因此特别适合应用于疲劳裂纹扩展分析。本文基于一维和二维裂纹权函数的统一形式，给出了利用权函数法进行应力强度因子计算时的数值积分方法，并基于半无限宽平板中角裂纹的应力强度因子的计算对数值积分方法进行验证，证实了该方法的有效性。

参考文献

[1] H. F. Bueckner, A novel principle for the computation of stressintensity factors. Z. angew. Math. Mech.50(9), 125-146 (1970).

[2] G. shen, G. Glinka, Universal features of weight functions forcracks in mode I. Engng Fracture Mech.40, 1135-1146 (1991).

[3] A. A. Moftakhar, G. Glinka. Calculation of stress intensityfactors by efficient integration of weight functions. Engng Fracture Mech. 43, 749-756 (1992).

[4] M. P. Stallybrass, A crack perpendicular to an elastic halfplane. Int. J. Engng Sci. 8, 351-362 (1970).

来源：FEM and FEA