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

doi:10.6052/0459-1879-20-360

非连续问题中单元分割的模板方法¹⁾

王理想^*^,^†, 文龙飞^*^,^†, 肖桂仲^*^,^†^,^**, 田荣^,^*^,^†^,²⁾

^*中物院高性能数值模拟软件中心, 北京 100088

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

^**南京理工大学, 南京 210094

A TEMPLATED METHOD FOR PARTITIONING OF SOLID ELEMENTS IN DISCONTINUOUS PROBLEMS¹⁾

Wang Lixiang^*^,^†, Wen Longfei^*^,^†, Xiao Guizhong^*^,^†^,^**, Tian Rong^,^*^,^†^,²⁾

^*CAEP Software Center for High Performance Numerical Simulation$,$ Beijing $100088$$,$ China

^†Institute of Applied Physics and Computational Mathematics$,$ Beijing $100088$$,$ China

^**Nanjing University of Science and Technology$,$ Nanjing $210094$$,$ China

通讯作者: 2) 田荣, 研究员, 主要研究方向: 计算力学与高性能计算. E-mail:tian_rong@iapcm.ac.cn

收稿日期: 2020-10-20 网络出版日期: 2021-03-18

基金资助:

1) 国家重点研发计划.  2016YFB0201002
国家重点研发计划.  2016YFB0201004
科学挑战专题.  TZ2018002

Received: 2020-10-20 Online: 2021-03-18

作者简介 About authors

摘要

扩展有限元法 (extended finite element method, XFEM) 因具有裂纹几何独立于模拟网格、裂纹扩展时无需网格重分重映、计算精度高等优点,成为裂纹分析的主流数值方法之一. 但该方法在工程实践中存在单元被裂纹分割的几何困难 —— 现有精确几何分割方法实现复杂、计算量大、鲁棒性差. 为克服这一困难, 本文提出一种基于单元水平集的模板分割方法, 用于非连续单元子剖分和数值积分. 首先, 遍历单元水平集值所有形态并建立标准单元分割模板库; 然后, 根据单元水平集值, 对非标准单元进行形态查询和模板插值; 最后, 套用标准单元分割模板实现单元高效分割和子剖分. 将该方法与常规XFEM、改进型XFEM进行结合,从而应用于孔洞、夹杂、裂纹等非连续问题分析中. 算例分析表明, 本文提出的模板分割方法具有较高计算精度. 由于不引入复杂几何操作, 该模板分割方法同时具有较高计算效率和鲁棒性, 故可为XFEM类方法在实际工程应用中提供有效支撑.

关键词： 子剖分 ; 水平集 ; 非连续性 ; 裂纹扩展 ; 常规XFEM ; 改进型XFEM

Abstract

The extended finite element method (XFEM) has been one of the privileged tools for crack analysis due to its significant advantages: (1) Independence of crack geometry on the simulation mesh; (2) no necessity of remeshing when a crack grows; and (3) high accuracy. However, the method is hindered in engineering practices by the partitioning difficulty of discontinuous elements, i.e. the geometric interaction between discontinuous interfaces and solid elements. Though current partitioning algorithms are geometrically exact, they are cumbersome to implement, computationally expensive, and insufficiently robust. To overcome these issues, a templated partitioning algorithm is proposed based on element level sets for subdivision and numerical integration of discontinuous elements. Firstly, a templated partitioning library for standard discontinuous elements is established by enumerating all the patterns of element level set values. Secondly, the pattern of a non-standard element to be partitioned is looked up and the sub-coordinates are interpolated based on the element level set values. Lastly, the non-standard element is efficiently partitioned into sub-triangles based on the standard element template. The algorithm is incorporated into the conventional XFEM and the improved XFEM for analysis of discontinuous problems, i.e. the problems with holes, inclusions, cracks and so forth. Numerical examples indicate that the proposed algorithm achieves favorable accuracy. Without cumbersome geometrical operations, the templated partitioning algorithm is also efficient and robust, thereby enabling itself to support the extended finite element methods in practical engineering problems.

Keywords： subdivision ; level set ; discontinuity ; crack propagation ; conventional XFEM ; improved XFEM

