CN110866350B - 一种裂纹尖端单元部分增强的扩展有限元方法 - Google Patents

Abstract

本发明涉及断裂力学领域，尤其涉及一种裂纹尖端单元部分增强的扩展有限元方法，S1，建立裂纹分析模型；S2，获得裂纹尖端单元部分增强的扩展有限元方法的矩阵方程；S3，施加边界条件和载荷，求解系统方程，获得整个模型的位移场；S4，在获得位移场的基础上，采用交叉相积分的方法，获得裂纹尖端处的应力强度因子。该发明的优点在于：裂纹尖端单元采用部分增强的位移场，对于裂纹撕裂部分采用增强位移场，模拟裂纹面的位移不连续性；对于未撕裂部分，采用有限元连续位移场，可以实现裂纹尖端单元位移场特征的模拟。

Description

一种裂纹尖端单元部分增强的扩展有限元方法

技术领域

本发明涉及断裂力学领域，尤其涉及一种裂纹尖端单元部分增强的扩展有限元方法。

背景技术

材料和结构破坏会引起事故、甚至造成巨大生命财产损失，材料和结构的破坏很多是由裂纹的存在导致，为了减少由于裂纹带来的材料或结构破坏，断裂力学应运而生。数值方法已经成为了断裂力学分析中不可或缺的工具，常用的断裂力学数值方法有奇异单元法、光滑有限元法、扩展有限元等方法。

在这些方法中，奇异单元法和光滑有限元通过重复节点的方法来模拟裂纹面位移不连续性。网格划分时，裂纹位于单元边界，且裂纹扩展时，网格需要重新划分，因此网格划分工作极为复杂和繁琐。在扩展有限元中，采用跳跃函数和尖端位移场函数来表征裂纹面的位移不连续和裂纹尖端位移奇异性。由于裂纹尖端位移场函数的使用会带来奇异项积分问题，因此不少学者在扩展有限元中不采用裂纹尖端位移场函数增强，而只采用跳跃函数(以下简称“XFEM-H”)。如果只采用跳跃函数，裂纹尖端需要位于单元边界上，这一限制条件给网格划分造成了极大的不便，特别是对于裂纹扩展时，裂纹形貌发生改变，要保证裂纹尖端始终位于单元边界，十分困难。

发明内容

针对采用跳跃函数增强的扩展有限元(以下简称“XFEM-H”)，在模拟裂纹时存在的不足，本发明提供了一种裂纹尖端单元部分增强的扩展有限元方法(以下简称“XFEM-P”)。在该方法中，裂纹尖端不需要位于单元边界，这给网格划分工作带来了极大的便利，同时该方法相较于XFEM-H可以提供更准确的计算结果和更高的计算效率。

为实现上述目的，本发明采用以下技术方案：

一种裂纹尖端单元部分增强的扩展有限元方法，裂纹尖端单元采用部分增强的位移场，该方法包括如下步骤：

S1，建立裂纹分析模型；

S2，获得裂纹尖端单元部分增强的扩展有限元方法的矩阵方程；

S3，施加边界条件和载荷，求解系统方程，获得整个模型的位移场；

S4，在获得位移场的基础上，采用交叉相积分的方法，获得裂纹尖端处的应力强度因子。

本发明的优点在于：

(1)本发明中，裂纹尖端单元采用部分增强的位移场，对于裂纹撕裂部分采用增强位移场，模拟裂纹面的位移不连续性；对于未撕裂部分，采用有限元连续位移场，可以实现裂纹尖端单元位移场特征的模拟。

(2)本发明中，裂纹尖端可以位于单元的任意处，消除了裂纹尖端必须位于单元边界上的限制，解决了跳跃函数增强扩展有限元中网格划分中的常见问题。

(3)本发明提出的断裂力学计算方法，比传统的XFEM-H方法提供更准确的计算结果。

附图说明

图1为本发明的流程示意图；

图2为本发明模型中的单元分类和跳跃函数增强节点示意图；

图3为裂纹面的外法线示意图。

图4为开裂单元和裂纹尖端单元示意图；

