非连续问题中单元分割的模板方法1)

王理想 文龙飞 肖桂仲,∗∗ 田荣,2)

∗(中物院高性能数值模拟软件中心,北京 100088)

†(北京应用物理与计算数学研究所,北京 100088)

∗∗(南京理工大学,南京 210094)

引言

非连续问题广泛存在于基础研究[1-4]和工程应用[5-8]中,如裂纹、孔洞、夹杂、多材料、流固耦合、相变、剪切带等.研究这类问题,具有重要的科学意义和广泛的应用价值.

传统有限元法模拟这类问题面临诸多困难,如网格划分需与非连续面对齐、非连续面更新(如裂纹扩展)时需网格重分重映、空间梯度变化较大处(如裂尖处)需进行网格加密等.

鉴于传统有限元法局限性,1990 年代开始,发展出一批基于单位分解(PU)[9]的数值方法,如单位分解有限元(PUFEM)[10]、广义有限元(GFEM)[11-15]、数值流形法(NMM)[16-20]、扩展有限元(XFEM)[21-25]等等.这类方法通过在基函数中引入先验的加强函数,可有效求解各类连续、非连续问题.

特别针对裂纹问题,1999 年美国西北大学Ted Belytschko 研究组提出XFEM[21-22],在学术界获得广泛关注和迅速发展.大型商业有限元软件纷纷添加XFEM 模块,标志该方法在工业界获得认可.

XFEM 通过引入C−1型局部加强函数,在函数空间逼近裂纹的非连续性.由此,裂纹几何独立于模拟网格,裂纹扩展时亦无需网格重分重映.通过引入解析性质的加强函数,其计算精度亦获得大幅提高.

XFEM 在模拟裂纹扩展方面取得巨大成功,但同时也面临两大困扰:总体刚度矩阵高度病态和动力学计算时额外自由度上能量无法正确传递[23-24].前者导致迭代法求解收敛缓慢甚至不收敛;后者导致动力学求解实施困难.有鉴于此,作者团队基于无额外自由度的裂尖插值格式,提出一种改进型XFEM(improved XFEM)[26-31],成功克服上述困难.

XFEM 类方法在求解非连续问题时,面临非连续单元积分问题——此类单元受到非连续面的分割,若使用常规高斯积分,精度将严重损失,故需进行特殊处理.目前有以下几类处理方法[23,32-33].

(1)高阶高斯积分法[34-35]:在单元内布置大量高斯积分点(如6×7[34],8×8[35]),判断积分点与裂纹面的位置关系,从而对积分点所在积分域进行积分.该方法实现简单,但精度较低[11].

(2)自适应积分法[11,36-37]:将初始非连续单元进行网格细分,在细分后的格子上自适应布置高斯积分点,然后进行非连续单元积分.该方法实现较简单,精度获得提高,但需要布置大量高斯积分点.

(3)子网格积分法[38-39]:通过子网格剖分,将非连续单元与裂纹进行几何分割,非连续单元积分在子网格上实现.由于该方法一般引入精确几何操作,精度较高,但实现复杂、鲁棒性差.

(4)矩量拟合(moment fitting)法[40-42]:通过一族函数(如幂函数)在非连续单元上的精确积分,构造一组非线性方程组,求解该方程组获得积分点位置和权重,再使用所求得的积分点和权重进行非连续单元积分.该方法不需要对单元进行分割,仅需引入少量积分点就可取得较高精度,但由于每个单元都需要求解非线性方程组,导致求解效率降低.

(5)特殊积分法:这类方法的思想是将原积分进行一定特殊处理后再进行积分,具有实现简便或精度高的优点,但通用性一般.这类方法包括:①等效法.例如:Ventura[43]将Heaviside 函数等效为二次多项式,Abedian 和Düster[44]将非连续函数等效为Legendre 多项式;由此把非连续积分转换为连续积分进行计算.②降维法.例如:Sudhakar 和Wall[45]通过高斯定理将体(面)积分转换为面(线)积分.③变换法.例如:Mousavi 和Sukumar[46]使用广义Duffy 变换,将在三角形(金字塔)上积分变换为在正方形(立方体)上积分.此外,还有其他特殊积分方法[47-48],不再一一赘述.

