CN104850683A - 基于弱形式求积元法计算材料裂纹尖端应力场系数的方法 - Google Patents

Abstract

本发明公开了一种基于弱形式求积元法计算材料裂纹尖端应力场系数的方法，它包括以下步骤：步骤1、根据分区广义变分原理，将含裂纹区域划分为势能区、余能区及其边界，建立分区广义变分方程；步骤2、建立势能区势能、余能区余能和势能区与余能区边界上的混合功的表达式，利用弱形式求积元法对这些表达式中的数值积分和微分进行离散近似；步骤3、运用变分驻值条件，得到含应力场系数的代数方程组；通过求解代数方程组即可直接得到应力场系数。本发明具有的优点是：计算推导过程更加直接、简明；对裂纹尖端的奇异性不需要进行特殊处理；通过增加单元内部节点数量提高计算精度，从而有效减少单元划分数量。

Description

基于弱形式求积元法计算材料裂纹尖端应力场系数的方法

技术领域

本发明属于含缺陷材料分析的技术领域，具体涉及一种材料裂纹尖端应力场系数的计算方法。

背景技术

材料裂纹尖端应力场系数是描述裂纹尖端应力强度及其分布情况的重要参数。上世纪50年代，M.L.Williams运用一系列特征展开式来表示裂纹尖端应力场，该展开式各项的系数即为裂纹尖端应力场系数。这些裂纹尖端应力场系数可以通过引入荷载和边界条件来求解得到。

在裂纹尖端应力场系数方法提出之初，由于只有一些几何形状和荷载条件、边界条件相对比较简单的问题可以得到解析解，因此其应用比较有限。后来随着计算机技术的快速发展，新的数值计算方法不断涌现，可以处理情况比较复杂的问题，因此逐渐成为求解裂纹尖端应力场系数的主要方法，并受到越来越多的关注和研究，其中应用最为广泛的当属有限元法(Finite Element Method)，因为它具有稳定高效、流程标准、便于程序化实现等优点。Karihaloo和Xiao通过运用有限元法，在裂纹尖端建立一种杂交裂纹单元(Hybrid Crack Element)，计算了多种边界条件下的应力场系数，计算结果具有较好的精度。Su和Feng、Su和Fok同样利用有限元法，在裂纹尖端建立分形有限单元(FractalFinite Element)，计算了多种情况下的应力场系数，也得到了比较好的计算精度和效率。Ayatollahi和Nejati则开发了一种超确定性有限元法(Finite Element Over-DeterministicMethod)。以上几种方法都基于有限元法，虽然能够得到比较准确的结果，但是由于都需要在裂纹尖端建立反映裂尖奇异性的特殊单元，因此造成理论推导比较繁琐、运用不够直接等缺点。近年来，新兴的无网格法(Meshless Method)得到快速发展，He等运用基于无网格伽辽金法(Element-Free Galerkin)的比例边界方法(Scaled Boundary Method)来计算应力场系数。但是，无网格法由于需要确定一些未知参数，如插值域的大小、背景积分域的大小等，因此具有计算量较大、效率不高等缺点。

发明内容

针对现有技术中存在的技术问题，本发明所要解决的技术问题就是提供一种基于弱形式求积元法(Weak Form Quadrature Element Method)计算材料裂纹尖端应力场系数的方法，它能准确、简单、直接地得到裂纹尖端应力场系数，从而为材料的工程应用提供参考。

本发明所要解决的技术问题是通过这样的技术方案实现的，它包括以下步骤：

步骤1、根据分区广义变分原理，将含裂纹区域划分为势能区、余能区及其边界，建立分区广义变分方程；

步骤2、建立势能区势能、余能区余能和势能区与余能区边界上的混合功的表达式，利用弱形式求积元法对这些表达式中的数值积分和微分进行离散近似；

步骤3、运用变分驻值条件，得到含应力场系数的代数方程组；通过求解代数方程组即可直接得到应力场系数。

与现有的技术相比，由于本发明将弱形式求积元法与分区广义变分原理相结合，使本发明具有以下显著的优点：

1、计算推导过程更加直接、简明。

2、对裂纹尖端的奇异性不需要进行特殊处理。

3、通过增加单元内部节点数量提高计算精度，从而有效减少单元划分数量。

附图说明

本发明的附图说明如下：

图1为含裂纹区域划分示意图；

