基于图像四叉树的改进型比例边界有限元法研究1)

孙立国 江守燕 杜成斌

(河海大学工程力学系,南京 210098)

引言

比例边界有限元法(scaled boundary finite element method,SBFEM)具有仅需离散结构域边界使得计算问题的维度降低一维、可存在悬挂节点而不会带来额外的计算问题等优点,该方法已广泛应用于结构的计算分析中[1-2],特别是在断裂损伤问题领域[3-11],SBFEM 具有其独特的优势,可自动地解析求出缝端的应力强度因子.

由于SBFEM 本身的固有优势,SBFE 网格可为任意的多边形单元[12],特别地,SBFEM 与四叉树网格的结合[13-17]使得前处理工作具有快速、高效的特点,粗细网格过渡十分方便,并能实现自动化的网格剖分,大大地减轻了网格剖分的负担.然而,上述这些研究在模拟界面演化问题(如裂纹扩展问题)时,虽然能够使得网格重剖分工作达到最小化,但局部的网格重剖分工作仍不可避免.

Natarajan 等[18]和大连理工大学李建波等[19-23]结合扩展有限元法和SBFEM 的优点,提出了扩展比例边界有限元法(X-SBFEM)的概念,该方法采用扩展有限单元模拟裂纹贯穿单元、SBFE 模拟裂尖单元,从而避免了传统扩展有限元法需要构造裂尖单元复杂改进函数以及需要采用特殊的数值积分技术进行裂尖单元刚度矩阵求解的弊端.江守燕等[24-25]进一步完善和改进了该方法,借助于扩展有限元法的基本思想,将扩展有限元法中所有改进单元转化为SBFE 的子域,被裂纹切割的改进单元无需引入特殊的富集函数,取得了较好的模拟效果.Huynh 等[26]通过在扩展有限元法中耦合SBFEM 的基本思想,缓解了裂纹生长模拟中网格细化过程的悬挂节点问题.

本工作在文献[24-25]的研究基础上,结合四叉树网格的优势,提出基于图像四叉树的改进型比例边界有限元法,通过结构域图像实现结构域网格的全自动剖分,且结构内不连续处局部范围的网格粗细过渡方便、快速.

1 基于图像的四叉树网格

四叉树生成的网格可以直接转换为SBFEM 子域,而不需要像传统有限元那样区分网格的角节点和悬挂节点.由于SBFEM 的半解析性质,它只需对结构边界进行离散化,径向具有解析解,只需比例中心对所有边界都有良好的可视性,即使边界大小不同,也不会影响数值精度.因此,四叉树网格可以完美地用作SBFEM 子域,基于图像的四叉树网格剖分技术可实现SBFEM 网格的快速、全自动剖分[27].下面具体介绍基于图像的四叉树网格剖分的关键步骤.

(1)获取待四叉树网格剖分的结构图像.图1(a)显示了一个像素为N×M的图像,图像内部包含一个夹杂,基体、夹杂以及两者之间的界面通过不同颜色区分.使用四叉树算法对其进行网格剖分,图像颜色矩阵QN×M用于存储像素的所有颜色.

(2)初始化四叉树网格.如图1(b)所示,通过将图像分割为多个边长为smax的单元来初始化四叉树.

(3)如果单元中最大和最小颜色强度之间的差异大于颜色阈值,则该单元被递归地进一步划分为四个大小相等的单元.该过程递归执行,直到所有单元满足均匀性标准或达到最小边缘长度.然后得到最终的四叉树网格,见图1(c).

在四叉树分解过程中使用了平衡四叉树分解,该分解仅产生16 种可能的单元模式,通过考虑旋转对称性,只有六种类型的唯一节点排列是可能的,见图1(d).在SBFEM 中,这些悬挂节点可以直接地视为新多边形单元的节点,无需任何特殊处理.

图1 基于图像的四叉树网格剖分Fig.1 Image-based quadtree cell generation

2 改进型比例边界有限元法

改进型SBFEM 的主要思想是: 在SBFEM 的理论框架下,基于水平集法表征材料的内部裂纹面,含不连续裂纹面的子域可通过节点水平集函数识别,裂纹扩展时无需进行网格重剖分,界面的几何特征通过SBFEM 子域的附加自由度表征.

2.1 裂纹的水平集描述

φ(x)Φk(x)(k=1,2)

如图2 所示,两个水平集和用于描述裂纹[28].通常假设裂纹尖端水平集函数Φk(x)(k=1,2)与裂纹面水平集函数 φ(x) 正交.函数φ(x)可以用符号距离函数表示为

图2 裂纹的水平集表征Fig.2 Crack description based on two level set functions