图5-a、图5-b、图5-c为裂纹尖端单元的三种不同区域划分方法；

图6为采用交叉相积分时的相交单元的示意图；

图7-a为含裂纹平板拉伸试样示意图；

图7-b为XFEM-P网格示意图

图7-c为XFEM-H网格示意图；

图8为本发明所提供方法的计算结果与XFEM-H和解析解的对比图；

具体实施方式

本实施例采用含裂纹的平板拉伸模型，如图7-a所示，该模型上端面施加拉伸载荷σ＝1Mpa，底边为脚支，左底脚为固支。高度H＝2mm，裂纹长度a＝0.3mm，板宽b＝1mm，材料属性：杨氏模量E＝1×10³Mpa，泊松比υ＝0.3。假设平板为平面应变状态。

如图1所示，对该实施例采用尖端单元部分增强的扩展有限元方法进行计算，主要步骤如下：

S1，建立裂纹分析模型；

具体地步骤如下：

S11，对几何区域进行网格划分，分成若干个单元；本方案中分成两种不同的网格划分方案，如图7-b所示，图节点布置为22×42，为XFEM-P所用网格图，图7-c节点布置为21×42为XFEM-H所用网格。

S12，在网格中确定裂纹面和裂纹尖端的位置；如图7-b和图7-c所示，在网格中分别确定裂纹面和裂纹尖端的位置。

S13，根据裂纹面和裂纹尖端所在位置，将网格中的单元分为裂纹贯穿整个单元的开裂单元、裂纹尖端所在区域的裂纹尖端单元以及无裂纹的普通单元；如图2所示。

S14，其中开裂单元的所有节点为跳跃函数增强节点，其它节点均为普通节点；跳跃函数增强节点对应的跳跃函数定义如下：

式中x是样点，x^*是裂纹面上离x最近的点，n是裂纹面的外法线方向，如图3所示。

S2，获得裂纹尖端单元部分增强的扩展有限元方法的矩阵方程；

具体步骤如下：

S21、分别获取开裂单元、裂纹尖端单元、普通单元的刚度矩阵K_i，i为单元编号：

获取开裂单元的的刚度矩阵K_i的具体步骤如下：

SA211、裂纹将开裂单元分成第一区域和第二区域，第一区域和第二区域分别为图4中的3'-3-5-4'区域和3'-4'-6-4区域，对应开裂单元的位移场由下式给出：

式中

为有限元位移场部分，N为普通节点集合，/>

为跳跃函数增强部分，N^c为增强节点集合，u(x)为位移场，u_m为节点位移，a_n为跳跃函数增强自由度，x_n是节点坐标，m和n为节点编号，N_m(x)和N_n(x)分别为m节点和n节点的节点形函数；

SA212、通过位移场获得相应的应变向量，如下式所示，

ε＝Βu

式中ε为应变向量，u为单元位移集合向量，由节点位移u_m和增强节点位移a_n构成，应变-位移矩阵Β如下式所示

当节点为跳跃函数增强节点时，B_l＝N_l(x)(H(x)-H(x_l))，x_l为节点坐标，当节点为普通节点时B_l＝N_l(x)，l为矩阵组装指数，N_node为节点数目，包括普通节点和跳跃函数增强节点；

SA213、获得开裂单元的刚度矩阵，如下式所示，

式中，N_dom为i单元包含的区域数，N_gau为各区域里面的高斯积分点数目，W_t为t积分点对应的权重，r和t为求和指数，Ω_i为单元积分区域；D为材料的刚度矩阵，对于平面应变问题，D可表示为：

式中E为杨氏模量，υ为泊松比；

获取裂纹尖端单元的的刚度矩阵的具体步骤如下：

SB211、在裂纹尖端处作裂纹的垂直线，分别与裂纹尖端单元边界相交于点1'和2'。

SB212、根据裂纹面和裂纹尖端在裂纹单元的相对位置，且根据裂纹面和裂纹尖端在裂纹尖端单元的位置，将对应的裂纹尖端单元分成撕裂区和完整区；分为以下三种情况：