上述几类方法在精度、效率、实现或通用性方面存在各自的优点和缺点.本文从工程实用化角度出发,提出一种基于水平集的非连续单元模板分割(单元积分)方法.该方法采取模板查询的方式进行非连续单元子剖分,在三角形子网格上布置高斯积分点.相比于精确子剖分方法,该方法避免复杂几何操作,可提高计算效率和鲁棒性.将该方法与常规XFEM、改进型XFEM 进行结合,应用于孔洞、夹杂、裂纹等非连续问题分析中,从而为XFEM 类方法在实际工程应用中提供有效支撑.

1 控制方程

1.1 强形式

如图1 所示非连续静力问题,应满足平衡方程

式中,∇为梯度算子,σ为应力张量,b为体力.

图1 二维非连续静力问题描述Fig.1 Illustration of statics with discontinuities in 2D

位移边界条件和应力边界条件分别如下所示

式中,u为位移,为在位移边界Γu处给定的位移值,为在应力边界Γt处给定的应力值,Γ0为无应力边界,Γc为裂纹面,Γh为孔洞面,n为各个面上单位外法向量.

对于线弹性问题,其本构关系如下

式中,D为弹性张量,ε为应变张量.

小变形下,应变ε和位移u满足如下几何关系

1.2 弱形式

根据虚功原理,平衡方程(1)与应力边界条件(3)可转化为如下弱积分形式

式中,δu为虚位移,δε为虚应变.

将几何关系(5)代入平衡方程弱积分式(6),可得

2 非连续界面水平集描述

本文使用水平集法[49]隐式描述非连续界面,如孔洞面、夹杂面、裂纹面、多材料界面等.非连续界面可分为强、弱两种:裂纹面属于强非连续界面,夹杂面和多材料界面属于弱非连续界面.根据不同类型加强形式,孔洞面可以是强非连续界面,也可以是弱非连续界面.本文不对其进行展开讨论.

2.1 非连续面水平集函数φ(x)

非连续面水平集函数φ(x)定义为求解域Ω 内任意给定一点x∈Ω 到非连续面Γd的含符号距离,其数学表达式为

式中,||·||为欧式范数,n为非连续面单位外法向量,sign(·)为符号函数.上式中的非连续面Γd可以是裂纹面Γc、夹杂面Γi、孔洞面Γh或多材料界面Γm.

2.2 裂纹尖端水平集函数ψ(x)

裂纹尖端水平集函数ψ(x)定义为求解域Ω 内任意给定一点x∈Ω 到裂纹尖端xt的含符号距离,其数学表达式为

式中,t为沿裂纹扩展方向的单位向量.针对裂纹所构造的水平集函数,如图2 所示.

图2 裂纹水平集函数构造(根据文献[33]重画)Fig.2 Construction of level set functions for a crack(redrawn after Ref.[33])

2.3 裂纹尖端极坐标系(r,θ)

裂纹可通过水平集函数φ(x)和ψ(x)进行描述.除此之外,还需要定义裂纹尖端极坐标系(r,θ),其数学表达式可通过水平集函数给出

3 扩展有限元

3.1 常规扩展有限元:弱非连续问题

对于弱非连续问题,基于常规XFEM 的位移场逼近可表示为

3.2 常规扩展有限元:强非连续问题

对于裂纹问题,除了考虑非连续面(裂纹面)阶跃加强之外,还需考虑裂尖奇异加强.此时,基于常规XFEM 的位移场逼近为

式中,I 为全部节点集合,J 为裂纹面加强节点集合,K 为裂尖加强节点集合;i,j,k,α 为序数;为常规自由度,为裂纹面加强额外自由度,为裂尖加强额外自由度;N(x)为标准有限元形函数,Fjump(x)为阶跃加强函数,Ftip(x)为裂尖奇异加强函数.Fjump(x)常使用如下Heaviside 函数

对于各向同性线弹性材料,Ftip(x)可取为

3.3 改进型扩展有限元:裂纹问题

