引 言

损伤与破坏问题的准确建模是计算力学的重点研究课题之一^[1,4]. 近场动力学方法(peridynamics, PD)是基于非局部积分思想发展的数值方法,该方法基于粒子间的长程作用建模,通过近场范围体现结构特征尺度,避免了不连续处变量空间导数不存在而导致的奇异性问题,适用于分析含缺陷、裂纹等不连续问题的力学模型^[5]. 作为一种新兴理论,PD 自问世起便备受关注,近年来许多专家学者开始从事这方面的研究^[6,15].

虽然PD在模拟裂纹扩展时有很大的优势. 但作为非局部理论,PD计算效率远低于有限元法(finite element method, FEM). 为了发挥PD处理非连续问题以及FEM计算效率高的优势,许多学者在PD与FEM结合方面做了很多研究^[16,25]. Macek等^[12]将模型划分为PD区域与FEM区域,对两个区域重叠区域中的PD物质点与FEM结点采用位移协调的方法, PD物质点的位移可通过FEM结点位移进行插值后确定,但此过程不可逆,故边界条件只能在FEM区域实施. 文献^[21,22,23]发展了混合函数方法,该方法在局部模型和非局部模型的重合区域采用混合函数实现以上两种模型的平滑过渡, 用局部模型的应变能密度与非局部模型的应变能密度混合得到耦合区域的应变能密度,从而得到耦合区域的本构模型. Shojaei等^[24]将FEM和PD耦合求解动态裂纹扩展问题,首先将求解区域划分为3个区域,FEM区域、PD区域以及耦合区域. 其中危险区域采用PD粒子进行离散,其他区域采用有限元节点计算,该方案的优势在于在耦合区域内不需要设置混合函数,实施简单.

近场动力学在求解位移场以及破坏问题方面已有广泛的应用,但很少有文章提及近场动力学在热传导问题中的应用. Bobaru等^[26,27]根据能量平衡推导出键型PD热传导方程,算例表明解析解与PD热传导方程得到的结果相吻合. 此外,文章分析了近场范围 δ 对结果收敛性的影响. 根据欧拉-拉格朗日方程,Oterkus等^[28]推导出态型PD热传导方程,并采用虚拟层实现温度边界条件的施加. Chen等^[29]通过形状因子整合将热传导内核函数,讨论了形状因子对PD中 δ 收敛性的影响. 王飞等^[30]将PD热传导方程发展到非均质材料,研究了PD 方法中内核参数对非均匀材料热传导数值解的影响. 综上所述,目前还没有文献报道PD与FEM耦合求解热传导问题.

本文提出了一种求解稳态热传导问题的PD与FEM耦合方法,该方案不设置PD与FEM的重合区域,故不需要引入混合函数, 实施简单方便. 将PD与FEM写入同一个总体热传导矩阵,再采用Guyan缩聚法,进一步减小计算量.

1 近场动力学与有限元法耦合方案

键型PD热传导方程^[28]可以写为