1)如图5-a所示，裂纹尖端垂直线1'-2'将裂纹尖端单元分为三角形区域1'-2'-4和五边形区域1-2-3-2'-1'，其中三角形区域1'-2'-4可由裂纹继续划分为三角形2'-3'-Pt和四边形区域1-Pt-3'-4区域，由于这两个区域被裂纹所撕裂，所以称为裂纹撕裂区。为了方便积分计算，将五边形区域1-2-3-2'-1'分为三角形1'-2'-3和四边形1-2-3-1'。

2)如图5-b所示，裂纹尖端垂直线1'-2'将裂纹尖端单元分为四角形区域1'-2'-2-1和1'-2'-3-4，其中1'-2'-3-4可由裂纹继续划分为四边形区域1-Pt-3'-4和2'-3-3'-Pt区域，为裂纹撕裂区。为了方便积分计算，将五边形区域1-2-3-2'-1'分为三角形1'-2'-3和四边形1-2-3-1'。

3)如图5-c所示，裂纹尖端垂直线1'-2'将裂纹尖端单元分为三角形1'-2'-2和五边形1'-2'-3-4-1，其中五边形区域1'-2'-3-4-1可由裂纹继续划分为四边形2'-3-3'-Pt和五边形区域1'-Pt-3'-4-1区域，为了方便积分计算，将五边形区域1'-Pt-3'-4-1分为三角形Pt-3'-4和四边形1-1'-Pt-4，其中Pt-3'-4，1-1'-Pt-4和2'-3-3'-Pt为裂纹撕裂区。

以上三种情况，对于撕裂区，其位移场由下式给出

对于完整区，其位移场由下式给出

SB213、根据其位移场，产生应变-位移矩阵Β，对于撕裂区的应变-位移矩阵Β包含普通节点项(N_l(x))和增强项(N_l(x)(H(x)-H(x_l)))，对于完整区，只包含普通节点项；

SB214、撕裂区和完整区叠加最终形成裂纹尖端单元的刚度矩阵；

获取普通单元的步骤为：

对于模型中，与裂纹没有任何交叉的普通单元，其位移场与以上未撕裂区域一致。根据其位移场，产生普通单元的应变-位移矩阵Β，由于该类型单元与裂纹不交叉，所以无需划分区域，即积分区域数为1。

S22、获得以上三种单元刚度矩阵的基础上，组装单元矩阵，形成系统总体刚度矩阵，如式所示：

式中K为系统刚度矩阵，N_elem为模型中的单元数。

S3，施加边界条件和载荷，求解系统方程，获得整个模型的位移场；具体步骤为：

S31、根据下式，获得节点力向量f；

式中f为节点力向量，b为体积力向量，Γ_t为受分布载荷边界，t_Γ为分布载荷向量,Ω为模型区域,N^T(x)为形函数N(x)的转置矩阵；

S32、获得系统方程，并求解获得整个模型的位移场，如下式所示

Ku＝f。

S4，在获得位移场的基础上，采用交叉相积分的方法，获得裂纹尖端处的应力强度因子。

步骤S4具体如下所述：

S41、以裂纹尖端为圆心，取单元平均尺寸的3倍为半径画圆，如图6所示，与圆相交的单元，用N^eff表示；

S42、求取所求问题和辅助场的交叉相积分,在求解过程中，需要借用辅助场；公式如下：

式中，变量的右上标括弧里面的数字1和2，分别表示所求问题和辅助场；I^(1,2)为所求问题与辅助场交叉相积分，n为求和指数，A_eff,s为有效单元的积分区域，下标i,j,k表示三维坐标分量，

为所求问题位移i方向分量对x导数，x_j为j方向坐标分量，A为积分面积，q在圆圈内的节点为1，圆圈外的节点为0；当j＝x时，δ_xj＝1；其他情况下δ_xj＝0；/>

分别为所求问题应力在i平面上沿着j和k方向上的分量，裂纹尖端应力σ由下式给出：

σ＝Dε

式中ε为应变向量，由ε＝Bu给出；