图2为含中心穿透裂纹矩形钢板示意图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明：

本发明包括以下步骤：

步骤1、根据分区广义变分原理，将含裂纹区域划分为势能区、余能区及其边界，建立分区广义变分方程

如图1所示含裂纹区域，可将其划分为两个子区域，分别为势能区和余能区。余能区是以裂纹尖端为圆心，半径为r的圆形区域，其自变量为应力；势能区是余下的外围区域，其自变量为位移。

根据分区广义变分原理，参见文献Long,Y.Q.,Sub-region generalized principles in elasticity.Shanghai Journal of Mechanics,1981.22:p.1-9.(龙驭球.弹性力学中的分区广义变分原理.上海力学,1981.22:p.1-9.)，该系统的能量泛函∏可表示为：

∏＝∏_P-∏_C+∏_PC     (1)

式(1)中，∏_P为势能区势能，∏_C为余能区余能，∏_PC为两个子区域边界上余能区应力在势能区位移上所做的混合功。

步骤2、建立势能区势能、余能区余能和势能区与余能区边界上的混合功的表达式，利用弱形式求积元法对这些表达式中的数值积分和微分进行离散近似；

步骤1)、势能区势能∏_P的推导

图1中所示的势能区可根据问题需要将其划分为若干个求积单元，势能区的势能为所有求积单元的势能之和，即：

${\mathrm{\Π}}_P=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Π}}_p^s---\left(2\right)$

式(2)中，n为求积单元个数，为单个求积单元的势能。

每一个求积单元的势能又可以表示为应变能U^s与外力势能之和，即：

${\mathrm{\Π}}_P^s=U^s+U_Q^s---\left(3\right)$

利用弱形式求积元法，参见文献Zhong,H.Z.and Yu,T.,A weak form quadrature elementmethod for plane elasticity problems.Applied Mathematical Modelling,2009.33(10):p.3801-3814.(钟宏志,喻畑.平面弹性问题的弱形式求积元法.应用数学建模,2009.33(10):p.3801-3814.)，势能可转换成如下形式：

${\mathrm{\Π}}_P=\frac{1}{2}{d}^T\mathrm{Kd}-d^{T}Q---\left(4\right)$

式(4)中，K为整体刚度矩阵、Q为整体荷载向量、d为整体位移向量，d^T为d的转置。

步骤2)、余能区余能∏_C的推导

余能区余能∏_C可以表示为：

${\mathrm{\Π}}_C=\frac{1}{2}{\∫\∫}_{{\mathrm{\Ω}}_C}{\mathrm{\σ}}^T{D}^{-1}\mathrm{\σrdrd\θ}---\left(5\right)$

式(5)中，Ω_C为余能区积分区域，σ为应力向量，σ^T为σ的转置，D^-1为材料弹性模量逆矩阵，r和θ分别为以裂纹尖端为原点的极坐标系下的极半径和极角。σ应力向量可以用William特征展开式表示为：

σ＝F·R·A     (6)

式(6)中，F为关于θ的特征函数矩阵，R为关于r的极半径函数矩阵，A为裂纹尖端应力场系数向量。

将式(6)代入式(5)可得到：

${\mathrm{\Π}}_C=\frac{1}{2}{\∫\∫}_{{\mathrm{\Ω}}_C}{A}^T{R}^T{F}^T{D}^{-1}\mathrm{FRArdrd\θ}---\left(7\right)$

式(7)中，A^T为A的转置，R^T为R的转置，F^T为F的转置，其余符号意义与式(5)、式(6)相同。

将式(7)中应力场系数向量A、A^T提取到积分号外，即可得到：

${\mathrm{\Π}}_C=\frac{1}{2}{A}^T\mathrm{MA}---\left(8\right)$

式(8)中，M为与裂纹尖端应力场系数向量对应的积分矩阵，其表达式为：

$M={\∫\∫}_{{\mathrm{\Ω}}_C}{R}^T{F}^T{D}^{-1}\mathrm{FRrdrd\θ}---\left(9\right)$

步骤3)、势能区与余能区边界上的混合功∏_PC的推导

势能区与余能区边界上的混合功是由余能区在边界上的应力在相应势能区的边界位移上所做的功：