在使用改进型扩展有限元(IXFEM)求解时,仅考虑裂纹问题,其位移场逼近形式为[26-31]

式中,I 为所有节点集合,J 为裂纹面加强节点集合,K 为裂尖加强节点集合,Kk⊂K 为节点k处裂尖加强影响域内节点集合;i,j,k,m为序数;为常规自由度,为裂纹面加强额外自由度;N(x)为标准有限元形函数,Fjump(x)为阶跃加强函数,与式(15)定义相同,(x)为裂尖局部加强函数

其中,p(x)为基向量,pk为p(xk)的缩写;A为矩量矩阵,为A−1的第一列,为的第一个元素;δ为Kronecker 符号.p(x)和A的表达式如下所示

式中,c为常数,一般取为2 ∼3;h为网格尺寸;Ftip(x)为裂尖奇异加强函数,其定义与式(16)相同.

观察式(17)与式(14)可看出,改进型XFEM 与常规XFEM 的区别仅在于裂尖加强项.改进型XFEM裂尖加强项中不再包含额外自由度,可有效解决常规XFEM 线性相关问题.

3.4 混合单元修正

混合单元为同时包含有限元节点与加强节点的单元(如图3 所示).为保证XFEM 算法精度和收敛性,混合单元形函数需加以修正[51].

图3 混合单元处理:处理之前(左);处理之后(右)Fig.3 Blending element treatment:Before(left);after(right)

引入如下斜坡函数[52]

则在混合单元Ωbld上,位移逼近具有如下加权形式

4 控制方程离散

通过Galerkin 法,并应用位移逼近式(12)或式(14)或式(17),可将平衡方程弱积分式(7)离散为如下线性方程组形式

式中,K为刚度矩阵,U为自由度向量,F为载荷向量.对于不同位移逼近,K,U,F形式略有不同.

若采用式(12)进行离散,式(23)具有如下形式

若采用式(14)进行离散,式(23)具有如下形式

若采用式(17)进行离散,式(23)具有如下形式

式(24)∼式(26)中K和F可统一表示为

其中,α,β=u,a或u,b,c或u,b,表示所采用的位移逼近方式;B为应变矩阵,N为形函数矩阵.

5 非连续单元分割方法

本方法的核心思想是,根据水平集在单元内的近似属性,假定裂纹在单元内平直分布,由此根据水平集值特定形态对非连续单元进行分割.该方法由于采用模板查询的方式,相比于精确分割,可避免复杂几何操作,从而提高计算效率和鲁棒性.

所分割的非连续单元分为两类:一类是被非连续面完全贯穿的单元,称为“切割单元”;一类是被非连续面(即裂纹面)部分贯穿的单元,称为“裂尖单元”.下面分别给出这两类单元的分割子剖分.

5.1 切割单元

5.1.1 三角形单元

首先,遍历三角形切割单元水平集值所有形态并建立标准单元分割模板库,如表1 和图4 所示.

表1 三角形切割单元水平集值符号形态Table 1 LSV sign patterns of triangular cut elements

图4 三角形切割单元分割形态Fig.4 The partitioning patterns of triangular cut elements

表1 中,case 表示所有水平集值符号的可能情况,共33=27 种;status 表示单元水平集值符号个数所组成的三元组,每个三元组表示一种状态,共10种;type 表示切割单元类型,共3 种.每种类型可分为(a)、(b)两种亚型,由sign 进行标识.三角形切割单元的3 种切割类型及其亚型,如图4 所示.

然后,根据三角形切割单元水平集值,对非标准单元进行形态查询(查询type 和sign)和模板插值.对棱ij上交点p的位置进行线性插值

其中,λ 可由水平集值插值得到

交点p的坐标可由i和j的坐标和水平集值给出

最后,套用标准单元分割模板实现单元分割和子三角剖分,如算法1 所述.序号ijk的具体排列根据不同切割单元类型预先排定,可通过查表得到.

5.1.2 四边形单元

首先,遍历四边形切割单元水平集值所有形态并建立标准单元分割模板库.四边形切割单元中,共34=81 种情况、15 种状态.

