文章快速检索     高级检索
  浙江大学学报(工学版)  2017, Vol. 51 Issue (7): 1324-1330  DOI:10.3785/j.issn.1008-973X.2017.07.008
0

引用本文 [复制中英文]

郁杨天, 章青, 顾鑫. 近场动力学与有限单元法的混合模型与隐式求解格式[J]. 浙江大学学报(工学版), 2017, 51(7): 1324-1330.
dx.doi.org/10.3785/j.issn.1008-973X.2017.07.008
[复制中文]
YU Yang-tian, ZHANG Qing, GU Xin. Hybrid model of peridynamics and finite element method under implicit schemes[J]. Journal of Zhejiang University(Engineering Science), 2017, 51(7): 1324-1330.
dx.doi.org/10.3785/j.issn.1008-973X.2017.07.008
[复制英文]

基金项目

国家自然科学基金资助项目(11372099,11132003);江苏省自然科学基金资助项目(BK20151493)

作者简介

郁杨天(1992—), 男, 硕士生, 从事灾变破坏力学的研究.
orcid/org/0000-0002-6397-4158.
Email: yuyangtian@hhu.edu.cn

通信联系人

章青, 男, 教授.
orcid/org/0000-0002-2819-5976.
Email: lxzhangqing@hhu.edu.cn

文章历史

收稿日期:2016-04-09
近场动力学与有限单元法的混合模型与隐式求解格式
郁杨天 , 章青 , 顾鑫     
河海大学 工程力学系, 江苏 南京 210098
摘要: 利用近场动力学方法(PD)在模拟不连续变形问题的独特优势和有限单元法(FEM)的计算效率,提出近场动力学与有限单元法混合建模的方法,并用于求解断裂力学问题.裂纹出现的区域采用改进的近场动力学微观弹脆性(PMB)模型进行离散,其他区域采用有限元离散,通过杆单元连接不同的子区域.在隐式求解体系下实现了两种方法的混合建模,该模型在求解静力学问题时无需引入阻尼项,有效提高了计算效率和计算精度.通过模拟计算简支梁的弹性变形和三点弯曲梁Ⅰ型裂纹的扩展过程,与理论解吻合良好,验证了提出的混合模型和求解方法的准确性和有效性.
关键词: 近场动力学    有限单元法(FEM)    混合模型    改进的近场动力学微观弹脆性模型    隐式求解    
Hybrid model of peridynamics and finite element method under implicit schemes
YU Yang-tian , ZHANG Qing , GU Xin     
Department of Engineering Mechanics, Hohai University, Nanjing 210098, China
Abstract: A hybrid model of peridynamics (PD) and finite element method (FEM) was proposed and applied to solve problems of fracture mechanics in order to combine the unique advantage of PD in solving discontinuities and the computational efficiency of FEM. The improved prototype microelastic brittle (PMB) model of peridynamics was utilized for the regions where material failure was expected. The region without failure was discretized by FEM. The truss element was introduced to bridge peridynamic subregions and finite element subregions. The hybrid model is based on the implicit schemes, and it need not consider a fictitious damping term in solving static problems. The computational efficiency and accuracy of the model were improved. The static elastic deformation of a simply supported beam and the propagation process of mode Ⅰ fracture in a three points bend beam were simulated to verify the accuracy and utility of the presented model. Results obtained by the model agreed well with the theoretical solutions.
Key words: peridynamics    finite element method (FEM)    hybrid model    improved prototype microelastic brittle model    implicit scheme    

有限单元法已在工程中得到了广泛应用, 但在求解裂纹扩展等不连续问题时, 通常须借助断裂准则以判断开裂位置及扩展方向, 在裂纹扩展后要不断进行网格重构, 计算结果具有网格依赖性[1].设置黏结单元[2-3]虽然可以避免计算过程的复杂性, 但必须预先判断裂纹扩展方向和所在区域.有限单元法在求解断裂破坏问题时所面临的局限性源自于连续性假设的理论基础与实际问题不连续的矛盾.近场动力学[4-6]作为一种新兴的基于非局部积分思想的数值方法, 避免了传统的局部微分方程求解不连续问题时的奇异性, 在破坏问题分析中具有独特的优势, 取得了许多新的成果[7-14], 已成为国际学术界研究的热点, 在岩石破裂和土体破坏分析等方面得到了成功应用[15-19].