${\mathrm{\Π}}_{\mathrm{PC}}={\∫}_{{\mathrm{\Γ}}_{\mathrm{PC}}}{T}^T\mathrm{Urd\θ}---\left(10\right)$

式(10)中，∏_PC为势能区与余能区积分边界，U为势能区的边界位移向量，T为余能区的边界应力向量，T^T为T的转置。余能区的边界应力向量T可表示为：

$T=\stackrel{\‾}{\mathrm{\σ}}\mathrm{\Θ}---\left(11\right)$

式(11)中，Θ为关于极角θ的函数向量，特征展开式表示的应力矩阵。应力矩阵的各个分量也可以分解为：

${\stackrel{\‾}{\mathrm{\σ}}}_{i,j}=F_{i,j}\·R\·A,i,j=\mathrm{1,2}---\left(12\right)$

式(12)中，F_ij为关于θ的特征函数矩阵，i,j分别表示矩阵分量的行、列坐标，R和A的意义与式(6)相同。

将式(12)代入式(11)可得到：

T＝HA    (13)

式(13)中A为裂纹尖端应力场系数向量，H为与裂纹尖端应力场系数向量对应的系数矩阵。

将式(13)代入式(10)可得到：

${\mathrm{\Π}}_{\mathrm{PC}}={\∫}_{{\mathrm{\Γ}}_{\mathrm{PC}}}{A}^T{H}^T\mathrm{Urd\θ}---\left(14\right)$

式(14)中，A^T为A的转置，H^T为H的转置，其余符号意义与式(10)相同。

将式(14)中势能区的边界位移向量U转换为整体位移向量d即可得到：

∏_PC＝A^TWd       (15)

式(15)中，A^T为A的转置，d为整体位移向量，W为与裂纹尖端应力场系数向量和整体位移向量对应的矩阵，可表示为：

$W={\∫}_{{\mathrm{\Γ}}_{\mathrm{PC}}}{H}^T\mathrm{rd\θ}---\left(16\right)$

步骤3、运用变分驻值条件，得到含应力场系数的代数方程组；通过求解代数方程组即可直接得到应力场系数

由分区广义变分原理驻值条件有：

δ∏＝δ(∏_P-∏_C+∏_PC)＝0          (17)

将式(4)，式(8)，式(15)代入式(16)，得到：

$\mathrm{\δ}(\frac{1}{2}{d}^T\mathrm{Kd}-d^{T}Q-\frac{1}{2}{A}^T\mathrm{MA}+A^T\mathrm{Wd})=0---\left(18\right)$

式(17)中，δ()表示对括号内的函数进行一阶变分，因此：

δd^T(Kd-Q+W^TA)+δA^T(-MA+Wd)＝0        (19)

由于δd^T，δA^T是任意的，因此可得到含裂纹尖端应力场系数的代数方程组为：

$\left[\begin{array}{cc}-M& W\\ W^T& K\end{array}\right]\left\{\begin{array}{c}A\\ d\end{array}\right\}=\left\{\begin{array}{c}0\\ Q\end{array}\right\}---\left(20\right)$

式(19)中，W^T为W的转置，其余符号意义和前面相同。

求解以上方程组，即可直接得到材料裂纹尖端的应力场系数。

将裂纹尖端应力场系数代入Williams特征展开式中，即可得到裂纹尖端的应力场分布，从而预测含缺陷材料的剩余强度。

实施例：

飞机、高铁已经成为人们远途出行的主要交通工具，它的安全性能决定了人们的生命财产安全。然而，用于建造飞机、高铁的金属材料，在制造、加工以及使用过程中，不可避免的会出现许多微小裂纹，这些微小裂纹降低了材料的工作性能，因此，研究含缺陷材料的力学性能，具有重要的现实意义。

如图2所示的一块宽度(2w)和高度(2h)均为20cm，厚度为1cm的矩形薄钢板，其中心含有长度(2a)为0.4cm的穿透裂纹(近似于无限大板含有限长度裂纹)，裂纹倾斜角度为β。已知，钢板的弹性模量E＝206GPa，泊松比μ＝0.25。钢板一端固定，另一端受到均匀拉应力σ，其大小为10KPa。本实施例，在线弹性条件下，计算裂纹尖端应力场系数。

步骤1、在裂纹两侧尖端分别划分半径为r(r＝0.5a)的圆形区域作为余能区，余下区域根据需要划分为20个求积单元作为势能区。