在不考虑一些特殊情况(如裂纹在单元内发生较大转折、两条裂纹同时切割一个单元等)的假设下,四边形切割单元可分为5 种类型,每个类型分为(a)、(b)两种亚型,如图5 所示.

图5 四边形切割单元分割形态Fig.5 Partitioning patterns of quadrilateral cut elements

然后,根据四边形切割单元水平集值,对非标准单元进行形态查询(查询type 和sign)和模板插值,其中模板插值与三角形切割单元中的实现一致.

最后,套用标准单元分割模板实现单元分割和子三角剖分,如算法2 所述.与三角形单元类似,四边形单元中的序号ijkl具体排列也已根据不同切割单元类型预先排定,可通过查表得到.

5.2 裂尖单元

仅考虑裂尖在单元内的情况.裂尖单元分割形态如图6 所示,分别对其进行子剖分,详见算法3.

值得注意的是,本文方法在实现过程中:(1)预先将若干标准单元进行子三角化,形成标准单元分割模板库,并在计算机程序中进行存储;这一存储过程发生在程序编译期,不占用计算时间.(2)在处理任意单元分割时,仅需查询存储在计算机中的标准单元分割模板库进行匹配,就可对任意单元进行分割;这样可避免复杂的几何操作,从而大大提高程序计算效率和鲁棒性.

图6 裂尖单元分割形态Fig.6 The partitioning patterns of tip elements

5.3 非连续单元积分

为便于区分,以下各式中,以x表示全局坐标,ξ表示有限元母单元上的自然坐标,表示子三角形母单元上的局部坐标.

切割单元、裂尖单元非连续积分可分别表示为

式(32)和式(33)中,Nsub表示子剖分三角形个数,Ngp表示子三角形上高斯点个数,二者上标“+”、“–”、“t”分别表示裂纹面上、下和裂尖材料域;|J|为全局坐标与自然坐标之间雅可比矩阵行列式;ξkl为高斯点自然坐标,Wkl为高斯点权重,分别可表示为

6 裂纹相关计算

6.1 应力强度因子

相互作用积分法计算应力强度因子(SIF),具有精度高、适应范围广的优点.故本文采用该方法进行计算,其定义式如下[53]

相互作用积分与真实场和附加场应力强度因子之间存在以下关系

式中,E′为等效杨氏模量.对于平面应变问题:E′=E/(1 −ν2);对于平面应力问题:E′=E.

在实际计算中,常在围线Γ0外再增设一回路C(如图7),将积分项乘以一光滑权函数q(x),并将线积分转化为等效面积分[54]

其中,q(x)在Γ 内取1,在C0外取0.A为由Γ,C0,C+和C−围成的闭合区域;mj为A的单位外法向量.

图7 相互作用积分域定义Fig.7 Interaction integral domain

6.2 裂纹扩展

本文采用最大周向拉应力强度因子理论[55]计算裂纹扩展方向.令以下临界断裂韧度对θ 偏导为

可得裂纹扩展角度

其中,θ0∈(−π,π).当KII=0 时,取θ0=0.

7 算例分析

7.1 孔洞问题

如图8(a)所示,在一边长l=2 m 的方板中心有一半径为a=0.2 m 的圆孔.板受x向均匀拉伸作用,拉应力σ0=1 Pa.板的弹性模量E=1000 Pa,泊松比ν=0.3.计算网格如图8(b)所示,网格数为40×40.假设该问题为平面应变问题.在极坐标系下,该问题应力场解析解可表示为[33]

图8 方板中心带一圆孔:(a)问题描述;(b)计算网格Fig.8 A square plate with a circular hole at its center:(a)Problem definition;(b)simulation mesh

图9 σx 应力云图对比:本文(左侧);解析解(右侧)Fig.9 σx stress contour:Present(left);exact(right)