近场动力学作为一种粒子类方法, 计算效率较低, 为了充分利用有限元与近场动力学各自的优势, 国外一些学者对近场动力学与有限元的混合建模方法进行研究.Macek等[20]将近场动力学中的物质点对用有限元中的杆单元来描述, 采用有限元分析软件ABAQUS中的镶嵌单元(embedded element), 实现了近场动力学与有限单元法的混合建模;为了避免镶嵌单元刚度过大, 镶嵌单元的弹性模量和密度必须取很小的数值.Kilic等[21]通过设置重叠区域来实现近场动力学与有限元的耦合, 在重叠区内, 采用有限元结点位移插值的方法确定近场动力学物质点的位移, 但无法通过物质点的位移求得有限元的结点位移.Liu等[22]通过引入界面单元实现了两种方法的耦合, 给出两种耦合方案, 但在界面单元处, 由近场动力学物质点的受力情况计算界面单元结点力的过程较繁琐.此外, 上述混合建模方法均建立在显式求解体系下, 在求解静力学问题时必须引入阻尼项[21], 而阻尼的大小直接影响计算的收敛速度[23], 致使计算效率不高.

本文建立近场动力学与有限单元法的混合模型, 并应用于断裂力学问题.对可能出现破坏的区域, 采用近场动力学方法离散, 其他区域采用有限单元法离散.通过杆单元连接PD和FEM子域, 在隐式求解体系下实现了两种方法的混合建模.该模型在求解静力学问题时无需引入阻尼项, 提高了计算效率.此外, 在PD子域, 采用改进的微观弹脆性(PMB)模型, 提高了计算精度.基于所建立的混合模型, 模拟计算了简支梁的弹性变形和三点弯曲梁Ⅰ型裂纹的扩展过程, 验证了该方法的准确性和有效性.

1 理论简介 1.1 近场动力学方法