为辅助位移i方向分量对x导数，辅助位移由下式给出，

式中

和/>

分别为辅助位移在1,2方向上的分量，(r,θ)为裂纹尖端局部极坐标，μ＝E/2(1+ν)为剪切模量,平面应力状态下κ＝(3-υ)/(1+υ)，平面应变状态下，κ＝3-4υ，

为辅助场的I型和II型应力强度因子，/>

根据链式法则求出，如下式所示；

式中

为辅助应变在i平面上沿着k方向上的分量，可由下式给出

为辅助应力在i平面上沿着j方向上的分量，可由下式给出

S43,由下式给出所求问题的I型和II型应力强度因子；

式中

为所求问题的I型和II型应力强度因子，平面应力状态下E^*＝E，平面应变状态下E^*＝E/(1-υ)，I^{(1,Mode I)}为所求问题与辅助场/>

时交叉相积分结果，I^{(1,Mode II)}为所求问题与辅助场/>

时的积分结果。

计算结果如图8所示，图中横坐标为模型中的节点数，纵坐标正则化的I型应力强度因子，图中XFEM-H表示由XFEM-H提供的应力强度因子，XFEM-P表示有XFEM-P提供的应力强度因子，解析解为该问题的精确解。从图中可以看出，在相同的网格密度下，XFEM-P给出的计算结果比XFEM-H更接近精确解，即本发明提供的方法精度更高。同时从网格可以看出，在XFEM-H中，裂纹尖端需要位于单元边界上，而XFEM-P中，裂纹尖端可以位于单元任意地方。

以上仅为本发明创造的较佳实施例而已，并不用以限制本发明创造，凡在本发明创造的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明创造的保护范围之内。

Claims (4)

1.一种裂纹尖端单元部分增强的扩展有限元方法，其特征在于，裂纹尖端单元采用部分增强的位移场，该方法包括如下步骤：

S1，建立裂纹分析模型；

S2，获得裂纹尖端单元部分增强的扩展有限元方法的矩阵方程；

S3，施加边界条件和载荷，求解系统方程，获得整个模型的位移场；

S4，在获得位移场的基础上，采用交叉相积分的方法，获得裂纹尖端处的应力强度因子；

步骤S3包括以下步骤：

S31、根据下式，获得节点力向量f；

式中f为节点力向量，b为体积力向量，Γ_t为受分布载荷边界，t_Γ为分布载荷向量,Ω为模型区域,N^T(x)为形函数N(x)的转置矩阵；

S32、获得系统方程，并求解获得整个模型的位移场，如下式所示

Ku＝f；

步骤S4具体如下所述：

S41、以裂纹尖端为圆心，取单元平均尺寸的若干倍为半径画圆，与圆相交的单元，用N^eff表示；

S42、求取所求问题和辅助场的交叉相积分,在求解过程中，需要借用辅助场；公式如下：

式中，变量的右上标括弧里面的数字1和2，分别表示所求问题和辅助场；I^(1,2)为所求问题与辅助场交叉相积分，s为求和指数，A_eff,s为有效单元的积分区域，下标i,j,k表示三维坐标分量，

为所求问题位移i方向分量对x导数，x_j为j方向坐标分量，A为积分面积，q在圆圈内的节点为1，圆圈外的节点为0；当j＝x时，δ_xj＝1；其他情况下δ_xj＝0；/>

分别为所求问题应力在i平面上沿着j和k方向上的分量，裂纹尖端应力σ由下式给出：

σ＝Dε

式中ε为应变向量，由ε＝Bu给出；

为辅助位移i方向分量对x导数，辅助位移由下式给出，

式中

和/>

分别为辅助位移在1,2方向上的分量，(r,θ)为裂纹尖端局部极坐标，μ＝E/2(1+ν)为剪切模量,平面应力状态下κ＝(3-υ)/(1+υ)，平面应变状态下，κ＝3-4υ，

为辅助场的I型和II型应力强度因子，/>

根据链式法则求出，如下式所示；

式中

为辅助应变在i平面上沿着k方向上的分量，可由下式给出