使用XFEM 联合本文非连续单元分割方法计算得到的应力云图,与解析解对比如图9 所示.从该图可看出,本文与解析解计算结果非常接近.特别是在孔洞周围,本文方法计算得到的应力集中因子为2.992,与理论值3.0 非常接近.进一步地,取y轴上a≤y≤L/2 区间范围内的应力,对比本文方法与解析解计算结果,如图10 所示.该图定量说明本文计算结果与解析解符合较好.该孔洞算例表明,本文发展的非连续单元分割方法在XFEM 求解孔洞问题上具有有效性.

图10 本文方法计算结果与解析解对比Fig.10 Comparison between present and exact solutions

为研究本文分割方法计算效率,在本算例中,同时使用本文发展的非连续单元分割方法以及德劳内三角剖分(Delaunay triangulation)方法,进行非连续单元分割.以下以Sub-DT 表示德劳内三角剖分方法;以Sub-LSV 表示本文非连续单元分割方法.

由于计算模型的对称性,本文仅给出1/4 孔洞与固体单元(共7 个)分割计算时间对比,如表2 所示.从该表可知,使用德劳内方法(Sub-DT)计算总时间为219µs,而使用本文方法(Sub-LSV)计算总时间仅为57 µs.该算例说明,本文非连续单元分割方法在计算效率上明显优于传统德劳内方法.

表2 使用不同分割方法所消耗计算时间Table 2 Computational costs with different partitioning methods

7.2 夹杂问题

如图11(a)所示,在一边长l=2 m 的方板中心有一半径为a=0.2 m 的圆孔.板受y向均匀压缩作用,压应力σ0=1 Pa.基质材料弹性模量E1=1000 Pa,泊松比ν1=0.3;夹杂材料弹性模量E2=3000 Pa,泊松比ν2=0.3.假设该问题为平面应变问题.

图11 方板中心带一夹杂:(a)问题描述;(b)有限元网格Fig.11 A square plate with a circular inclusion at its center:(a)Problem definition;(b)finite element mesh

XFEM 计算网格仍采用图8(b)中的40×40 规则网格.该问题无解析解,为验证方法可行性,同时使用有限元法(FEM)进行计算,网格如图11(b).该网格具有1685 个节点,3288 个三角形单元.

使用XFEM 联合本文非连续单元分割方法计算得到uy位移云图,与FEM 对比,如图12 所示.XFEM与FEM 计算的uy最大值均为1.809 mm,二者一致.取直线y=x上−L/2≤x≤L/2 区间范围内ux和uy位移,对比二者结果,如图13 所示,图中曲线两两重合.进一步分析可知,二者总位移平均偏差为0.04%.该夹杂算例表明,本文发展的非连续单元分割方法在XFEM 求解夹杂问题上同样具有有效性.

图12 uy 位移云图对比:本文(左侧);有限元(右侧)Fig.12 uy displacement contour:Present(left);FEM(right)

7.3 稳定裂纹问题

如图14(a)所示,一边裂纹板受x向剪切作用,剪应力τ0=1 Pa.板宽b=7 m、高2h=16 m,裂纹长a=3.5 m.板弹性模量E=1000 Pa,泊松比ν=0.3.计算网格如图14(b)所示,网格数为19×39.假设为平面应变问题.该问题中裂尖处应力强度因子的参考解为[56]

在本算例中,同时使用本文发展的非连续单元分割方法(Sub-LSV)以及德劳内三角剖分方法(Sub-DT),进行非连续单元分割.

图13 本文方法计算结果与有限元对比Fig.13 Comparison between present and FEM solutions

图14 边裂纹板受单向剪切:(a)问题描述;(b)计算网格Fig.14 Plate with an edge crack under unidirectional shear:(a)Problem definition;(b)simulation mesh

分别联合Sub-LSV,Sub-DT 和XFEM,IXFEM 计算裂尖处应力强度因子,如表3 所示.一方面,通过对比Sub-LSV 和Sub-DT 两种分割方法可知,不论是应用于XFEM 还是应用于IXFEM,Sub-LSV 均可取得接近于Sub-DT 的计算精度.Sub-LSV 精度略低的原因是:其分割裂尖单元为5 个子三角形,少于Sub-DT分割的6 个子三角形,导致精度略有降低.另一方面,通过对比XFEM 和IXFEM 可知,IXFEM 具有明显更高的计算精度.该稳定裂纹算例表明,本文发展的非连续单元分割方法在XFEM 和IXFEM 求解稳态裂纹问题上均具有有效性.