式中,x是考察点的坐标;x*是考察点在裂纹表面上的投影点;n是垂直于裂纹表面的外法向向量;sign(x)是符号函数.函数 Φk(x)(k=1,2) 可以定义为

式中,xk是第k个裂纹尖端的坐标,t是第k个裂纹尖端的单位切向矢量.

2.2 改进子域的识别

基于水平集函数,可以确定四叉树网格与裂纹的切割关系.如图3 所示,可以观察到两种类型的关系: (1)四叉树单元被裂纹完全切割,即域 Ωcc;(2)四叉树单元被裂纹部分切割,单元包含一个裂纹尖端,即域 Ωcp.改进子域识别的具体流程如下.

图3 四叉树网格与裂纹间的切割关系Fig.3 Relationship between the quadtree cell and a crack

步骤1: 遍历求解域内所有节点,根据式(1)和式(2)求出所有节点的两个水平集函数值.

步骤2: 针对某一个含有N个节点的四叉树单元,根据式(2)计算的每个节点的波前水平集函数值,判断裂尖点是否位于该单元内.若该单元包含裂尖点,则属于域 Ωcp.若不包含裂尖点,则转至步骤3 进一步判断该单元是否属于域 Ωcc.

步骤3: 若单元属于域 Ωcc,需要满足条件: 其中M(0 ＜M＜N)个节点的裂纹面水平集函数 φ(x)＜0,则剩余的 (N-M) 个节点的裂纹面水平集函数φ(x)≥0 .若不满足前述条件,则该单元为常规单元.

步骤4: 对于求解域内的所有四叉树单元,重复步骤2、步骤3,即可识别出所有改进单元.

2.3 改进的比例边界有限元子域

如图4(a) 所示,完全由裂纹切割的单元视为SBFEM 的两个独立子域.与常规SBFEM 不同的是,通过在裂纹与单元边界的交点处引入额外的节点,将原四叉树单元一分为二.在本研究中,SBFE 子域将通过在交叉点处引入额外的自由度而得到富集,单元并没有物理分区,这样就无需改变网格以遵从几何界面,尤其是当裂纹扩展时无需进行局部的网格重剖分.

如图4(b)所示,对于包含裂纹尖端的单元,选择裂纹尖端周围的若干单元形成超级单元.根据裂尖所在单元找出裂尖单元周边的单元,如周边单元无法形成一个方形域,则继续寻找,直至包围裂尖的单元可形成一个方形域为止,超级单元作为一个SBFEM子域.超级单元内部的节点位移可通过SBFE 位移场近似获得.被裂纹切割的单元边界并没有物理分开,而是增加了额外的自由度来表示不连续性.

图4 改进的SBFE 子域Fig.4 Enriched SBFE subdomain

2.4 改进子域刚度矩阵的求解

如图5 所示,具有4 个角节点和2 个悬挂节点的四叉树单元(SBFE 子域)的边被划分为多个线单元,即L1,L2,L3,L4,L5,L6.单元局部坐标系ξ-η的 ξ 方向由比例中心指向子域边界,η 方向沿着子域环向边界.比例中心O位于四边形单元的中心,对所有边界均具有较好的可见性.值得注意的是,与有限元不同,在SBFEM 中,悬挂节点不会带来额外的计算困难,它只是将一条线单元分成两条.

图5 比例边界有限元的四叉树网格Fig.5 SBFE quadtree cell

SBFEM 位移模式可以表示为

式中,N(η) 为形函数矩阵,u(ξ) 为径向的位移函数.

由位移模式式(3),可得应变场为

线性算子矩阵L可表示为

式中,|J| 为雅克比矩阵行列式,可表示为

根据胡克定律,应力场可表示为

式中,D为本构矩阵且

根据虚功原理,可得SBFE 方程为

载荷向量为

系数矩阵 (Ei,i=0,1,2) 为

为求解式(12),首先变换该式为一个一阶微分方程

式中,哈密顿矩阵Z为

通过Schur 分解,可得

其中,Schur 矩阵S和变换V可表达为

式中,Sn和Sp为与哈密顿矩阵Z的负特征值和正特征值对应的上三角矩阵;Vu和分别为有限域和无限域的模态位移;Vq和分别为与有限域和无限域相应的模态力.对于有限域问题,仅Sn在比例中心处产生有限位移,u(ξ) 可以表示为

式中,c为依赖于边界条件的积分常数.根据节点位移ub=u(ξ=1),积分常数可表达为

SBFE 子域的边界节点力为

因此,SBFE 子域的刚度矩阵可通过下式求解得到