PDF (849KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

王理想, 文龙飞, 肖桂仲, 田荣. 非连续问题中单元分割的模板方法¹⁾. 力学学报[J], 2021, 53(3): 823-836 DOI:10.6052/0459-1879-20-360

Wang Lixiang, Wen Longfei, Xiao Guizhong, Tian Rong. A TEMPLATED METHOD FOR PARTITIONING OF SOLID ELEMENTS IN DISCONTINUOUS PROBLEMS¹⁾. Chinese Journal of Theoretical and Applied Mechanics[J], 2021, 53(3): 823-836 DOI:10.6052/0459-1879-20-360

引言

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

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

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

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

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

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

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

(1) 高阶高斯积分法^[34-35]: 在单元内布置大量高斯积分点 (如6$\times$7^[34], 8$\times$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所示非连续静力问题, 应满足平衡方程

(1)$\nabla \cdot \boldsymbol{\sigma}+\boldsymbol{b}=\mathbf{0}, \quad \text { in } \Omega$

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 二维非连续静力问题描述

Fig. 1 Illustration of statics with discontinuities in 2D

式中, $\nabla $为梯度算子, $\boldsymbol{\sigma }$为应力张量, ${b}$为体力.

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

(2)$\boldsymbol{u}=\overline{\boldsymbol{u}}, \text { on } \Gamma_{\mathrm{u}}$

(3)$\left.\begin{array}{c}\boldsymbol{\sigma} \cdot \boldsymbol{n}=\overline{\boldsymbol{t}}, \text { on } \Gamma_{\mathrm{t}} \\\boldsymbol{\sigma} \cdot \boldsymbol{n}=\mathbf{0}, \text { on } \Gamma_{0} \\\boldsymbol{\sigma} \cdot \boldsymbol{n}=\mathbf{0}, \text { on } \Gamma_{\mathrm{c}} \\\boldsymbol{\sigma} \cdot \boldsymbol{n}=\mathbf{0}, \text { on } \Gamma_{\mathrm{h}}\end{array}\right\}$

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

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

(4)$\boldsymbol{\sigma }={D}:\boldsymbol{\varepsilon}$

式中, ${D}$为弹性张量, $\boldsymbol{\varepsilon }$为应变张量.

小变形下, 应变$\boldsymbol{\varepsilon }$和位移${u}$满足如下几何关系

(5)$\boldsymbol{\varepsilon }=\nabla_{s} {u}=\frac{1}{2}\left({\nabla {u}+(\nabla {u})^{T}} \right)$

1.2 弱形式

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

(6)$\int_\varOmega \delta \boldsymbol{\varepsilon }:{\boldsymbol{\sigma }} {d}\varOmega -\int_\varOmega {\delta {u}\cdot {b}} {d}\varOmega-\int_{{\varGamma }_{t} } {\delta {u}\cdot \bar{t}} {d}{\varGamma }=0$

式中, $\delta {u}$为虚位移, $\delta \boldsymbol{\varepsilon }$为虚应变.

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

(7)$\int_{\varOmega} (\nabla_{s} \delta {u}):{D}:(\nabla _{s} {u}) {d}\varOmega-\int_{\varOmega} \delta {u}\cdot {b}{d}\varOmega -\\\int_{{\varGamma }_{t} } \delta {u}\cdot \bar{t} {d}{\varGamma} = 0$

2 非连续界面水平集描述

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

2.1 非连续面水平集函数$\varphi (x)$

非连续面水平集函数$\varphi(x)$定义为求解域$\varOmega$内任意给定一点$x\in \varOmega $到非连续面$\varGamma_{{d}}$的含符号距离, 其数学表达式为

(8)$\varphi (x)=\mathop{\min}\limits_{{\bar{{x}}}\in \varGamma_{{d}} } \left\| {x-{\bar{{x}}}} \right\|{sign}\left[ {{n}\cdot (x-{\bar{{x}}})} \right]$

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

2.2 裂纹尖端水平集函数$\psi (x)$

裂纹尖端水平集函数$\psi (x)$定义为求解域$\varOmega$内任意给定一点$x\in \varOmega $到裂纹尖端$x_{t}$的含符号距离, 其数学表达式为

(9)$\psi (x)=t\cdot (x-x_{t} )$

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

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 裂纹水平集函数构造 (根据文献[33]重画)

Fig. 2 Construction of level set functions for a crack (redrawn after Ref.^[33])

2.3 裂纹尖端极坐标系($r$, $\theta )$

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

(10)$r=r(\boldsymbol{x})=\sqrt{\varphi^{2}(\boldsymbol{x})+\psi^{2}(\boldsymbol{x})}$

(11)$\theta=\theta(\boldsymbol{x})=\arctan \left(\frac{\varphi(\boldsymbol{x})}{\psi(\boldsymbol{x})}\right)$

3 扩展有限元

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

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

(12)${u}^{h}(x)=\sum\limits_{i\in {\cal I}} {N_{i}(x)\bar{u}_{i} } +\\\sum\limits_{j\in {\cal J}} {N_{j} (x)\left( {F_{kink} (x)-F_{kink} (x_{j} )} \right){\bar{{a}}}_{j} }$

式中, ${\cal I}$为全部节点集合, ${\cal J}$为非连续面加强节点集合; $i$和$j$为序数; $\bar{u}$为常规自由度, ${\bar{{a}}}$为非连续面加强额外自由度; $N(x)$为标准有限元形函数, $F_{kink}(x)$为弱非连续加强函数,可取为^[50]

(13)$F_{kink} (x)=\sum\limits_{j\in {\cal J}} {N_{j} (x)\left| {\varphi_{j} } \right|} -\left| {\sum\limits_{j\in {\cal J}} {N_{j} (x)\varphi_{j} } } \right|$

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

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

(14)${u}^{h}(x)=\sum\limits_{i\in {\cal I}} {N_{i} (x)\bar{u}_{i} } +\\\sum\limits_{j\in {\cal J}} {N_{j} (x)\left( {F_{{jump}} (x)-F_{{jump}} (x_{j} )} \right){\bar{{b}}}_{j} } +\\\sum\limits_{k\in {\cal K}} {N_{k} (x)\sum\limits_{\alpha =1}^4 {\left( {F_{{tip}}^{\alpha } (x)-F_{{tip}}^{\alpha } (x_{k} )} \right){\bar{{c}}}_{k}^{\alpha } } }$

式中, ${\cal I}$为全部节点集合, ${\cal J}$为裂纹面加强节点集合, ${\cal K}$为裂尖加强节点集合; $i$, $j$, $k$, $\alpha $为序数;$\bar{u}$为常规自由度, ${\bar{{b}}}$为裂纹面加强额外自由度, ${\bar{{c}}}$为裂尖加强额外自由度; $N(x)$为标准有限元形函数, $F_{jump}(x)$为阶跃加强函数, $F_{tip}(x)$为裂尖奇异加强函数. $F_{jump}(x)$常使用如下Heaviside函数

(15)$F_{{jump}} (x)=H(\varphi (x))=\left\{ {{\begin{array}{l@{\quad }l} 1, & {{if\ \ }\varphi (x)\geqslant 0} \\ 0, & {{if\ \ }\varphi (x)<0} \\ \end{array} }} \right.$

对于各向同性线弹性材料, $F_{tip}(x)$可取为

(16)$\left\{ {F_{{tip}}^{\alpha } (x)} \right\}=\left\{ {F_{{tip}}^{\alpha } \left( {x(r,\theta )} \right)} \right\}= \Bigg\{ \sqrt r \cos \frac{\theta }{2}, \sqrt r \sin \frac{\theta }{2},\\\sqrt r \sin \frac{\theta }{2}\sin \theta ,\sqrt r \cos \frac{\theta }{2}\sin \theta \Bigg\}$

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

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

(17)${u}^{h}(x)=\sum\limits_{i\in {\cal I}/{\cal K}} {N_{i} (x)\bar{u}_{i} } +\\\sum\limits_{j\in {\cal J}} {N_{j} (x)\left( {F_{{jump}} (x)-F_{{jump}} (x_{j} )}\right){\bar{{b}}}_{j} } +\\\sum\limits_{k\in {\cal K}} {\left( {\sum\limits_{m\in {\cal K}_{k} }{N_{m} (x)\phi_{k}^{m} (x)} } \right)\bar{u}_{k} }$

式中, ${\cal I}$为所有节点集合, ${\cal J}$为裂纹面加强节点集合, ${\cal K}$为裂尖加强节点集合, ${\cal K}_{k} \subset {\cal K}$为节点$k$处裂尖加强影响域内节点集合; $i$, $j$, $k$, $m$为序数; $\bar{u}$为常规自由度, ${\bar{{b}}}$为裂纹面加强额外自由度; $N(x)$为标准有限元形函数, $F_{jump}(x)$为阶跃加强函数, 与式(15)定义相同, $\phi_{k}^{m} (x)$为裂尖局部加强函数

(18)$\phi_{k}^{m} ({x})={p}^{T}({x})\left( {{A}^{-1}{p}_{k} -\frac{{A}_{(1)}^{-1} {A}_{(1)}^{-{T}} {p}_{k} }{{A}_{11}^{-1} }+\frac{{A}_{(1)}^{-1} \delta _{mk} }{{A}_{11}^{-1} }} \right)$

其中, $p(x)$为基向量, $p_{k}$为$p(x_{k})$的缩写; $A$为矩量矩阵, ${A}_{(1)}^{-1}$为${A}^{-1}$的第一列, ${A}_{11}^{-1} $为${A}_{(1)}^{-1}$的第一个元素; $\delta $为Kronecker符号. $p(x)$和$A$的表达式如下所示

(19)$\boldsymbol{p}(\boldsymbol{x})=\left[1, \frac{\boldsymbol{x}^{\mathrm{T}}-\boldsymbol{x}_{k}^{\mathrm{T}}}{c h}, F_{\mathrm{tip}}^{\alpha}(\boldsymbol{x})-F_{\mathrm{tip}}^{\alpha}\left(\boldsymbol{x}_{k}\right)\right]^{\mathrm{T}}$

(20)$\boldsymbol{A}=\sum_{l \in \mathcal{K}_{k}} \boldsymbol{p}_{l} \boldsymbol{p}_{l}^{\mathrm{T}}$

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

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

3.4 混合单元修正

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

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 混合单元处理: 处理之前 (左); 处理之后 (右)

Fig. 3 Blending element treatment: Before (left); after (right)

引入如下斜坡函数^[52]

(21)$R({x})=\sum\limits_{{x}\in \varOmega_{{bld}} } {N({x})}$

则在混合单元$\varOmega_{{bld}} $上, 位移逼近具有如下加权形式

(22)${u}^{h}({x})=\left( {1-R({x})} \right){u}_{{std}}^{h} ({x})+R({x}){u}_{{enr}}^{h} ({x})$

式中, ${u}_{{std}}^{h} ({x})$, ${u}_{{enr}}^{h}({x})$分别为标准有限元位移逼近和加强单元位移逼近. 对于上式具体形式,XFEM和IXFEM可分别参考文献[51]和文献[26,27,28,29,30,31].

4 控制方程离散

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

(23)${KU}={F}$

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

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

(24)$\left[ {{\begin{array}{c@{\quad }c} {{K}_{uu} } & {{K}_{ua} } \\ {{K}_{au} } & {{K}_{aa} } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {\bar{u}} \\ {{\bar{{a}}}} \\ \end{array} }} \right\}=\left\{ {{\begin{array}{*{20}c} {{F}_{u} } \\ {{F}_{a} } \\ \end{array} }} \right\}$

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

(25)$\left[ {{\begin{array}{c@{\quad }c@{\quad }c} {{K}_{uu} } & {{K}_{ub} } & {{K}_{uc} } \\ {{K}_{bu} } & {{K}_{bb} } & {{K}_{bc} } \\ {{K}_{cu} } & {{K}_{cb} } & {{K}_{cc} } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {\bar{u}} \\ {{\bar{{b}}}} \\ {{\bar{{c}}}} \\ \end{array} }} \right\}=\left\{ {{\begin{array}{*{20}c} {{F}_{u} } \\ {{F}_{b} } \\ {{F}_{c} } \\ \end{array} }} \right\}$

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

(26)$\left[ {{\begin{array}{c@{\quad }c} {{K}_{uu} } & {{K}_{ub} } \\ {{K}_{bu} } & {{K}_{bb} } \\ \end{array} }} \right]\left\{ {{\begin{array}{*{20}c} {\bar{u}} \\ {{\bar{{b}}}} \\ \end{array} }} \right\}=\left\{ {{\begin{array}{*{20}c} {{F}_{u} } \\ {{F}_{b} } \\ \end{array} }} \right\}$

式(24)$\sim \!$式(26)中${K}$和${F}$可统一表示为

(27)$\boldsymbol{K}_{i j}^{\alpha \beta}=\int_{\Omega}\left(\boldsymbol{B}_{i}^{\alpha}\right)^{\mathrm{T}} \boldsymbol{D} \boldsymbol{B}_{j}^{\beta} \mathrm{d} \Omega$

(28)$\boldsymbol{F}_{i}^{\alpha}=\int_{\Gamma_{\mathrm{t}}}\left(\boldsymbol{N}_{i}^{\alpha}\right)^{\mathrm{T}} \overline{\boldsymbol{t}} \mathrm{d} \Gamma+\int_{\Omega}\left(\boldsymbol{N}_{i}^{\alpha}\right)^{\mathrm{T}} \boldsymbol{b} \mathrm{d} \Omega$

其中, $\alpha $, $\beta = 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

新窗口打开| 下载CSV

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 三角形切割单元分割形态

Fig. 4 The partitioning patterns of triangular cut elements

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

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

(29)$\boldsymbol{\xi}_{p} -\boldsymbol{\xi}_{i} =\lambda (\boldsymbol{\xi}_{j} -\boldsymbol{\xi}_{i} )$

其中, $\lambda $可由水平集值插值得到

(30)$\lambda =\frac{\varphi_{p} -\varphi_{i} }{\varphi_{j} -\varphi_{i} }{=}\frac{-\varphi_{i} }{\varphi_{j} -\varphi_{i} }$

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

(31)$\boldsymbol{\xi}_{p} =\frac{\varphi_{j} \boldsymbol{\xi}_{i} -\varphi_{i} \boldsymbol{\xi}_{j} }{\varphi_{j} -\varphi_{i} }$

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

算法1 三角形切割单元子三角化

输入: 三角形切割单元水平集值

输出: 三角形切割单元子三角化

1: if type $==$ 1 and sign $==$ a

2: $\quad \varOmega_{e}^{+} =\{\Delta ijk\}, \quad \varOmega_{e}^{-} =\{\emptyset \}$

3: if type $==$ 1 and sign $==$ b

4: $\quad \varOmega_{e}^{+} =\{\emptyset \}, \quad \varOmega_{e}^{-} =\{\Delta ijk\}$

5: if type $==$ 2 and sign $==$ a

6: $\quad \varOmega_{e}^{+} =\{\Delta ijq,\Delta iqp\}, \quad \varOmega_{e}^{-} =\{\Delta pqk\}$

7: if type $==$ 2 and sign $==$ b

8: $\quad \varOmega_{e}^{+} =\{\Delta pqk\}, \quad \varOmega_{e}^{-} =\{\Delta ijq,\Delta iqp\}$

9: if type $==$ 3 and sign $==$ a

10: $\quad \varOmega_{e}^{+} =\{\Delta ijp\}, \quad \varOmega_{e}^{-} =\{\Delta ipk\}$

11: if type $==$ 3 and sign $==$ b

12: $\quad \varOmega_{e}^{+} =\{\Delta ipk\}, \quad \varOmega_{e}^{-} =\{\Delta ijp\}$

5.1.2 四边形单元

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

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

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 四边形切割单元分割形态

Fig. 5 Partitioning patterns of quadrilateral cut elements

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

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

算法2 四边形切割单元子三角化

输入: 四边形切割单元水平集值

输出: 四边形切割单元子三角化

1: if type $==$ 1 and sign $==$ a

2: $\quad \varOmega_{e}^{+} =\{\Box ijkl\}, \quad \varOmega_{e}^{-} =\{\emptyset \}$

3: if type $==$ 1 and sign $==$ b

4: $\quad \varOmega_{e}^{+} =\{\emptyset \}, \quad \varOmega_{e}^{-} =\{\Box ijkl\}$

5: if type $==$ 2 and sign $==$ a

6: $\quad \varOmega_{e}^{+} =\{\Delta ijp,\Delta pjq,\Delta qjk\}, \quad \varOmega_{e}^{-} =\{\Delta pql\}$

7: if type $==$ 2 and sign $==$ b

8: $\quad \varOmega_{e}^{+} =\{\Delta pql\}, \quad \varOmega_{e}^{-} =\{\Delta ijp,\Delta pjq,\Delta qjk\}$

9: if type $==$ 3 and sign $==$ a

10: $\quad \varOmega_{e}^{+} =\{\Delta ijp,\Delta pjk\}, \quad \varOmega_{e}^{-} =\{\Delta pkl\}$

11: if type $==$ 3 and sign $==$ b

12: $\quad \varOmega_{e}^{+} =\{\Delta pkl\}, \quad \varOmega_{e}^{-} =\{\Delta ijp,\Delta pjk\}$

13: if type $==$ 4 and sign $==$ a

14: $\quad \varOmega_{e}^{+} =\{\Delta ijk\}, \quad \varOmega_{e}^{-} =\{\Delta ikl\}$

15: if type $==$ 4 and sign $==$ b

16: $\quad \varOmega_{e}^{+} =\{\Delta ikl\}, \quad \varOmega_{e}^{-} =\{\Delta ijk\}$

17: if type $==$ 5 and sign $==$ a

18: $\quad \varOmega_{e}^{+} =\{\Delta ijp,\Delta pjq\}, \quad \varOmega_{e}^{-} =\{\Delta pql,\Delta lqk\}$

19: if type $==$ 5 and sign $==$ b

20: $\quad \varOmega_{e}^{+} =\{\Delta pql,\Delta lqk\}, \quad \varOmega_{e}^{-} =\{\Delta ijp,\Delta pjq\}$

5.2 裂尖单元

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

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 裂尖单元分割形态

Fig. 6 The partitioning patterns of tip elements

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

算法3 裂尖单元子三角化

输入: 裂尖单元水平集值、裂尖坐标

输出: 裂尖单元子三角化

1: if type $==$ 1 and elem $==$ triangle

2: $\quad \varOmega_{e}^{{tip}} =\{\Delta itp,\Delta ijt,\Delta tjk,\Delta ptk\}$

3: if type $==$ 2 and elem $==$ triangle

4: $\quad \varOmega_{e}^{{tip}} =\{\Delta ijt,\Delta tjk,\Delta kit\}$

5: if type $==$ 1 and elem $==$ quadrilateral

6: $\quad \varOmega_{e}^{{tip}} =\{\Delta itp,\Delta ijt,\Delta tjk,\Delta tkl,\Delta ptl\}$

7: if type $==$ 2 and elem $==$ quadrilateral

8: $\quad \varOmega_{e}^{{tip}} =\{\Delta ijt,\Delta tjk,\Delta tkl,\Delta itl\}$

5.3 非连续单元积分

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

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

(32)$\begin{aligned}\boldsymbol{K}_{i j}^{\alpha \beta} &=\int_{\Omega_{e}^{\text {cut }}(x)}\left[\boldsymbol{B}_{i}^{\alpha}(\boldsymbol{x})\right]^{\mathrm{T}} \boldsymbol{D}\left[\boldsymbol{B}_{j}^{\beta}(\boldsymbol{x})\right] \mathrm{d} \Omega_{(\boldsymbol{x})}=\\& \int_{\Omega_{e}^{+}(\xi)}\left[\boldsymbol{B}_{i}^{\alpha}(\xi)\right]^{\mathrm{T}} \boldsymbol{D}\left[\boldsymbol{B}_{j}^{\beta}(\xi)\right]|\boldsymbol{J}(\xi)| \mathrm{d} \Omega_{(\xi)}+\\& \int_{\Omega_{e}^{-}(\xi)}\left[\boldsymbol{B}_{i}^{\alpha}(\xi)\right]^{\mathrm{T}} \boldsymbol{D}\left[\boldsymbol{B}_{j}^{\beta}(\xi)\right]|\boldsymbol{J}(\xi)| \mathrm{d} \Omega_{(\xi)}=\\& \sum_{\mathrm{sub}}^{+} \sum_{\mathrm{g}}^{+}\left[\boldsymbol{B}_{i}^{\alpha}\left(\xi_{k l}\right)\right]^{\mathrm{T}} \boldsymbol{D}\left[\boldsymbol{B}_{j}^{\beta}\left(\xi_{k l}\right)\right]\left|\boldsymbol{J}\left(\xi_{k l}\right)\right| W_{k l}+\\& \sum_{k=1}^{\alpha} \\& \sum_{k=1}^{\mathcal{S u b}} \sum_{l=1}^{-}\left[\boldsymbol{B}_{i}^{\alpha}\left(\xi_{k l}\right)\right]^{-} \boldsymbol{D}\left[\boldsymbol{B}_{j}^{\beta}\left(\xi_{k l}\right)\right]\left|\boldsymbol{J}\left(\xi_{k l}\right)\right| W_{k l}\end{aligned}$

(33)$\begin{aligned}\boldsymbol{K}_{i j}^{\alpha \beta}=& \int_{\Omega_{e}^{\mathrm{ip}}(\boldsymbol{x})}\left[\boldsymbol{B}_{i}^{\alpha}(\boldsymbol{x})\right]^{\mathrm{T}} \boldsymbol{D}\left[\boldsymbol{B}_{j}^{\beta}(\boldsymbol{x})\right] \mathrm{d} \Omega_{(\boldsymbol{x})}=\\& \int_{\Omega_{e}^{\text {tip }}(\xi)}\left[\boldsymbol{B}_{i}^{\alpha}(\xi)\right]^{\mathrm{T}} \boldsymbol{D}\left[\boldsymbol{B}_{j}^{\beta}(\xi)\right]|\boldsymbol{J}(\xi)| \mathrm{d} \Omega_{(\xi)}=\\& \sum_{k=1}^{\mathrm{S}} \sum_{l=1}^{\mathrm{sub}}\left[\boldsymbol{B}_{i}^{\alpha}\left(\xi_{k l}\right)\right]^{\mathrm{T}} \boldsymbol{D}\left[\boldsymbol{B}_{j}^{\beta}\left(\xi_{k l}\right)\right]\left|\boldsymbol{J}\left(\xi_{k l}\right)\right| W_{k l}\end{aligned}$

式(32)和式(33)中, ${\cal N}_{{sub}} $表示子剖分三角形个数, ${\cal N}_{{gp}} $表示子三角形上高斯点个数,二者上标“$+$”、“--”、“t”分别表示裂纹面上、下和裂尖材料域; $\left| J \right|$为全局坐标与自然坐标之间雅可比矩阵行列式; $\boldsymbol{\xi}_{kl}$为高斯点自然坐标, $W_{kl} $为高斯点权重, 分别可表示为

(34)$\boldsymbol{\xi}_{k l}=\sum_{n=1}^{3} \bar{N}_{n}\left(\overline{\boldsymbol{\xi}}_{k l}\right) \boldsymbol{\xi}_{n}^{\mathrm{sub}}$

(35)$W_{k l}=\left(2 A_{k l}\right) w_{k l}$

式(34)和式(35)中, $\bar{N}(\overline{\boldsymbol{\xi}})$为子三角形形函数, $\overline {\xi}_{kl}$为子三角形上高斯点局部坐标, $\boldsymbol{s}_{n}^{\mathrm{sub}}$为子三角形顶点自然坐标, $w_{kl} $为子三角形上积分点权重, $A_{kl} $为 (自然坐标下) 子三角形面积.

6.1 应力强度因子

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

(36)$I^{(1,2)}=\int_\varGamma \left( {W^{(1,2)}\delta _{1j} -\sigma _{ij}^{(1)} \frac{\partial u_{i}^{(2)} }{\partial x_{1} }-\sigma_{ij}^{(2)} \frac{\partial u_{i}^{(1)} }{\partial x_{1} }} \right)n_{j} \ {d}\varGamma$

式中, $\Big(\sigma_{ij}^{(1)}$, $\varepsilon_{ij}^{(1)}$, $u_{i}^{(1)}\Big)$, $\Big(\sigma_{ij}^{(2)}$, $\varepsilon_{ij}^{(2)}$, $ u_{i}^{(2)}\Big)$分别为真实场、辅助场, $W^{(1,2)}=\sigma_{ij}^{(1)}\varepsilon_{ij}^{(2)} =\sigma_{ij}^{(2)} \varepsilon_{ij}^{(1)}$为相互作用应变能密度, $\varGamma$为包含裂尖的围线, $n_{j}$为$\varGamma$单位外法向量.

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

(37)$I^{(1,2)}=\frac{2}{{E}'}\left( {K_{{I}}^{(1)} K_{{I}}^{(2)} +K_{{II}}^{(1)} K_{{II}}^{(2)} } \right)$

式中, ${E}'$为等效杨氏模量. 对于平面应变问题: ${E}'={E/{(1-\nu^{2})}}$;对于平面应力问题: ${E}'=E$.

在式(37)中分别令$K_{{I}}^{(2)} =1$, $K_{{II}}^{(2)}=0$和$K_{{I}}^{(2)} =0$, $K_{{II}}^{(2)} =1$,可得$K_{{I}}^{(1)} =I^{(1,2)}_{mode\ I}{{{E}'}/2}$, $K_{{II}}^{(1)}=I^{(1,2)}_{mode\ II}{{{E}'}/2}$.

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

(38)$I^{(1,2)}=\int_C {\left[ {W^{(1,2)}\delta _{1j} -\sigma_{ij}^{(1)} \frac{\partial u_{i}^{(2)} }{\partial x_{1} }-\sigma_{ij}^{(2)} \frac{\partial u_{i}^{(1)} }{\partial x_{1} }} \right]qm_{j}{d}C} =\\ \int_A {\left[ {W^{(1,2)}\delta _{1j} -\sigma_{ij}^{(1)} \frac{\partial u_{i}^{(2)} }{\partial x_{1} }-\sigma_{ij}^{(2)} \frac{\partial u_{i}^{(1)} }{\partial x_{1} }} \right]q_{,j}{d}A}$

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 相互作用积分域定义

Fig. 7 Interaction integral domain

其中, $q({x})$在$\varGamma$内取1, 在$C_0$外取0.$A$为由$\varGamma$, $C_0$, $C_{+}$和$C_{-}$围成的闭合区域; $m_{j}$为$A$的单位外法向量.

6.2 裂纹扩展

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

(39)$K_{\theta {C}} =\cos \frac{\theta }{2}\left( {K_{{I}} \cos ^{2}\frac{\theta }{2}-\frac{3K_{{II}} }{2}\sin \theta } \right)$

可得裂纹扩展角度

(40)$\theta_{0} =2\arctan \frac{1}{4}\left[ {\frac{K_{{I}} }{K_{{II}} }-{sign}(K_{{II}} )\sqrt {\left( {\frac{K_{{I}} }{K_{{II}} }} \right)^{2}+8} } \right]$

其中, $\theta_{0} \in (-\pi ,\pi )$. 当$K_{{II}} =0$时, 取$\theta_{0}=0$.

7 算例分析

7.1 孔洞问题

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

(41)$\left. \begin{array}{l} \sigma_{r} =\dfrac{\sigma_{0} }{2}\left( {1-\dfrac{a^{2}}{r^{2}}} \right)+\dfrac{\sigma_{0} }{2}\left( {1-\dfrac{a^{2}}{r^{2}}} \right)\left( {1-3\dfrac{a^{2}}{r^{2}}} \right)\cos 2\theta \\[5mm] \sigma_{\theta } =\dfrac{\sigma_{0} }{2}\left( {1+\dfrac{a^{2}}{r^{2}}} \right)-\dfrac{\sigma_{0} }{2}\left( {1+3\dfrac{a^{4}}{r^{4}}} \right)\cos 2\theta \\[5mm] \tau_{r\theta } =-\dfrac{\sigma_{0} }{2}\left( {1-\dfrac{a^{2}}{r^{2}}} \right)\left( {1+3\dfrac{a^{2}}{r^{2}}} \right)\sin 2\theta \\ \end{array} \right\}$

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 方板中心带一圆孔: (a) 问题描述; (b) 计算网格

Fig. 8 A square plate with a circular hole at its center: (a) Problem definition; (b) simulation mesh

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

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 $\sigma_{x}$应力云图对比: 本文 (左侧); 解析解 (右侧)

Fig. 9 $\sigma_{x}$ stress contour: Present (left); exact (right)

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 本文方法计算结果与解析解对比

Fig. 10 Comparison between present and exact solutions

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

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

表2 使用不同分割方法所消耗计算时间

Table 2 Computational costs with different partitioning methods

新窗口打开| 下载CSV

7.2 夹杂问题

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

图11

新窗口打开| 下载原图ZIP| 生成PPT

图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$\times$40规则网格. 该问题无解析解,为验证方法可行性, 同时使用有限元法 (FEM) 进行计算, 网格如图11(b).该网格具有1685个节点, 3288个三角形单元.

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

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 $u_{y}$位移云图对比: 本文 (左侧); 有限元 (右侧)

Fig. 12 $u_{y}$ displacement contour: Present (left); FEM (right)

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13 本文方法计算结果与有限元对比

Fig. 13 Comparison between present and FEM solutions

7.3 稳定裂纹问题

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

(42)$K_{{I}}^{{ref}} =34.0{Pa}\cdot {m}^{{1/2}},\ \ K_{{II}}^{{ref}} =4.55 \ {Pa}\cdot {m}^{{1/2}}$

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14 边裂纹板受单向剪切: (a) 问题描述; (b) 计算网格

Fig. 14 Plate with an edge crack under unidirectional shear: (a) Problem definition; (b) simulation mesh

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

分别联合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

新窗口打开| 下载CSV

7.4 扩展裂纹问题

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

图15

新窗口打开| 下载原图ZIP| 生成PPT

图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

新窗口打开| 下载原图ZIP| 生成PPT

图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, 泊松比$\nu =0.0$.板底端固支, 顶端受$y$轴正向拉伸载荷$\sigma_{0} =1$ Pa作用.计算网格如图17(b)所示, 网格数为19$\times$39$\times$10.

图17

新窗口打开| 下载原图ZIP| 生成PPT

图17 三维贯穿型边裂纹模型: (a) 问题描述; (b) 计算网格

Fig. 17 Three-dimensional cut-through edge crack model: (a) Problem definition; (b) simulation mesh

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

图18

新窗口打开| 下载原图ZIP| 生成PPT

图18 $u_{y}$位移场对比: (a) 高阶高斯积分法; (b) 本文分割方法

Fig. 18 Comparison of $u_{y}$ displacement fields: (a) High-order Gauss integration method; (b) presented partitioning method

8 结论

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

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

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

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

王勃, 张阳博, 左宏, 等.

压应力对压剪裂纹扩展的影响研究

力学学报, 2019,51(3):845-851