表3 使用不同方法计算得到的应力强度因子Table 3 Stress intensity factors from different approaches

7.4 扩展裂纹问题

如图15(a)所示,为一带边裂纹双臂梁模型,长l=6 m,宽w=2 m,裂纹长a=2.05 m.双臂梁弹性模量E=1000 Pa,泊松比ν=0.3.梁右端固支,左端受非对称集中力载荷P1=1 N,P2=1.01 N.计算网格如图15(b)所示,网格数为25×75.假设裂纹扩展18 步,每一步扩展长度为0.1 m.

图15 双臂梁带一边裂纹:(a)问题描述;(b)计算网格Fig.15 Double cantilever beam with an edge crack:(a)Problem definition;(b)simulation mesh

分别使用Sub-LSV,Sub-DT 联合IXFEM 计算裂纹扩展路径,如图16 所示.从该图可看出,使用IXFEM+Sub-LSV 计算得到的裂纹扩展路径,与使用IXFEM+Sub-DT 计算得到的裂纹扩展路径具有高度重合性.该算例进一步说明,对于扩展裂纹问题,本文发展的非连续单元分割方法同样具有有效性.

图16 裂纹扩展路径对比Fig.16 Comparison of crack propagation paths

7.5 三维裂纹问题

本文发展的非连续单元模板分割方法,可推广至三维情形.但由于三维情况下,弯曲裂纹切割单元需要进行等参变换等诸多算法研究,因此另文叙述.此处,预先给出一三维边裂纹算例,以供参考.

如图17(a)所示,为一三维贯穿型边裂纹板,该板宽w=7 m,高h=16 m,厚t=4 m,裂纹长a=3.5 m.板弹性模量E=1000 Pa,泊松比ν=0.0.板底端固支,顶端受y轴正向拉伸载荷σ0=1 Pa 作用.计算网格如图17(b)所示,网格数为19×39×10.

分别使用高阶高斯积分法和本文切割方法,并联合XFEM 进行求解.计算所得到的uy位移场对比,如图18 所示.使用XFEM+高阶高斯积分法,位移最大值为73.92 mm,最小值为−4.978 mm.使用XFEM+本文分割方法,位移最大值为73.91 mm,最小值为−4.976 mm.两种方法的位移最大值偏差为0.01%,位移最小值偏差为0.04%.该算例说明,对于三维裂纹问题,本文发展的非连续单元分割方法同样具有有效性.

图17 三维贯穿型边裂纹模型:(a)问题描述;(b)计算网格Fig.17 Three-dimensional cut-through edge crack model:(a)Problem definition;(b)simulation mesh

图18 uy 位移场对比:(a)高阶高斯积分法;(b)本文分割方法Fig.18 Comparison of uy displacement fields:(a)High-order Gauss integration method;(b)presented partitioning method

8 结论

本文提出一种基于单元水平集的模板分割方法,用于非连续单元子剖分和数值积分.首先,遍历单元水平集值所有形态并建立标准单元分割模板库;然后,根据单元水平集值,对非标准单元进行形态查询和模板插值;最后,套用标准单元分割模板实现单元高效分割和子剖分.

将该方法与常规XFEM、改进型XFEM 进行结合,从而应用于孔洞、夹杂、裂纹等非连续问题分析中.算例分析表明,本文提出的分割方法具有较高的计算精度.由于不引入复杂几何操作,该模板分割方法同时具有较高计算效率和鲁棒性,可为XFEM 类方法在非连续问题中提供有效支撑.

作为方法预研,本文仅给出二维非连续单元分割算法,用以探究模板分割可行性.事实上,笔者已初步将该算法推广至三维情形.在三维情况下,对于弯曲裂纹切割单元需要进行等参变换等诸多算法研究,所以将另文叙述.本文分割方法同时还可用于界面与计算网格不对齐的其他数值方法,如数值流形法(NMM)、有限格子法(FCM)等等.