如图6(a)所示,四叉树单元被裂纹完全切割,切割面为位移不连续界面,在局部坐标系下,该单元包含四个角节点(节点编号: 1,2,3,4)和两个悬挂节点(节点编号: 5,6),由于被裂纹切割,该单元还包括四个虚节点(节点编号: 7,8,9,10),将四叉树单元(子域)分割成两个次子域SS1 和SS2.在两个次子域SS1 和SS2 内,分别采用式(25)计算这两个次子域的刚度矩阵和,然后集成四叉树单元(子域)的刚度矩阵,其形式如下

图6 改进单元节点编号示意图(*表示比例中心)Fig.6 Schematic diagram of node number of enriched element(* means scale center)

式中,kuu表示与真实节点(即角节点和悬挂节点)相应的刚度矩阵贡献;kaa表示与虚节点(裂纹与四叉树单元边界的交叉点)相应的刚度矩阵贡献;kua和kau为真实节点与虚节点的耦合贡献.

如图6(b)所示,四叉树单元被裂纹面部分切割,取围绕裂尖一圈的若干单元组成一个超级单元,超级单元的比例中心取为裂尖位置,该超级单元包括两个虚节点(节点编号: 17,18),采用式(25)计算超级单元的刚度矩阵.

2.5 改进型比例边界有限元法的控制方程

由2.4 节的讨论可知,无论何种改进类型的单元,单元刚度矩阵均具有式(26)的形式,因此,可得到改进型SBFEM 的控制方程为

约78%的学生表示遇到过此类困难。随机抽样访谈中，部分受访学生反映曾因文化差异在涉外场合导致过沟通障碍，也有为此闹过笑话的经历。

式中,u为常规节点(角节点和悬挂节点)的位移未知量,a为附加的虚节点的位移未知量.使用这种方法引入SBFEM 的附加节点未知量a可以解释为裂纹与单元边界交点处的实际位移.

3 裂纹生长准则

在SBFEM 中,I 型和II 型应力强度因子KI和KII可根据裂纹尖端 θ=0o处的奇异应力和求得,即

在包含裂纹尖端的超级单元中,当r→0 时,式(20)包含对角块,奇异应力可表达为

图7 为裂尖局部坐标系,通过插值可以得到裂纹尖端 θ=0o处的奇异应力模态值.奇异应力模态(η(θ))在沿裂纹子域边界的高斯积分点使用式(30)计算. θ=0o时的应力模态可以通过插值得到.本研究中使用样条插值来获得更精确的应力.

图7 裂尖局部坐标系Fig.7 Local coordinate system at the crack tip

根据关系式 ξ=r/L0,应力强度因子可通过下式计算[3]最大周向应力准则用于确定裂纹扩展方向,该准则给出了以下裂纹扩展方向

等效应力强度因子Ke可通过下式计算

4 数值算例

4.1 中心斜裂纹问题

本算例为复合型应力强度因子的求解问题,如图8 所示,双向拉伸板中心含一倾斜裂纹,裂纹倾角为 θ,裂纹半长度为a,该板在无限远处受单轴拉伸载荷 σ 作用.该问题I 和II 型应力强度因子的解析解为[29]

图8 含中心斜裂纹的拉伸板Fig.8 Tension plate with an inclined centre crack

数值计算时,取拉伸板边长L=16 m,裂纹半长度a=0.5 m,拉伸载荷 σ=1.0 kPa,按平面应变问题考虑,为消除刚体位移,板的下边缘竖向约束,板的左下角水平向约束,如图9 所示,四叉树网格单元总数为1276,节点总数为1365,弹性模量取为1.0 MPa,泊松比取为0.2,这里取的弹性模量和泊松比数值仅为了验证论文建议方法计算复合型应力强度因子的正确性,弹性模量和泊松比的取值不同对计算结果影响很小.

图9 四叉树网格Fig.9 Quad-tree cells

图10 给出了采用建议的改进型SBFEM 计算得到的裂纹不同倾斜角度下拉伸板的I 型和II 型应力强度因子值,并给出该问题的解析解作为比较,从图中可看出,数值解与解析解吻合较好,从而验证了数值方法的有效性.四叉树网格可实现网格的快速过渡,在含有裂纹的局部区域可采用较密的网格.

图10 拉伸板应力强度因子随倾角变化曲线Fig.10 Stress intensity factors of the tensile plate with different crack inclination angles

4.2 含预制裂纹的重力坝模型