步骤2、根据本发明，利用计算机语言FORTRAN编写程序计算能量泛函中的M,W(W^T),K和Q。

步骤3、对能量泛函变分后得到含应力场系数的代数方程组；求解代数方程组，直接得到裂纹尖端应力场系数A。

表1中，列出了本发明所获得的结果、精确解和常用算法的结果。

表1钢板裂纹尖端应力场系数

表1中

精确解是指无限大板的精确解。

[1]一种基于有限元法的超确定性方法。参见Ayatollahi,M.R.and Nejati,M.,Anover-deterministic method for calculation of coefficients of crack tip asymptotic field fromfinite element analysis.Fatigue&Fracture of Engineering Materials&Structures,2010.34(3):p.159-176.(Ayatollahi,M.R.与Nejati,M.，一种基于有限单元法求解裂纹尖端渐近场系数的超确定性方法，工程材料及结构的疲劳与断裂,2010.34(3):p.159-176)。

[2]一种基于有限元法的超级奇异单元法。参见Tsang,D.K.L.and Oyadiji,S.O.,Super singular element method for two-dimensional crack analysis.Proceedingsof the Royal Society A-Mathematical Physical and Engineering Sciences,2008.464(2098):p.2629-2648.(Tsang,D.K.L.与Oyadiji,S.O.二维裂纹分析的超奇异元法.英国皇家学会数学、物理与工程科学学报)。

[3]一种基于有限元法的杂交裂纹单元法。参见Xiao,Q.Z.and Karihaloo,B.L.,Implementation of hybrid crack element on a general finite element mesh and in combinationwith XFEM.Computer Methods in Applied Mechanics and Engineering,2007.196(13-16):p.1864-1873.(Xiao,Q.Z.与Karihaloo,B.L.，一般有限元网格结合XFEM的杂交裂纹单元的实现.应用力学和工程技术中的计算机方法，2007.196(13-16):p.1864-1873.)。

对比β＝0°、15°、30°三种不同情况下的计算结果，可以看出：本发明所得到的裂纹尖端应力场系数更趋近于精确解，而且，本发明在同一种单元划分系统下，计算三种情况下的裂纹尖端应力场系数，具有非常好的便捷性和稳定性。

Claims (5)

1.基于弱形式求积元法计算材料裂纹尖端应力场系数的方法，其特征是：包括以下步骤：

步骤1、根据分区广义变分原理，将含裂纹区域划分为势能区、余能区及其边界，建立分区广义变分方程；

步骤2、建立势能区势能、余能区余能和势能区与余能区边界上的混合功的表达式，利用弱形式求积元法对这些表达式中的数值积分和微分进行离散近似；

步骤3、运用变分驻值条件，得到含应力场系数的代数方程组；通过求解代数方程组即可直接得到应力场系数。

2.根据权利要求1所述的基于弱形式求积元法计算材料裂纹尖端应力场系数的方法，其特征是，在步骤2中，势能区的势能Π_p为：

${\mathrm{\Π}}_p=\frac{1}{2}{d}^T\mathrm{Kd}-d^{T}Q$

式中，K为整体刚度矩阵，Q为整体荷载向量，d为整体位移向量，d^T为d的转置。

3.根据权利要求2所述的基于弱形式求积元法计算材料裂纹尖端应力场系数的方法，其特征是，在步骤2中，余能区的余能Π_C为：

${\mathrm{\Π}}_C=\frac{1}{2}{A}^T\mathrm{MA}$

式中，A为裂纹尖端应力场系数向量，A^T为A的转置，M为与裂纹尖端应力场系数向量对应的积分矩阵。

4.根据权利要求3所述的基于弱形式求积元法计算材料裂纹尖端应力场系数的方法，其特征是，在步骤2中，势能区与余能区边界上的混合功Π_PC为：

Π_PC＝A^TWd

式中，W为与裂纹尖端应力场数向量和整体位移向量对应的矩阵。

5.根据权利要求4所述的基于弱形式求积元法计算材料裂纹尖端应力场系数的方法，其特征是，在步骤3中，含应力场系数的代数方程组为：

$\left[\begin{array}{cc}-M& W\\ W^T& K\end{array}\right]\left\{\begin{array}{c}A\\ d\end{array}\right\}=\left\{\begin{array}{c}0\\ Q\end{array}\right\}.$