$\rho c_v \dot {\varTheta }\left( {{\pmb x},t} \right) = \int_H f_h \Big( \varTheta \left( {\pmb x}',t \right) - \varTheta \left( {{\pmb x},t} \right),{\pmb x}' - {\pmb x},t \Big ) d V_{x}' + \rho S_{\rm b} \left( {{\pmb x},t} \right) (1) $

式中, $\rho $ 代表材料密度,$c_v $为材料比热容,$H$表示以$\delta $为半径的粒子${\pmb x}$的邻域,$f_{\rmh}$代表热流密度函数,$S_{\rm b}$为内热源. $t$和$\varTheta $分别表示时间和温度,$V_{x}'$代表粒子${x}'$所占的体积.

对于不计内热源稳态热传导方程可写为

$ \int_H {f_{\rm h} \left( {\varTheta \left( {{\pmb x}',t} \right) - \varTheta \left( {{\pmb x},t} \right),{\pmb x}' - {\pmb x},t} \right)} d V_{x}'=0 (2) $

将式(2)写为离散形式

$ \sum_{x}' {f_{\rm h} \left( {\varTheta \left( {{\pmb x}',t} \right) - \varTheta \left( {{\pmb x},t} \right),{\pmb x}' - {\pmb x},t} \right)} V_{x}'=0 (3) $

根据线性化的键型近场动力学理论^[16,31], f h 可写为

$ f_{\rm h} = \dfrac{\kappa }{\left| \xi \right|}\left( {\varTheta \left( {{\pmb x}'} \right) - \varTheta \left( {\pmb x} \right)} \right) (4)$

式中,$\left| \xi \right| = \left| {{\pmb x}' - {\pmb x}} \right|$,$ \kappa $代表微热传导系数. 为了与长程力更好地匹配, κ 不取常数,而是取一个与距离有关的量^[32].

$ \left.\!\!\begin{align} \kappa = \dfrac{2k}{A\delta ^2}\left[ {1 - \left( {\dfrac{\left| \xi \right|}{\delta }} \right)^2} \right]^2 , \ \hbox{one-dimensional} \\ \kappa = \dfrac{6k}{\pi h\delta ^3}\left[ {1 - \left( {\dfrac{\left| \xi \right|}{\delta }} \right)^2} \right]^2 , \ \hbox{two-dimensional} \\ \kappa = \dfrac{6k}{\pi \delta ^4}\left[ {1 - \left( {\dfrac{\left| \xi \right|}{\delta }} \right)^2} \right]^2 , \ \hbox{three-dimensional} \end{align}\!\!\right\} (5) $

式中,A代表横截面积,h表示平面厚度,k为热传导系数, δ 代表以近场作用半径.

对于处于一种材料中的两个PD粒子,微热传导系数 κ 为常数. 处于两种材料界面附近的PD粒子会不可避免地与另外一种材料中的PD粒子相互作用,如图1所示.

对于处于两种材料中的两个PD粒子的微热传导系数 κ 用下式表达^[8]

$ \kappa = \dfrac{l_1 + l_2 }{{l_1 } / {\kappa _1 } + {l_2 } / {\kappa _2 }} (6) $

式中,$l_1 $代表粒子${\pmb x}'$与${\pmb x}$的距离$\left| \xi \right|$在材料1中的部分(图中实线部分),其对应的微热传导系数为$\kappa _{1} $;$l_{2} $代表粒子${\pmb x}'$与${\pmb x}$的距离$\left| \xi \right|$在材料2中的部分(图中虚线部分),其对应的微热传导系数为$\kappa_{2} $. 从图2不难看出$\left| \xi \right| = l_1 + l_2 $,$\kappa_{1} $和 $\kappa _{2} $分别由对应材料的热导率$k_{1}$, $k_{2} $代入式(5)得到.

对于一维热传导问题,FEM单元热传导矩阵可以写为

$ {\pmb K}_{\rm FEM} = \dfrac{kA}{l}\left[ \!\! \begin{array}{cc} 1 & { - 1} \\ { - 1} & 1 \end{array} \!\! \right] = a\left[\!\! \begin{array}{cc} 1 & { - 1} \\ { - 1} & 1 \end{array} \!\! \right] (7) $

式中,$a = \dfrac{kA}{l}$.

根据式(4)和式(5),键型PD中粒子 x ' 与 x 的作用可写为矩阵形式

$ {\pmb K}_{\rm PD} = \dfrac{\kappa }{l}V_x V_{x'} \left[ \!\!\begin{array}{cc} 1 & { - 1} \\ { - 1} & 1 \end{array} \!\! \right] = b\left[\!\!\begin{array}{cc} 1 & { - 1} \\ { - 1} & 1 \end{array} \!\! \right] (8) $

式中,$b = \dfrac{\kappa }{l}V_x V_{x'} $,$ V_{x}$代表粒子${\pmb x}$所占的空间,$V_{x'} $代表粒子${\pmbx}'$所占的空间,对于一维情况,$V_x = A\Delta x$,$A$代表一维模型的横截面积.

对于一维模型,FEM与PD的耦合方法如图2所示. 图中圆点代表FEM节点,方块代表PD粒子. 直线代表有限元节点的作 用,曲线代表PD粒子的作用. 其中,FEM节点与其左右两边的节点(包括FEM节点和PD粒子)以有限元方式相互作用, PD粒子与其域内的节点(包括FEM节点和PD粒子) 以PD方式相互作用. 对于 δ = 3 Δ x 的情况,无论近场 动力学粒子在何位置,均与其周围的6个粒子相互作用. 对于图2(a)所示的PD粒子 b 1 ,与其左边的3个FEM节 点 a 2 , a 3 , a 4 以及右边的3个PD粒子 b 2 , b 3 , b 4 以PD方式作用;同理,对于图2(b)所示 的PD粒子 b 2 ,与其周围的2个FEM节点 a 3 , a 4 以及4个PD粒子 b 1 , b 3 , b 4 , b 5 以PD方 式作用.图2中FEM节点 a 4 与其左边的FEM节点 a 4 以及右边的PD粒子 b 1 以FEM方式作用. 对于图2所示的一维 模型,总体热传导矩阵可写为

$\left[\!\!\begin{array}{ccccccccccccc} a & -a & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots \\ -a & {2a} & -a & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ 0 & -a & {2a} & -a & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ 0 & 0 & -a & {2a} & -a & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\ 0 & -b/3 & -b/2 & -b & {8b /3} & -b & -b/2 & -b/3 & 0 & 0 & 0 & 0 & \\ 0 & 0 & -b/3 & -b/2 & -b & {8b / 3} & -b & -b/2 & -b/3 & 0 & 0 & 0 & \\ 0 & 0 & 0 & -b/3 & - b/ 2 & -b & 8b /3 & -b & -b/2 & -b/3 & 0 & 0 & \cdots \\ 0 & 0 & 0 & 0 & -b/3 & -b/2 & -b & {8b / 3} & -b & -b/2 & - b/ 3 & 0 & \\ 0 & 0 & 0 & 0 & 0 & - b /3 & -b/2 & -b & {8b/ 3} & -b & -b/2 & -b/3 & \\ \cdots & & & & & & \cdots & & & \cdots & & & \cdots \\ \cdots & & & & & & \cdots & & & & \cdots & & \cdots \\ \cdots & & & & & & \cdots & & & & & \cdots & \cdots \\ \cdots & & & & & & \cdots & & & & & & \cdots \end{array}\!\! \right]\cdot \left[\!\! \begin{array}{c} {\varPhi _{1} } \\ {\varPhi _{2} } \\ {\varPhi _{3} } \\ {\varPhi _{4} } \\ {\varPhi _{5} } \\ {\varPhi _{6} } \\ {\varPhi _{7} } \\ {\varPhi _{8} } \\ {\varPhi _{9} } \\ \cdots \\ \cdots \\ \cdots \\ \cdots \\ \end{array} \!\! \right] = \left[\!\! \begin{array}{c} {P_{1} } \\ {P_{2} } \\ {P_{3} } \\ {P_{4} } \\ {P_{5} } \\ {P_{6} } \\ {P_{7} } \\ {P_{8} } \\ {P_{9} } \\ \cdots \\ \cdots \\ \cdots \\ \cdots \end{array}\!\! \right] (9)$

从式(9)可以看出,耦合后的总体热传导矩阵不再对称,式中第1至第4行为FEM区域对总体热传导矩阵的贡献,第5至第9行为PD区域对总体热传导矩阵的贡献. 通过式(9)可以清晰的看出,近场动力学区域的宽度比有限元区域的宽. 式(9)第4行对应图2中FEM节点 a 4 ,该节点与两边的粒子采用FEM方式相互作用,故该行只有3个元素不为0. 式(9)第5行对应图2中PD粒子 b 1 ,其与周围的粒子采用PD方式相互作用,故该行有7个元素不为0. 此外,右端列向量代表节点的热流.

与一维情况类似,二维情况下FEM与PD的耦合方法如图3所示.

两区域交界附近的PD粒子 c 与其域内的节点(包括FEM节点和PD粒子) 以PD方式相互作用;FEM节点 d (本文中,二维FEM均采用四节点线性热传导单元^[33]) 与其周围的节点(包括FEM 节点和PD 粒子) 以有限元方式相互作用. 将PD 与FEM 写入同一个总体热传导矩阵如下

$ \left[\!\! \begin{array}{ccccc} & \vdots & & \vdots & \\ \cdots & \sum_j^N {k_{cj}^{11} } & \cdots & {k_{cd}^{12} } & \cdots \\ & \vdots & & \vdots & \\ \cdots & {k_{dc} } & \cdots & {\sum_i^M {k_{dd} } } & \cdots \\ & \vdots & & \vdots & \end{array} \!\! \right] \cdot \left\{ \!\!\begin{array}{c} \vdots \\ {\varTheta _c } \\ \vdots \\ {\varTheta _d } \\ \vdots \end{array} \!\! \right\} = \left\{ \!\! \begin{array}{c} \vdots \\ {P_c } \\ \vdots \\ {P_d } \\ \vdots \end{array} \!\! \right\} (10)$

式中,$P_c $和$P_d $代表节点的热流.

$ P_c = -\sum_i^M\dfrac {1} {2}q S (11) $

其中, q 代表边界热流密度, S 代表单元边界长度, M 代表FEM节点 d 所在单元的总和.

P d = - q V d (12)

其中 嚃 为粒子 d 所占有的体积. N 代表PD粒子 c 作用域内的粒子总数, k cj 11 和 k cd 12 代表PD方式的作用

$ {\pmb K}_{\rm PD}^{cd} = \dfrac{\kappa }{l}V_c V_d \left[ \!\! \begin{array}{cc} 1 & { - 1} \\ { - 1} & 1 \end{array} \!\! \right]{ = }\left[\!\! \begin{array}{cc} {k_{cd}^{11} } & {k_{cd}^{12} } \\ {k_{cd}^{21} } & {k_{cd}^{22} } \end{array} \!\! \right] (13) $

其中, k dc 和 k dd 为FEM单元热传导矩阵中的元素^[33]

式(10) 中的整体热传导矩阵由FEM 的单元热传导矩阵与PD 粒子热传导矩阵装配得到,其中PD 粒子热传导矩阵的装配方法与FEM 的相同.

将式(10) 写为分块格式

$ \left[\!\! \begin{array}{cc} {k_{ii} } & {k_{io} } \\ {k_{oi} } & {k_{oo} } \end{array} \!\! \right]\left[\!\! \begin{array}{c} {\varTheta _i } \\ {\varTheta _o } \end{array} \!\! \right] = \left[\!\! \begin{array}{c} {P_i } \\ {P_o } \end{array} \!\! \right] (14)$

式中,下标i代表从自由度,下标o代表主自由度,采用Guyan缩聚法^[34],式(14)可写为

$ k_{c0} \varTheta _o = P_{c0} (15) $

式(15)为缩聚后的有限元方程,$k_{c0} $的阶数与$k_{oo} $的相同,其中

$k_{c0} { = }k_{oo} - k_{oi} k_{ii}^{ - 1} k_{io} (16)$

$P_{c0} { = }P_o - k_{ii}^{ - 1} k_{io} P_i (17)$

通过式(15) 可以得到主自由度的温度值. 从自由度的温度可以由式(14) 第1 行得到

$ \Theta _i { = } - k_{ii}^{ - 1} \left( {P_i - k_{io} \varTheta _o } \right) (18) $

采用Guyan法将PD和FEM耦合求解热传导问题的计算效率进一步提高.

2 数值算例

2.1 一维情况下热传导算例

如图4所示的一维热传导模型,模型长 l z = 100 m ,其中左半部分( x < l z / 2 ) 由材料1组成,其热传导系数为 k 1 = 2 W / m ⋅ K ;右半部分 ( x > l z 2 ) 由材料2组成,其热传导系数为 k 2 = 1 W / m ⋅ K .PD区域的微热传导系数由式(5)和式(6)确定. 将模型划分为101个节点,其中节点31至节点70为PD粒 子(用正方形表示),其他节点为FEM节点(用圆形表示). PD粒子邻域半径 δ = 3 Δ x .对于Guyan法,通常把对接边界节点的自由度划分为主自由度,其他节点的自由度划分为从自由度. 本算例中,取边界节点以及边界附近的节点自由度为主自由度,其他节点自由度为从自由度. 取节点1 ~ 3,99 ~ 101为主自由度,其他节点为从自由度,如图4所示.

分别取两类边界条件,对比不同边界下解析解与本文耦合方法的计算结果.

第一类边界条件

Θ x = 0 = 0 K , Θ x = l z = 100 K (19)

第二类边界条件

Θ x = 0 = 0 K , q x = l z = 1 W / m 2 (20)

其中, q 代表热流密度.

本文耦合方法与解析解的计算结果如图5所示. 从图中可以看出,采用本文耦合方法得到的结果与解析解吻合很好.

表1为一维温度场模型正则化的绝对误差以及相对误差. 绝对误差由下式计算

$ \left| {e_{\rm abs} } \right| = \sqrt {\sum_{i = 1}^{N_{1} } {e_{{\rm abs},i}^2 } } (21) $

式中, N 1 代表节点总数, e abs , i 为节点 i 的绝对误差, e abs , i = Θ * - Θ , Θ * 为本文耦合方法计算得到的温度值, Θ 为解析解结果.

相对误差的表达式为

$ \left| {e_{\rm rel} } \right| = \sqrt {\sum_{i = 1}^{N_{1} } {e_{{\rm rel},i}^2 } } (22) $

式中, e rel , i 为节点 i 的绝对误差, e rel , i = Θ * - Θ / Θ .只对没有施加边界条件的节点计算误差,从表1可以看出,两种边界条件的绝对误差与相对误差都接近于0.

2.2 二维情况下热传导算例

针对二维热传导模型如图6所示,模型取长 a = 100 m,宽 b = 100 m的方板,方板厚度 h = 1 m. 模型上半部分( k 1 = 2 W / m ⋅ K )由材料1组成,其热传导系数 y < b / 2 ) ;下半部分( k 2 = 1 W / m ⋅ K 由材料2组成,其热传导系数为 0.3 a # x 2264 ; x # x 2264 ; 0.7 a .PD区域的微热传导系数由式(5)和式(6)确定. 如图6(a),将模型分为两部分, 其中区域 0.4 b # x 2264 ; y # x 2264 ; 0.6 b , k oi k ii - 1 k io 采用PD粒子描述,其他区域采用FEM节点描述. 本例中,除了取边界节点自由度为主自由度以外,PD区域粒子的自由度也取为主自由度,如图6(b)所示. 这样做的目的在于,在做破坏分析时,破坏只在发生在PD区域,其他节点的热传导矩阵(刚度矩阵)以及节点的热流(载荷矩阵) 在计算过程中不会发生变化,式(16)中的 k ii - 1 k io P i 和式(17)中的 k oo 只需在第一步需要计算,其他计算步骤只需要更新式(16)中的 P o 和式(17)中的 δ .

接下来讨论本文耦合算法的收敛性. 对 m = 4 收敛性分析,取 Δ x = 2 ,分别取 δ = 8 m ( Δ x = 1 m)、 δ = 4 m ( Δ x = 0.5 m)和 δ = 2 m ( δ m). 边界条件按式(23)选取

$ \left. \varTheta \right|_{x = 0} = \left. \varTheta \right|_{x = a} = \left. \varTheta \right|_{y = 0} = 0 {\rm K} , \ \ \left. \varTheta \right|_{y = b} = 100 {\rm K} (23) $

采用本文耦合方法得到的结果与ABAQUS计算得到的FEM结果对比.图7为不同 dT / dy ) 值时温度梯度( δ 的误差沿横向的变化曲线. 通过图7可以看出,本文的耦合方案在 δ 取不同值的情况下均有很高的计算精度(最大误差不超过1%). 网格密度越小,本文耦合方案的计算精度越高. 误差最大值处于两种材料界面附近的位置,这说明文献^[4]中的算法(式(6))在两种材料界面附近的计算精度比较低.

对于 δ = 4 收敛性分析,取 Δ x = 2 m,分别取 Δ x = 1 m, Δ x = 0.5 m和 Δ x = 2 m,对于 m = 2 m,每个PD粒子与周围12个粒子以PD方式相互作用 ( Δ x = 1 );对于 m = 4 m,每个PD粒子与周围48个粒子以PD方式相互作用 ( Δ x = 0.5 );对于 m = 8 m,每个PD粒子与周围196个粒子以PD方式相互作用( δ ).边界条件仍按式(23)选取.图8为不同 Δ x = 1 值时温度梯度的误差沿横向的变化曲线. 通过图8可以看出,对 于 Δ x = 0.5 m的情况,本文耦合方案的误差最小. Δ x = 1 m时,误差反而大于 Δ x = 1 m的情况,这是因为与 Δ x = 0.5 m的情况相比, Δ x = 2 m时两种材 料界面附近的PD粒子要与更多另外一种材料中的粒子相互作用,从而增大了误差. m m时,虽然材料界面附近的PD粒子与另外一种材料中较少的粒子相互作用,但其对应的 m 值仅为2,致使其误差最大.

通过以上 m 收敛性分析以及 m 收敛性分析,可以得到结论,为了提高材料界面处的计算精度, m 值不宜取太大,另一方面,为了保 证PD的计算精度, m 值不能太小. 通过上述算例分析, m 取4是一个恰当的选择. 在 Δ x = 1 值恰当的前提下,提高网格密度可以提高计算精度. 在以下二维算例中,取 m = 4 m, δ = 4 , {Invalid MML} m. 边界条件分别按式(23)和式(24)选取,其中式(24)为第二类边界条件.

$ \left. \varTheta \right|_{x = 0} = \left. \varTheta \right|_{x = a} = \left. \varTheta \right|_{y = 0} = 0 {\rm K} , \ \ \left. q \right|_{y = b} = 1 {\rm W} /{\rm m}^{2} (24) $

分别采用经典有限元理论和本文耦合方法对本节的模型进行求解(均采用显式计算). 模拟结果如图9所示,其中左边一组为 有限元模拟结果(ABAQUS),右边一组为本文耦合方法的结果. 由图9可知,本文耦合方法与有限元法的模拟结果吻合很好. 图形的上半部分温度梯度较大,下半部分温度 梯度很小. 对于第一类边界条件,采用两种方法计算的温度最大值均为100 K;对于第二类边界条件,采用有限元法计算得到温度最大值为190.83 K,采用本文方法得到温度最大值为190.98 K,相对误差为0.08%.

2.3 含裂纹情况下二维热传导算例

在2.2节二维热传导模型的基础上,在模型中间开一个裂纹(如图10所示),裂纹位置由下式表达

$ 0.25a < x < 0.75a , \ \ y = 0.5b (25) $

边界条件仍按式(19)、式(20)选取. PD区域引入初始裂纹有两种方法,一种为破坏经过裂纹线上的所有键, 如图11(a)所示. 另一种方法为移除裂纹线上的粒子同时破坏经过裂纹线上的所有键^[35],如图11(b)所示,本文采用第二种形式.

分别采用传统FEM理论、本文耦合方案和PD理论对含裂纹二维温度场进行模拟,采用式(23) (第一类边 界条件)的模拟结果如图12所示,采用式(24) (第二类边界条件)的模拟结果如图13所示. 从图12和图13可以看出,在以上两类边界条件下,裂纹上下表面的温度梯度最大,这是因为裂纹的存在阻断了热量的传播. 由于裂纹的阻隔以及左右边界的散热,热量几乎不能传到裂纹下侧.

对比图9与图12,图13可知,模型中的裂纹对温度场分布造成了很大的影响. 此外,本文耦合方案的结果与传统FEM理论、PD理论的结果相吻合.

需要指出,本文的耦合方案以及PD模拟结果由编制的MATLAB程序实现,程序采用显式计算,对以上两个算例均采用1000个 增量步进行计算,计算时间如表2所示

从表2可以看出,本文耦合方案和PD模型相比可以在保证计算精度的前提下将计算效率提高5倍左右,这充分说明了本文耦合方案的在提高计算效率上的优势.

本文计算采用计算机的硬件条件如下：

处理器: i5-4590 CPU @3.30 GHz

内存: 16.00 GB

操作系统: Windows 10 Enterprise 64 bit

3 结 论

本文给出了有效处理含裂纹结构热传导问题的一种新的有限元与近场动力学耦合方法. 该方法结合了近场动力学处理含裂纹问题和有限单元法处理连续问题的优势,并应用Guyan缩聚法进一步减小计算量. 通过数值算例表明,本文耦合方法得到的温度场结果准确. 利用该耦合方案可以进一步拓展到热力耦合条件下含裂纹材料和结构的裂纹扩展问题.

The authors have declared that no competing interests exist.

---

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子

[1]	马天宝 , 任会兰 , 李健 等 . 爆炸与冲击问题的大规模高精度计算 . 力学学报, 2016 , 48 ( 3 ): 599 - 608 [本文引用: 1] ( Ma Tianbao , Ren Huilan , Li Jian , et al . Large scale high precision computation for explosion and impact problems . Chinese Journal of Theoretical and Applied Mechanics, 2016 , 48 ( 3 ): 599 - 608 (in Chinese)) [本文引用: 1]
[2]	吕海宝 , 梁军 , 郭旭 等 . 第七届全国固体力学青年学者学术研讨会报告综述 . 力学学报, 2017 , 49 ( 1 ): 223 - 230 ( Lü Haibao , Liang Jun , Guo Xu , et al . Review of the seventh national symposium on solid mechanics for young scholars . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 1 ): 223 - 230 (in Chinese))
[3]	孙可明 , 张树翠 . 含层理页岩气藏水力压裂裂纹扩展规律解析分析 . 力学学报, 2016 , 48 ( 5 ): 1229 - 1237 ( Sun Keming , Zhang Shucui . Hydraulic fracture propagation in shale gas bedding reservoir analytical analysis . Chinese Journal of Theoretical and Applied Mechanics, 2016 , 48 ( 5 ): 1229 - 1237 (in Chinese))
[4]	路德春 , 李萌 , 王国盛 等 . 静动组合载荷下混凝土率效应机理及强度准则 . 力学学报, 2017 , 49 ( 4 ): 940 - 952 [本文引用: 2] ( Lu Dechun , Li Meng , Wang Guosheng , et al . Study on strain rate effect and strength criterion of concrete under static-dynamic coupled loading . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 4 ): 940 - 952 (in Chinese)) [本文引用: 2]
[5]	Silling SA , Lehoucq RB . Peridynamic theory of solid mechanics . Advances in Applied Mechanics, 2010 , 44 : 73 - 168 [本文引用: 1]
[6]	Silling SA . Origin and effect of nonlocality in a composite . Journal of Mechanics of Materials and Structures, 2014 , 9 ( 2 ): 245 - 258 [本文引用: 1]
[7]	Gu X , Zhang Q , Huang D , et al . Wave dispersion analysis and simulation method for concrete SHPB test in peridynamics . Engineering Fracture Mechanics, 2016 , 160 : 124 - 137
[8]	Madenci E , Oterkus E . Peridynamic Theory and Its Application . New York: Springer, 2014 : 203 - 244 [本文引用: 1]
[9]	章青 , 顾鑫 , 郁杨天 . 冲击载荷作用下颗粒材料动态力学响应的近场动力学模拟 . 力学学报, 2016 , 48 ( 1 ): 56 - 63 ( Zhang Qing , Gu Xin , Yu Yangtian . Peridynamics simulation for dynamic response of granular materials under impact loading, Chinese Journal of Theoretical and Applied Mechanics, 2016 , 48 ( 1 ): 56 - 63 (in Chinese))
[10]	黄丹 , 章青 , 乔丕忠 等 . 近场动力学方法及其应用 . 力学进展, 2010 , 40 ( 4 ): 448 - 459 ( Huang Dan , Zhang Qing , Qiao Pizhong , et al . A review on peridynamics method and its application . Advance in Mechanics, 2010 , 40 ( 4 ): 448 - 459 (in Chinese))
[11]	胡祎乐 , 余音 , 汪海 . 基于近场动力学理论的层压板损伤分析方法 . 力学学报, 2013 , 45 ( 4 ): 624 - 628 ( Hu Yile , Yu Yin , Wang Hai . Damage analysis method for laminates based on peridynamic theory . Chinese Journal of Theoretical and Applied Mechanics, 2013 , 45 ( 4 ): 624 - 628 (in Chinese))
[12]	Hu YL , Madenci E . Peridynamics for fatigue life and residual strength prediction of composite laminate . Composite Structures, 2017 , 160 : 169 - 184 [本文引用: 1]
[13]	Hu YL , De Carvalho NV , Madenci E . Peridynamic modeling of delamination growth in composite laminates . Composite Structures, 2015 , 132 : 610 - 620
[14]	Chen Z , Bakenhus D , Bobaru F . A constructive peridynamic kernel for elasticity . Computer Methods in Applied Mechanics and Engineering, 2016 , 311 : 356 - 373
[15]	Yaghoobi A , Chorzepa MG . Meshless modeling framework for fiber reinforced concrete structures . Computers & Structures, 2015 , 161 : 43 - 54 [本文引用: 1]
[16]	Macek RW , Silling SA . Peridynamics via finite element analysis . Finite Elements in Analysis and Design, 2007 , 43 ( 15 ): 1169 - 1178 [本文引用: 2]
[17]	Kilic B , Madenci E . Coupling of peridynamic theory and the finite element method . Journal of Mechanics of Materials and Structures, 2010 , 5 ( 5 ): 707 - 733
[18]	Liu W , Hong JW . A coupling approach of discretized peridynamics with finite element method . Computer Methods in Applied Mechanics and Engineering, 2012 , 245 : 163 - 175
[19]	Galvanetto U , Mudric T , Shojaei A , et al . An effective way to couple FEM meshes and peridynamics grids for the solution of static equilibrium problems . Mechanics Research Communications, 2016 , 76 : 41 - 47
[20]	Seleson P , Beneddine S , Prudhomme S . A force-based coupling scheme for peridynamics and classical elasticity . Computational Materials Science, 2013 , 66 : 34 - 49
[21]	Lubineau G , Azdoud Y , Han F , et al . A morphing strategy to couple non-local to local continuum mechanics . Journal of the Mechanics and Physics of Solids, 2012 , 60 ( 6 ): 1088 - 1102 [本文引用: 1]
[22]	Azdoud Y , Han F , Lubineau G . A morphing framework to couple non-local and local anisotropic continua . International Journal of Solids and Structures, 2013 , 50 ( 9 ): 1332 - 1341 [本文引用: 1]
[23]	Han F , Lubineau G , Azdoud Y , et al . A morphing approach to couple state-based peridynamics with classical continuum mechanics . Computer Methods in Applied Mechanics and Engineering, 2016 , 301 : 336 - 358 [本文引用: 1]
[24]	Shojaei A , Mudric T , Zaccariotto M , et al . A coupled meshless finite point/Peridynamic method for 2D dynamic fracture analysis . International Journal of Mechanical Sciences, 2016 , 119 : 419 - 431 [本文引用: 1]
[25]	Bie YH , Cui XY , Li ZC . A coupling approach of state-based peridynamics with node-based smoothed finite element method . Computer Methods in Applied Mechanics and Engineering, 2018 , 331 : 675 - 700 [本文引用: 1]
[26]	Bobaru F , Duangpanya M . The peridynamic formulation for transient heat conduction . International Journal of Heat and Mass Transfer, 2010 , 53 ( 19 ): 4047 - 4059 [本文引用: 1]
[27]	Bobaru F , Duangpanya M . A peridynamic formulation for transient heat conduction in bodies with evolving discontinuities . Journal of Computational Physics, 2012 , 231 ( 7 ): 2764 - 2785 [本文引用: 1]
[28]	Oterkus S , Madenci E , Agwai A . Peridynamic thermal diffusion . Journal of Computational Physics, 2014 , 265 : 71 - 96 [本文引用: 2]
[29]	Chen Z , Bobaru F . Selecting the kernel in a peridynamic formulation: A study for transient heat diffusion . Computer Physics Communications, 2015 , 197 : 51 - 60 [本文引用: 1]
[30]	王飞 , 马玉娥 , 郭妍宁 . 近场动力学中内核参数对非均匀材料热传导数值解的影响研究 . 西北工业大学学报, 2017 , 35 ( 2 ): 203 - 207 [本文引用: 1] ( Wang Fei , Ma Yu’e , Guo Yanning . Effects of kernel parameters of peridynamic theory on heat conduction numerical solution for non-homogeneous material . Journal of Northwestern Polytechnical University, 2017 , 35 ( 2 ): 203 - 207 (in Chinese)) [本文引用: 1]
[31]	Silling SA . Reformulation of elasticity theory for discontinuities and long-range forces . Journal of the Mechanics and Physics of Solids, 2000 , 48 ( 1 ): 175 - 209 [本文引用: 1]
[32]	Cheng Z , Zhang G , Wang Y , et al . A peridynamic model for dynamic fracture in functionally graded materials . Composite Structures, 2015 , 133 : 529 - 546 [本文引用: 1]
[33]	孔祥谦 . 热应力有限单元法分析. 上海 : 上海交通大学出版社 , 1999 : 136 - 139 [本文引用: 2] ( Kong Xiangqian. Analysis of thermal stress finite element method. Shanghai : Shanghai Jiao Tong University Press , 1999 : 136 - 139 (in Chinese)) [本文引用: 2]
[34]	邱吉宝 , 向树红 , 张正平 . 计算结构动力学. 合肥 : 中国科学技术大学出版社 , 2009 : 352 - 354 [本文引用: 1] ( Qiu Jibao , Xiang Shuhong , Zhang Zhengping. Computational Structure Dynamics. Hefei : University of Science and Technology of China Press , 2009 : 352 - 354 (in Chinese)) [本文引用: 1]
[35]	Ha YD , Bobaru F . Studies of dynamic crack propagation and crack branching with peridynamics . International Journal of Fracture, 2010 , 162 ( 1 ): 229 - 244 [本文引用: 1]