图 1所示, 假设在某一时刻t, 空间域R中任一物质点x与其邻近一定范围内的其他物质点(x′∈R:$\left\| \mathit{\boldsymbol{{x}'}}-\mathit{\boldsymbol{x}} \right\|$≤δ)存在相互作用力f.根据牛顿第二定律, 可得

图 1 物质点间的相互作用 Fig. 1 Interaction between material points
$ \begin{array}{l} \rho \mathit{\boldsymbol{\ddot u}}\left( {\mathit{\boldsymbol{x}},t} \right) = \int_{{{\rm{H}}_x}} {f\left( {\mathit{\boldsymbol{u}}\left( {\mathit{\boldsymbol{x'}},\mathit{\boldsymbol{t}}} \right) - \mathit{\boldsymbol{u}}\left( {\mathit{\boldsymbol{x}},\mathit{\boldsymbol{t}}} \right),\mathit{\boldsymbol{x'}} - \mathit{\boldsymbol{x}}} \right){\rm{d}}{V_{x'}}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;b\left( {\mathit{\boldsymbol{x}},t} \right). \end{array} $ (1)

式中:Hx为物质点x的近场范围, δ为近场范围尺寸, ρ为物质点的材料密度, u为物质点的位移, b为外荷载密度.

物质点间的相互作用力函数f称为本构力函数.它包含了材料的本构信息, 一般根据材料特性(如弹性、黏弹性、弹塑性、均匀各向同性、各向异性等)与长程作用力空间分布特征构造本构力函数[24-26].微观弹脆性(prototype microelastic brittle, PMB)材料[27]的本构力函数f

$ f\left( {\mathit{\boldsymbol{\eta }},\mathit{\boldsymbol{\xi }}} \right) = \frac{{\mathit{\boldsymbol{\xi + \eta }}}}{{\left\| {\mathit{\boldsymbol{\xi + \eta }}} \right\|}}\mathit{\boldsymbol{c}}\left( {\mathit{\boldsymbol{\xi }},\delta } \right)s\mu \left( {t,\mathit{\boldsymbol{\xi }}} \right). $ (2)

式中:ξ为相对位置, ξ=x′-xη为相对位移, η=u′-u; c(ξ, δ)为微观模量函数;s为物质点对的相对伸长率,

$ s = \frac{{\left\| {\mathit{\boldsymbol{\xi + \eta }}} \right\| - \left\| \mathit{\boldsymbol{\xi }} \right\|}}{{\left\| \mathit{\boldsymbol{\xi }} \right\|}}; $ (3)

μ为一标量函数,

$ \mu \left( {t,\mathit{\boldsymbol{\xi }}} \right) = \left\{ \begin{array}{l} 1,\;\;\;\;\;s\left( {t',\mathit{\boldsymbol{\xi }}} \right) < {s_0},0 \le t' \le t;\\ 0,\;\;\;\;其他. \end{array} \right. $ (4)

其中, s0为临界伸长率, 可以通过材料的抗拉强度ft和弹性模量E定义,

$ {s_0} = \frac{{{f_{\rm{t}}}}}{E}. $ (5)

s超过s0时, 点对间不再发生作用.

在近场动力学理论中, 定义标量函数φ(x, t)以反映物质点x的损伤:

$ \varphi \left( {\mathit{\boldsymbol{x}},t} \right) = 1 - \frac{{\int_H {\mu \left( {\mathit{\boldsymbol{x}},t,\mathit{\boldsymbol{\xi }}} \right){\rm{d}}{V_\xi }} }}{{\int_H {{\rm{d}}{V_\xi }} }}. $ (6)

式中:0≤φ≤1, φ=0表示材料未损伤, φ=1表示该点完全损伤, 不再与其他点发生相互作用.

1.2 改进的PMB模型

式(2) 中的微观模量函数c(ξ, δ)可以表示为

$ c\left( {\mathit{\boldsymbol{\xi }},\delta } \right) = c\left( {0,\delta } \right)g\left( {\mathit{\boldsymbol{\xi }},\delta } \right). $ (7)

式中:g(ξ, δ)为核函数, 反映物质点对间长程力的强度随两点间距离变化的规律.

在以往的PMB模型中, 核函数通常取

$ g\left( {\mathit{\boldsymbol{\xi }},\delta } \right) = \left\{ \begin{array}{l} 1,\;\;\;\;\left\| \mathit{\boldsymbol{\xi }} \right\| \le \delta ;\\ 0,\;\;\;\;\left\| \mathit{\boldsymbol{\xi }} \right\| > \delta . \end{array} \right. $ (8)

将物质点对的微观模量设为常数, 忽略物质点对间距离对微观模量函数的影响, 从而导致应用PMB模型计算弹性变形时存在较大的误差.

在改进的PMB模型[23, 28]中, 取核函数的表达式为

$ g\left( {\mathit{\boldsymbol{\xi }},\delta } \right) = \left\{ \begin{array}{l} {\left( {1 - {{\left( {\frac{{\left\| \mathit{\boldsymbol{\xi }} \right\|}}{\delta }} \right)}^2}} \right)^2},\;\;\;\;\;\;\left\| \mathit{\boldsymbol{\xi }} \right\| \le \delta ;\\ \;\;\;\;\;\;\;\;\;0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\| \mathit{\boldsymbol{\xi }} \right\| > \delta . \end{array} \right. $ (9)

该核函数能够反映长程力的空间分布规律, 计算精度高.根据近场动力学应变能密度与连续介质力学应变能密度相等的原则, 可以导出改进的PMB模型中c(0, δ)的表达式:

$ c\left( {0,\delta } \right) = \left\{ \begin{array}{l} \frac{{72E}}{{{\rm{ \mathsf{ π} }}{\delta ^4}}},\;\;\;\;\;\;\;三维;\\ \frac{{315E}}{{{\rm{8 \mathsf{ π} }}{\delta ^3}}},\;\;\;\;\;\;平面应力;\\ \frac{{210E}}{{{\rm{5 \mathsf{ π} }}{\delta ^3}}},\;\;\;\;\;\;平面应变. \end{array} \right. $ (10)
1.3 有限元方程

有限元分析的支配方程为

$ \mathit{\boldsymbol{Ku = F}}. $ (11)

式中:K为整体劲度矩阵, u为整体位移列阵, F为整体等效结点荷载列阵,

$ \left. \begin{array}{l} \mathit{\boldsymbol{K = }}\sum\limits_e {\mathit{\boldsymbol{C}}_e^{\rm{T}}\mathit{\boldsymbol{k}}{\mathit{\boldsymbol{C}}_e}} ,\\ \mathit{\boldsymbol{F = }}\sum\limits_e {\mathit{\boldsymbol{C}}_e^{\rm{T}}{\mathit{\boldsymbol{F}}^e}} . \end{array} \right\} $ (12)

其中, Ce为选择矩阵; k为单元劲度矩阵, $\mathit{\boldsymbol{k = }}\int_{{V^e}} {{\mathit{\boldsymbol{B}}^{\rm{T}}}} \mathit{\boldsymbol{DB}}{\rm{d}}V$, 其中B为应变转化矩阵, D为弹性矩阵;Fe为单元等效结点荷载列阵, ${\mathit{\boldsymbol{F}}^e} = \int_{{V^e}} {{\mathit{\boldsymbol{N}}^{\rm{T}}}} \mathit{\boldsymbol{r}}{\rm{d}}V\mathit{\boldsymbol{ + }}\int_{{S^e}} {{\mathit{\boldsymbol{N}}^{\rm{T}}}\mathit{\boldsymbol{\bar r}}{\rm{d}}S} $, 其中N为形函数矩阵, r为体积力, r为物体表面的面力.

2 数值计算格式 2.1 近场动力学隐式求解格式

将物体均匀离散, 取物质点间距为|Δx|, 运动方程(1) 的离散形式可以表示为

$ \rho \mathit{\boldsymbol{\ddot u}}_i^n = \sum\limits_j {\mathit{\boldsymbol{f}}\left( {\mathit{\boldsymbol{u}}_j^n - \mathit{\boldsymbol{u}}_i^n,{\mathit{\boldsymbol{x}}_j} - {\mathit{\boldsymbol{x}}_i}} \right){V_j} + \mathit{\boldsymbol{b}}_i^n} . $ (13)

式中:n为时间步长编号;Vjj处物质点体积, 在三维情况下Vj=|Δx|3, 在二维情况下Vj=|Δx|2.

对于静力问题, 令${\mathit{\boldsymbol{\ddot u}}}$=0, 则近场动力学平衡方程的离散形式为

$ \sum\limits_j {\mathit{\boldsymbol{f}}\left( {\mathit{\boldsymbol{u}}_j^n - \mathit{\boldsymbol{u}}_i^n,{\mathit{\boldsymbol{x}}_j} - {\mathit{\boldsymbol{x}}_i}} \right){V_j} + \mathit{\boldsymbol{b}}_i^n} = 0. $ (14)

式中:n表示第n级荷载.

将PMB模型的本构力函数式(2) 写成矩阵形式:

$ \mathit{\boldsymbol{f}} = {\mathit{\boldsymbol{K}}_{\rm{P}}}\mathit{\boldsymbol{u}} \cdot \mu \left( s \right). $ (15)

式中:KP为微观模量函数的矩阵形式.

建立局部直角坐标系x′y′z′, 并使x′轴与物质点对的长度方向平行, 如图 2所示.

图 2 局部坐标系下物质点对间的相互作用 Fig. 2 Interactions between material points in local coordinates

在局部坐标系下, 物质点对i、j的作用力和位移分别为

$ \begin{array}{l} \mathit{\boldsymbol{f'}} = {\left[ {{f_{x'i}},\;\;\;{f_{y'i}},\;\;{f_{z'i}},\;\;{f_{x'j}},\;\;\;{f_{y'j}},\;\;{f_{z'j}}} \right]^{\rm{T}}},\\ \mathit{\boldsymbol{u' = }}{\left[ {{{u'}_i},\;\;{{v'}_i},\;\;{{w'}_i},\;\;{{u'}_j},\;\;{{v'}_j},\;\;{{w'}_j}} \right]^{\rm{T}}}. \end{array} $ (16)

存在如下的关系:

$ \mathit{\boldsymbol{f'}} = {{\mathit{\boldsymbol{K'}}}_{\rm{P}}}\mathit{\boldsymbol{u'}} \cdot \mu \left( s \right). $ (17)

式中:

$ {{\mathit{\boldsymbol{K'}}}_{\rm{P}}} = \frac{c}{{\left\| \mathit{\boldsymbol{\xi }} \right\|}}\left[ {\begin{array}{*{20}{c}} 1&0&0&{ - 1}&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ { - 1}&0&0&1&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0 \end{array}} \right], $ (18)

对局部坐标系下的物质点对的微观模量矩阵Kp′进行坐标转换.令物质点对的长度方向(x′轴)与整体坐标xyz轴的余弦分别为lmn, 如图 3所示.

图 3 整体坐标系下的物质点对 Fig. 3 Material points in global coordinates

图 3可知,

$ \left. \begin{array}{l} l = \cos \alpha = \frac{{{x_j} - {x_i}}}{{\left\| \mathit{\boldsymbol{\xi }} \right\|}},\\ m = \cos \beta = \frac{{{y_j} - {y_i}}}{{\left\| \mathit{\boldsymbol{\xi }} \right\|}},\\ n = \cos \gamma = \frac{{{z_j} - {z_i}}}{{\left\| \mathit{\boldsymbol{\xi }} \right\|}}. \end{array} \right\} $ (19)

通过坐标转换, 可得整体坐标系下物质点对微观模量函数的矩阵形式:

$ {\mathit{\boldsymbol{K}}_{\rm{P}}} = \frac{c}{{\left\| \mathit{\boldsymbol{\xi }} \right\|}}\left[ {\begin{array}{*{20}{c}} {{l^2}}&{}&{}&{}&{}&{}\\ {lm}&{{m^2}}&{}&{}&{{\rm{sym}}}&{}\\ {\mathit{ln}}&{mn}&{{n^2}}&{}&{}&{}\\ { - {l^2}}&{ - lm}&{ - \mathit{ln}}&{{l^2}}&{}&{}\\ { - lm}&{ - {m^2}}&{ - mn}&{lm}&{{m^2}}&{}\\ { - \mathit{ln}}&{ - mn}&{ - {n^2}}&{\mathit{ln}}&{mn}&{{n^2}} \end{array}} \right]. $ (20)

于是, 式(14) 可以写成

$ \sum\limits_j {{\mathit{\boldsymbol{K}}_{\rm{P}}}\left( {\left| {{\mathit{\boldsymbol{x}}_j} - {\mathit{\boldsymbol{x}}_i}} \right|} \right)\left( {\mathit{\boldsymbol{u}}_j^n - \mathit{\boldsymbol{u}}_i^n} \right){V_j} + \mathit{\boldsymbol{b}}_i^n = {\bf{0}}.} $ (21)

可得与有限元方法求解静力问题形式相同的近场动力学求解方程.

2.2 近场动力学与有限元的混合模型

图 4所示, 将所考虑的物体划分为近场动力学子域和有限元子域进行离散.

图 4 有限元(FE)和近场动力学(PD)子区域 Fig. 4 Finite element (FE) and peridynamic (PD) regions

在两种区域的交界面上, 采用杆单元连接近场动力学中的物质点与有限元中的结点, 实现近场动力学与有限单元法的混合建模.如图 5所示, 交界面上的有限元结点不仅与所在单元的其他结点发生作用, 而且通过杆单元与其近场范围δ内的物质点相互作用.

图 5 结点与物质点连接示意图 Fig. 5 Interactions between FE nodes and PD points through trusses

参照式(9)、(10), 杆单元的劲度可以取为

$ k\left( {\xi ',\delta } \right) = \left\{ \begin{array}{l} \frac{{72E}}{{{\rm{ \mathsf{ π} }}{\delta ^4}}}{\left( {1 - {{\left( {\frac{{\xi '}}{\delta }} \right)}^2}} \right)^2},\;\;\;\;\;\;\;三维;\\ \frac{{315E}}{{{\rm{8 \mathsf{ π} }}{\delta ^3}}}{\left( {1 - {{\left( {\frac{{\xi '}}{\delta }} \right)}^2}} \right)^2},\;\;\;\;\;\;平面应力;\\ \frac{{210E}}{{{\rm{5 \mathsf{ π} }}{\delta ^3}}}{\left( {1 - {{\left( {\frac{{\xi '}}{\delta }} \right)}^2}} \right)^2},\;\;\;\;\;\;平面应变. \end{array} \right. $ (22)

式中:ξ′为有限元结点与物质点的相对距离.

在整体坐标系下, 杆单元劲度矩阵的形式与式(20) 相同.具体计算时, 只需用杆单元劲度k和有限元结点与物质点的相对距离ξ′替换式(20) 中的微观模量c和物质点对间距离‖ξ‖, 再形成总体的整体劲度矩阵和荷载列阵进行求解.

针对上述计算方法, 采用Fortran语言编写了计算程序, 主要的计算步骤如图 6所示.

图 6 PD与FEM混合建模计算流程图 Fig. 6 Flowchart for hybrid model of PD and FEM
3 算例分析 3.1 简支梁的弹性变形

考虑图 7所示的简支梁, 在梁的上表面中点受集中力P=10 kN的作用, 材料的弹性模量为30 GPa, 泊松比为0.3.在计算中, 将梁划分为两个有限元子域和一个近场动力学子域, 取物质点间距Δx=2.5 mm, 近场范围δ=4Δx=10 mm, 有限元网格采用10 mm×10 mm的四结点矩形单元.

图 7 简支梁模型示意图 Fig. 7 Geometry of specimen

为了分析比较, 对简支梁全部采用近场动力学模型和有限元模型进行计算.其中有限元采用ANSYS软件, 单元类型为SOLID四边形单元, 尺寸为10 mm×10 mm.

表 1所示为采用不同方法得到的中性轴跨中挠度计算结果.计算机型号为Intel Core i3 CPU @2.53 GHz, 内存为4.0 GB.表中,hm为跨中挠度,er为相对误差,tc为计算耗时.图 8给出利用混合模型计算得到的简支梁中性轴竖向位移h与解析解、有限单元法、近场动力学法计算结果的对比.图中,L为中性轴位置.可见, 采用混合模型方法计算得到的中性轴跨中挠度相对误差为0.31%, 与有限元解的相对误差0.14%相当, 并且计算耗时仅为近场动力学用时的一半, 验证了提出的近场动力学与有限元混合模型的正确性和有效性.

表 1 跨中挠度计算结果 Table 1 Results of mid-span deflection
图 8 中性轴竖向位移曲线 Fig. 8 Vertical displacement curves of neutral axis

在弹性问题的计算中, 虽然近场动力学方法的计算精度和计算效率均低于有限元法, 但近场动力学方法主要是用于有限元难以求解的材料和结构的破坏问题.采用提出的近场动力学与有限元混合模型, 精度与有限元相当, 计算效率得到了很大的提高, 为近场动力学方法应用于工程结构破坏问题的计算分析提供了一种有效的途径.

3.2 含初始裂缝三点弯曲梁的破坏分析

考虑图 9所示的三点弯曲梁, 跨中位置设置一预制裂缝, 初始裂缝长度a=80 mm, 材料的弹性模量为30 GPa, 泊松比为0.3, 抗拉强度为1.65 MPa.计算中, 将梁划分为两个有限元子域和一个近场动力学子域, 取Δx=2.5 mm, δ=4Δx=10 mm, 有限元单元采用大小为10 mm×10 mm的四结点矩形单元, 加载方式为位移加载.

图 9 三点弯曲梁几何尺寸 Fig. 9 Geometry of specimen

三点弯曲梁Ⅰ型裂纹张开口位移(crack-mouth-opening displacement, CMOD)的线弹性断裂力学(LEFM)的解答[30-31]

$ {\rm{CMOD = }}\frac{{4\sigma a}}{E} \cdot {V_1}\left( \alpha \right). $ (23)

式中:

$ \begin{array}{*{20}{c}} {{V_1}\left( \alpha \right) = 0.76 - 2.28\alpha + 3.87{\alpha ^2} - 2.04{\alpha ^3} + \frac{{0.66}}{{{{\left( {1 - \alpha } \right)}^2}}},}\\ {\alpha = a/D,\sigma = \frac{{3PL}}{{2{D^2}}}.} \end{array} $

图 10给出计算得到的荷载-CMOD曲线, 与线弹性断裂力学的结果有较好的一致性.图 11给出不同时刻梁的破坏特征.从计算结果可以发现, 荷载达到6.4 kN时结构出现损伤, 认为裂纹开始稳定扩展;荷载继续增加, 当荷载达到7.5 kN时, 裂纹开始失稳扩展.随着裂纹的不断扩展, 所需荷载逐渐减小.

图 10 荷载-CMOD曲线 Fig. 10 Curves of load-CMOD
图 11 梁破坏过程的模拟结果 Fig. 11 Failure process simulation of beam
4 结语

本文提出近场动力学与有限元混合建模的方法求解材料和结构的破坏问题.该方法在可能出现破坏的区域, 采用改进的近场动力学PMB模型离散, 其他区域采用有限元模型离散, 通过杆单元来连接各子区域.在隐式求解体系下实现了近场动力学与有限元的混合建模和数值求解.在求解静力学问题时, 无需引入阻尼项, 提高了计算效率和计算精度.简支梁的弹性变形和三点弯曲梁Ⅰ型裂纹的扩展过程的计算结果表明, 采用提出的PD-FEM混合模型和静力求解方法能够模拟裂纹扩展问题, 减少近场动力学的计算域, 节约计算时间, 为解决工程结构破坏问题提供了一种有效的分析方法.

近场动力学突破了传统连续介质力学方法求解破坏问题时的奇异性和困难, 在各类材料和结构的静动力裂纹扩展、冲击破坏等许多力学问题中均得到成功应用.近场动力学处于理论体系完善和初步应用阶段, 还有很大的研究空间, 在基于状态的近场动力学理论和方法、复杂环境下本构力函数的构建、非均匀离散技术、高效的数值求解体系、软件研制和工程应用等方面亟需加大研究力度.

参考文献
[1] MUROTANI K, YAGAWA G, CHOI J B. Adaptivefinite elements using hierarchical mesh and its application to crack propagation analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2013, 253: 1–14. DOI:10.1016/j.cma.2012.07.024
[2] CAMACHO G T, ORTIZ M. Computational modeling of impact damage in brittle materials[J]. International Journal of Solids and Structures, 1996, 33(20): 2899–2938.
[3] XU X P, NEEDLEMAN A. Numerical simulations of fast crack growth in brittle solids[J]. Journal of the Mechanics and Physics of Solids, 1994, 42(9): 1397–1434. DOI:10.1016/0022-5096(94)90003-5
[4] SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces[J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209. DOI:10.1016/S0022-5096(99)00029-0
[5] SILLING S A, EPTON M, WECKNER O, et al. Peridynamic states and constitutive modeling[J]. Journal of Elasticity, 2007, 88(2): 151–184. DOI:10.1007/s10659-007-9125-1
[6] SILLING S A. Linearized theory of peridynamic states[J]. Journal of Elasticity, 2010, 99(1): 85–111. DOI:10.1007/s10659-009-9234-0
[7] SILLING S A, BOBARU F. Peridynamic modeling of membranes and fibers[J]. International Journal of Nonlinear Mechanics, 2005, 40(2): 395–409.
[8] GERSTLE W, SAU N, SILLING S. Peridynamic modeling of concrete structures[J]. Nuclear Engineering and Design, 2007, 237(12): 1250–1258.
[9] ASKARI E, BOBARU F, LEHOUCQ R B, et al. Peridynamics for multiscale materials modeling [C]//Journal of Physics: Conference Series. [S. l.]: IOP Publishing, 2008, 125(1): 012078.
[10] XU J, ASKARI A, WECKNER O, et al. Peridynamic analysis of impact damage in composite laminates[J]. Journal of Aerospace Engineering, 2008, 21(3): 187–194. DOI:10.1061/(ASCE)0893-1321(2008)21:3(187)
[11] 沈峰, 章青, 黄丹, 等. 冲击荷载作用下混凝土结构破坏过程的近场动力学模拟[J]. 工程力学, 2012(supple.1): 12–15.
SHEN Feng, ZHANG Qing, HUANG Dan, et al. Peridynamics modeling of failure process of concrete structure subjected to impact loading[J]. Engineering Mechanics, 2012(supple.1): 12–15.
[12] 胡祎乐, 余音, 汪海. 基于近场动力学理论的层压板损伤分析方法[J]. 力学学报, 2013, 45(4): 624–628.
HU Yi-le, YU Yin, WANG Hai. Damage analysis method for laminates based on peridynamic theory[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(4): 624–628. DOI:10.6052/0459-1879-12-368
[13] 谷新保, 周小平. 含圆孔拉伸板的近场动力学数值模拟[J]. 固体力学学报, 2015, 36(5): 376–383.
GU Xin-bao, ZHOU Xiao-ping. The numerical simulation of tensile plate with circular hole using peridynamic theory[J]. Chinese Journal of Solid Mechanics, 2015, 36(5): 376–383.
[14] 郁杨天, 章青, 顾鑫. 含单边缺口混凝土梁破坏的近场动力学模拟[J]. 工程力学, 2016, 33(12): 80–85.
YV Yang-tian, ZHANG Qing, GU Xin. Impact failure simulation of a single-edged notched concrete beam based on peridynamics[J]. Engineering Mechanics, 2016, 33(12): 80–85.
[15] LAI X, REN B, FAN H, et al. Peridynamics simulations of geomaterial fragmentation by impulse loads[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39(12): 1304–1330. DOI:10.1002/nag.v39.12
[16] LAI X, LIU L S, LIU Q W, et al. Slope stabilityanalysis by peridynamic theory [C]//Applied Mechanics and Materials. [S. l.]: Trans Tech Publications, 2015, 744: 584-588.
[17] ZHOU X P, GU X B, WANG Y T. Numerical simulations of propagation, bifurcation and coalescence of cracks in rocks[J]. International Journal of Rock Mechanics and Mining Sciences, 2015, 80: 241–254. DOI:10.1016/j.ijrmms.2015.09.006
[18] REN B, FAN H, BERGEL G L, et al. A peridynamics-SPH coupling approach to simulate soil fragmentation induced by shock waves[J]. Computational Mechanics, 2015, 55(2): 287–302. DOI:10.1007/s00466-014-1101-6
[19] FAN H, BERGEL G L, LI S. A hybrid peridynamics-SPH simulation of soil fragmentation by blast loads of buried explosive[J]. International Journal of Impact Engineering, 2016, 87: 14–27. DOI:10.1016/j.ijimpeng.2015.08.006
[20] MACEK R W, SILLING S A. Peridynamics via finite element analysis[J]. Finite Elements in Analysis and Design, 2007, 43(15): 1169–1178. DOI:10.1016/j.finel.2007.08.012
[21] KILIC B, MADENCI E. Coupling of peridynamic theory and the finite element method[J]. Journal of Mechanics of Materials and Structures, 2010, 5(5): 707–733. DOI:10.2140/jomms
[22] LIU W, HONG J W. A coupling approach of discretized peridynamics with finite element method[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 245: 163–175.
[23] HUANG D, LU G, QIAO P. An improved peridynamic approach for quasi-static elastic deformation and brittle fracture analysis[J]. International Journal of Mechanical Sciences, 2015, 94: 111–122.
[24] AZIZI M A B. The peridynamic model of viscoelastic creep and recovery[J]. Multidiscipline Modeling in Materials and Structures, 2015, 11(4): 579–597. DOI:10.1108/MMMS-03-2015-0017
[25] MADENCI E, OTERKUS S. Ordinary state-based peridynamics for plastic deformation according to von Mises yield criteria with isotropic hardening[J]. Journal of the Mechanics and Physics of Solids, 2015, 86: 192–219.
[26] GHAJARI M, IANNUCCI L, CURTIS P. A peridynamic material model for the analysis of dynamic crack propagation in orthotropic media[J]. Computer Methods in Applied Mechanics and Engineering, 2014, 276(7): 431–452.
[27] SILLING S A, ASKARI E. A meshfree method based on the peridynamic model of solid mechanics[J]. Computers and Structures, 2005, 83(17): 1526–1535.
[28] HUANG D, LU G, WANG C, et al. An extended peridynamic approach for deformation and fractureanalysis[J]. Engineering Fracture Mechanics, 2015, 141: 196–211. DOI:10.1016/j.engfracmech.2015.04.036
[29] 铁摩辛柯. 弹性理论[M]. 北京: 高等教育出版社, 1990: 100-107.
[30] JENQ Y S, SHAH S P. Mixed-mode fracture of concrete[J]. International Journal of Fracture, 1988, 38(2): 123–142.
[31] JOHN R, SHAH S P. Mixed-mode fracture of concrete subjected to impact loading[J]. Journal of Structural Engineering, 1990, 116(3): 585–602. DOI:10.1061/(ASCE)0733-9445(1990)116:3(585)