如图11,某混凝土重力坝比例模型含1 条初始裂纹,试件的尺寸和边界条件如图所示.数值计算时,该重力坝模型被离散成四叉树网格,包括4140 个单元和4392 个节点,按平面应变问题考虑.重力坝模型的材料参数取值参照文献[30],弹性模量E=35.7 GPa、泊松比 ν=0.1、密度 ρ=2400 kg/m3、断裂能G=0.184 N/mm、临界应力强度因子KIC=81.456 N·mm1/2、裂纹生长长度取为150 mm.

重力坝模型所受载荷包括自重和水载荷.水载荷的模拟简化成在图11 中四个加载点处施加一定大小的集中力,施加的集中力大小比例如图所示.初始裂纹的位置如图11 所示,深度为150 mm.图12给出了裂纹的最终扩展路径,同时与文献[30]的数值和试验结果进行比较,二者一致性较好.

图11 含初始裂纹的重力坝模型(单位: mm)Fig.11 Gravity dam model with initial crack (unit: mm)

图12 重力坝模型的裂纹开裂扩展路径Fig.12 Crack propagation path of gravity dam model

4.3 含双裂纹双孔洞的薄板模型

如图13 所示,一尺寸为 20 mm×10 mm 带有半径r=2 mm 的双孔洞薄板受双向集中拉伸载荷作用,板的左、右边缘处分别有长a=1 mm 的裂纹.数值计算时,该薄板模型被离散成四叉树网格,包括2456 个单元和3000 个节点.按平面应力问题考虑,模型的材料参数取值参照文献[13],弹性模量E=200 GPa、泊松比ν=0.3、密度ρ=7800 kg/m3、临界应力强度因子KIC=1500 N·mm1/2、裂纹生长长度取为0.5 mm.

图13 双孔板结构内多裂纹扩展问题(单位: mm)Fig.13 Multiple crack propagation in a thin plate with two holes(unit: mm)

图14 给出了预测的多裂纹扩展路径,当开始施加外载荷后,两条裂纹都向邻近的孔洞位置扩展,随后裂纹逐渐改变扩展角度,以几乎水平的方向向薄板中部扩展,当两裂纹都扩展到试件中部以后,由于裂尖的应力集中效应,两条裂纹彼此相向扩展,模拟结果与文献[13]采用的自适应网格的SBFEM 数值结果进行了比较,二者一致性较好.

图14 含双孔板问题的预测多裂纹开裂扩展路径Fig.14 Predicted multiple crack paths for the plate with two holes

4.4 L 形混凝土面板

Winkler 等[31]对L 形混凝土面板进行了实验测试,许多学者已经对其进行了模拟,以验证数值模型.面板的几何尺寸如图15 所示.面板承受垂直集中载荷F.混凝土面板的细观结构,使用具有代表性的直径为45 mm、32.5 mm 和17.5 mm 的三段式圆形骨料级配,骨料面积分数为40%.L 形面板的细观结构如图16 所示.要生成四叉树网格,具体图像的大小为1024 × 1024 像素,最小边缘长度设置为4 × 4 像素.L 形面板的材料特性见表1.初始裂纹生长长度取为5 mm.

表1 L 型板的材料参数Table 1 Material properties for the L-shaped panel

图15 L 型板的示意图(单位: mm)Fig.15 Schematic diagram of an L-shaped panel (unit: mm)

图16 L 型板的细观结构Fig.16 Mesostructure of L-shaped panel

图17 给出了当前数值结果计算得到的裂纹开裂路径并与试验结果进行了对比,图18 给出了数值计算得到的载荷-加载点位移曲线并与试验结果进行了对比,从图中可以看出,当前建议的基于四叉树网格和改进型比例边界有限元法的模拟结果均与试验结果吻合较好.

图17 开裂路径Fig.17 Crack path

图18 载荷-加载点位移曲线Fig.18 Load-displacement curves at loading point

5 结论

SBFEM 是一种半解析的数值方法,由于其本身的固有优势,可以方便使用悬挂节点而无需任何特殊处理,且SBFEM 可以在求解过程中自动解析地求出裂尖的应力强度因子,SBFEM 与四叉树可完美结合并实现网格的全自动剖分.论文提出的改进型SBFEM 继承了传统SBFEM 的所有优点,通过引入虚节点的思想,将裂纹与四叉树单元边界交叉点作为虚节点,虚节点的自由度作为附加自由度处理,并采用水平集函数表征裂纹面,在模拟裂纹扩展问题时无需进行裂尖的局部网格重剖分工作.

最后,通过三个数值算例验证了建议方法的数值精度,提出的基于四叉树网格的改进型SBFEM可准确求出裂纹的复合型应力强度因子;且能够较好地模拟出裂纹的开裂扩展路径,裂纹扩展过程无需进行网格重剖分.