为辅助应力在i平面上沿着j方向上的分量，可由下式给出

S43,由下式给出所求问题的I型和II型应力强度因子；

式中

为所求问题的I型和II型应力强度因子，平面应力状态下E^*＝E，平面应变状态下E^*＝E/(1-υ)，I^{(1,Mode I)}为所求问题与辅助场/>

时交叉相积分结果，I^{(1,Mode II)}为所求问题与辅助场/>

时的积分结果。

2.根据权利要求1所述的一种裂纹尖端单元部分增强的扩展有限元方法，其特征在于，步骤S1包括以下步骤：

S11，对几何区域进行网格划分，分成若干个单元；

S12，在网格中确定裂纹面和裂纹尖端的位置；

S13，根据裂纹面和裂纹尖端所在位置，将网格中的单元分为裂纹贯穿整个单元的开裂单元、裂纹尖端所在区域的裂纹尖端单元以及无裂纹的普通单元；

S14，其中开裂单元的所有节点为跳跃函数增强节点，其它节点均为普通节点；跳跃函数增强节点对应的跳跃函数定义如下：

式中x是样点，x^*是裂纹面上离x最近的点，n是裂纹面的外法线方向。

3.根据权利要求2所述的一种裂纹尖端单元部分增强的扩展有限元方法，其特征在于，步骤S2包括以下步骤：

S21、分别获取开裂单元、裂纹尖端单元、普通单元的刚度矩阵K_i，i为单元编号：

S22、获得以上三种单元刚度矩阵的基础上，组装单元矩阵，形成系统总体刚度矩阵，如式所示：

式中K为系统刚度矩阵，N_elem为模型中的单元数。

4.根据权利要求3所述的一种裂纹尖端单元部分增强的扩展有限元方法，其特征在于，获取开裂单元的刚度矩阵K_i的具体步骤如下：

SA211、裂纹将开裂单元分成第一区域和第二区域，对应开裂单元的位移场由下式给出：

式中

为有限元位移场部分，N为普通节点集合，/>

为跳跃函数增强部分，N^c为增强节点集合，u(x)为位移场，u_m为节点位移，a_n为跳跃函数增强自由度，x_n是节点坐标，m和n为节点编号，N_m(x)和N_n(x)分别为m节点和n节点的节点形函数；

SA212、通过位移场获得相应的应变向量，如下式所示，

ε＝Βu

式中ε为应变向量，u为单元位移集合向量，由节点位移u_m和增强节点位移a_n构成，应变-位移矩阵Β如下式所示

当节点为跳跃函数增强节点时，B_l＝N_l(x)(H(x)-H(x_l))，x_l为节点坐标，当节点为普通节点时B_l＝N_l(x)，l为矩阵组装指数，N_node为节点数目，包括普通节点和跳跃函数增强节点；

SA213、获得开裂单元的刚度矩阵，如下式所示，

式中，N_dom为i单元包含的区域数，N_gau为各区域里面的高斯积分点数目，W_t为t积分点对应的权重，r和t为求和指数，Ω_i为单元积分区域；D为材料的刚度矩阵，对于平面应变问题，D可表示为：

式中E为杨氏模量，υ为泊松比；

获取裂纹尖端单元的刚度矩阵的具体步骤如下：

SB211、在裂纹尖端处作裂纹的垂直线，分别与裂纹尖端单元边界相交于点1'和2'；

SB212、根据裂纹面和裂纹尖端在裂纹尖端单元的位置，将对应的裂纹尖端单元分成撕裂区和完整区；

对于撕裂区，其位移场由下式给出

对于完整区，其位移场由下式给出

SB213、根据其位移场，产生应变-位移矩阵Β，对于撕裂区的应变-位移矩阵Β包含普通节点项N_l(x)和增强项N_l(x)(H(x)-H(x_l))，对于完整区，只包含普通节点项；

SB214、撕裂区和完整区叠加最终形成裂纹尖端单元的刚度矩阵；

获取普通单元的步骤为：

SC21、位移场与上述完整区一致，根据其位移场，产生普通单元的应变-位移矩阵Β。

Images