刘硕¹, 方国东^1,*,, 王兵¹, 付茂青¹, 梁军^2,*,
1 哈尔滨工业大学特种环境复合材料技术国防科技重点实验室,哈尔滨 150001
2 北京理工大学先进结构技术研究院,北京 100081

STUDY OF THERMAL CONDUCTION PROBLEM USING COUPLED PERIDYNAMICS AND FINITE ELEMENT METHOD

LiuShuo¹, FangGuodong^1,*,, WangBing¹, FuMaoqing¹, LiangJun^2,*,
1 Science and Technology on Advanced Composites in Special Environments Key Laboratory, Harbin Institute of Technology, Harbin 150001, China
2 Institute of Advanced Structure Technology, Beijing Institute of Technology,Beijing 100081, China
中图分类号:O34

通讯作者:^*通讯作者：方国东，讲师，主要研究方向：复合材料建模分析及性能评价. E-mail: fanggd@hit.edu.cn;梁军，教授，主要研究方向：先进复合材料力学. E-mail: liangjun@bit.edu.cn^*通讯作者：方国东，讲师，主要研究方向：复合材料建模分析及性能评价. E-mail: fanggd@hit.edu.cn;梁军，教授，主要研究方向：先进复合材料力学. E-mail: liangjun@bit.edu.cn^*通讯作者：方国东，讲师，主要研究方向：复合材料建模分析及性能评价. E-mail: fanggd@hit.edu.cn;梁军，教授，主要研究方向：先进复合材料力学. E-mail: liangjun@bit.edu.cn
收稿日期:2017-10-9
接受日期:2018-02-24
网络出版日期:2018-02-10
版权声明:2018《力学学报》编辑部《力学学报》编辑部 所有
基金资助:国家自然科学基金资助项目(11672089,11732002,11325210,11421091).

展开

摘要
求解含裂纹等不连续问题一直是计算力学的重点研究课题之一,以偏微分方程为基础的连续介质力学方法处理不连续问题时面临很大的困难. 近场动力学方法是一种基于积分方程的非局部理论,在处理不连续问题时有很大的优越性. 本文提出了求解含裂纹热传导问题的一种新的近场动力学与有限元法的耦合方法. 结合近场动力学方法处理不连续问题的优势以及有限元方法计算效率高的优势,将求解区域划分为两个区域,近场动力学区域和有限元区域. 包含裂纹的区域采用近场动力学方法建模,其他区域采用有限元方法建模. 本文提出的耦合方案实施简单方便,近场动力学区域与有限元区域之间不需要设置重叠区域. 耦合方法通过近场动力学粒子与其域内所有粒子(包括近场动力学粒子和有限元节点)以非局部方式连接,有限元节点与其周围的所有粒子以有限元方式相互作用. 将有限元热传导矩阵和近场动力学粒子相互作用矩阵写入同一整体热传导矩阵中,并采用Guyan缩聚法进一步减小计算量. 分别采用连续介质力学方法和近场动力学方法对一维以及二维温度场算例进行模拟,结果表明,本文的耦合方法具有较高的计算精度和计算效率. 该耦合方案可以进一步拓展到热力耦合条件下含裂纹材料和结构的裂纹扩展问题.

关键词：近场动力学;耦合;有限元法;热传导
Abstract
To accurately model discontinuous problems with cracks is one important topic in computational mechanics. It is very difficult to solve discontinuous problems using continuum mechanics methods based on partial differential equations. However, peridynamics (PD), a non-local theory based on integral equations, has great advantages in solving these problems. In this paper, a new method is proposed to solve heat conduction problems with cracks using coupled PD and finite element method (FEM). This method has both the advantage of the computational efficiency of FEM and the advantage of PD in solving discontinuous problems. The computational domain can be partitioned into two regions, PD region and FEM region. The region containing the crack is modeled by PD, and the other region is modeled by FEM. Application of the coupling scheme proposed in this paper is simple and convenient, since there is no need to introduce an overlapping region between PD region and FEM region. As for the coupling approach, the PD particle is connected non-locally to all particles (PD particles and finite element nodes) within its horizon, whereas the finite element node interacts with other nodes in the finite element manner. The heat conduction matrices of FEM and the matrices of the interaction between PD particles are combined to be a global heat conduction matrix. The Guyan reduction method is used to further reduce the computational cost. The temperature fields of a one-dimensional bar and a two-dimensional plate obtained by the coupling approach are compared with classical solutions. Results show that the proposed coupling method is accurate and efficient. The coupling scheme can be extended to solve crack propagation problems with the thermo-mechanical load.

Keywords：peridynamics;coupling;finite element method;heat conduction

-->0
PDF (15242KB)元数据多维度评价相关文章收藏文章
本文引用格式导出EndNoteRisBibtex收藏本文-->
刘硕, 方国东, 王兵, 付茂青, 梁军. 近场动力学与有限元方法耦合求解热传导问题[J]. 力学学报, 2018, 50(2): 339-348 https://doi.org/10.6052/0459-1879-17-332
Liu Shuo, Fang Guodong, Wang Bing, Fu Maoqing, Liang Jun. STUDY OF THERMAL CONDUCTION PROBLEM USING COUPLED PERIDYNAMICS AND FINITE ELEMENT METHOD[J]. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(2): 339-348 https://doi.org/10.6052/0459-1879-17-332

引 言

损伤与破坏问题的准确建模是计算力学的重点研究课题之一^[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], fh可写为
$ 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所示.
显示原图|下载原图ZIP|生成PPT
图1界面附近的PD粒子
-->Fig. 1PD particles close to interface
-->

对于处于两种材料中的两个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)得到.
显示原图|下载原图ZIP|生成PPT
图2一维情况下耦合区域节点图
-->Fig. 2Nodes in coupling region for one-dimensional case
-->

对于一维热传导问题,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粒子 b1,与其左边的3个FEM节 点 a2, a3, a4以及右边的3个PD粒子 b2, b3, b4以PD方式作用;同理,对于图2(b)所示 的PD粒子 b2,与其周围的2个FEM节点 a3, a4以及4个PD粒子 b1, b3, b4, b5以PD方 式作用. 图2中FEM节点 a4与其左边的FEM节点 a4以及右边的PD粒子 b1以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节点 a4,该节点与两边的粒子采用FEM方式相互作用,故该行只有3个元素不为0. 式(9)第5行对应图2中PD粒子 b1,其与周围的粒子采用PD方式相互作用,故该行有7个元素不为0. 此外,右端列向量代表节点的热流.
与一维情况类似,二维情况下FEM与PD的耦合方法如图3所示.
显示原图|下载原图ZIP|生成PPT
图3二维情况下耦合区域节点图
-->Fig. 3Nodes in coupling region for two-dimensional case
-->

两区域交界附近的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所在单元的总和.
Pd=-qVd(12)
其中 嚃为粒子 d所占有的体积. N代表PD粒子 c作用域内的粒子总数, kcj11和 kcd12代表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) $
其中, kdc和 kdd为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所示的一维热传导模型,模型长 lz=100m,其中左半部分( x<lz/2)由材料1组成,其热传导系数为 k1=2W/m?K;右半部分 (x>lz2)由材料2组成,其热传导系数为 k2=1W/m?K. PD区域的微热传导系数由式(5)和式(6)确定. 将模型划分为101个节点,其中节点31至节点70为PD粒 子(用正方形表示),其他节点为FEM节点(用圆形表示). PD粒子邻域半径 δ=3Δx. 对于Guyan法,通常把对接边界节点的自由度划分为主自由度,其他节点的自由度划分为从自由度. 本算例中,取边界节点以及边界附近的节点自由度为主自由度,其他节点自由度为从自由度. 取节点1 ~3,99 ~101为主自由度,其他节点为从自由度,如图4所示.
显示原图|下载原图ZIP|生成PPT
图4一维热传导模型
-->Fig. 4One-dimensional heat conduction model
-->

分别取两类边界条件,对比不同边界下解析解与本文耦合方法的计算结果.
第一类边界条件
Θx=0=0K,Θx=lz=100K(19)
第二类边界条件
Θx=0=0K,qx=lz=1W/m2(20)
其中, q代表热流密度.
本文耦合方法与解析解的计算结果如图5所示. 从图中可以看出,采用本文耦合方法得到的结果与解析解吻合很好.
显示原图|下载原图ZIP|生成PPT
图 5一维温度场模型模拟结果
-->Fig. 5Simulation results for one-dimensional temperature field model
-->

Table 1
表 1
表 1一维温度场模型正则化的误差
Table 1Regularization error of one dimensional temperature field model


新窗口打开
表1为一维温度场模型正则化的绝对误差以及相对误差. 绝对误差由下式计算
$ \left| {e_{\rm abs} } \right| = \sqrt {\sum_{i = 1}^{N_{1} } {e_{{\rm abs},i}^2 } } (21) $
式中, N1代表节点总数, eabs,i为节点 i的绝对误差, eabs,i=Θ*-Θ, Θ*为本文耦合方法计算得到的温度值, Θ为解析解结果.
相对误差的表达式为
$ \left| {e_{\rm rel} } \right| = \sqrt {\sum_{i = 1}^{N_{1} } {e_{{\rm rel},i}^2 } } (22) $
式中, erel,i为节点 i的绝对误差, erel,i=Θ*-Θ/Θ. 只对没有施加边界条件的节点计算误差,从表1可以看出,两种边界条件的绝对误差与相对误差都接近于0.

2.2 二维情况下热传导算例

针对二维热传导模型如图6所示,模型取长 a=100m,宽 b=100m的方板,方板厚度 h=1m. 模型上半部分( k1=2W/m?K)由材料1组成,其热传导系数 y<b/2);下半部分( k2=1W/m?K由材料2组成,其热传导系数为 0.3a#x2264;x#x2264;0.7a. PD区域的微热传导系数由式(5)和式(6)确定. 如图6(a),将模型分为两部分, 其中区域 0.4b#x2264;y#x2264;0.6b, koikii-1kio采用PD粒子描述,其他区域采用FEM节点描述. 本例中,除了取边界节点自由度为主自由度以外,PD区域粒子的自由度也取为主自由度,如图6(b)所示. 这样做的目的在于,在做破坏分析时,破坏只在发生在PD区域,其他节点的热传导矩阵(刚度矩阵)以及节点的热流(载荷矩阵) 在计算过程中不会发生变化,式(16)中的 kii-1kioPi和式(17)中的 koo只需在第一步需要计算,其他计算步骤只需要更新式(16)中的 Po和式(17)中的 δ.
显示原图|下载原图ZIP|生成PPT
图6FEM/PD区域划分及子结构划分示意图
-->Fig. 6FEM/PD region division and substructure division diagram
-->

接下来讨论本文耦合算法的收敛性. 对 m=4收敛性分析,取 Δx=2,分别取 δ=8m ( Δx=1m)、 δ=4m ( Δx=0.5m)和 δ=2m ( δ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))在两种材料界面附近的计算精度比较低.
显示原图|下载原图ZIP|生成PPT
图7不同\delta值时耦合模型温度场模拟结果
-->Fig. 7Simulation results of temperature field of coupling model with different m values
-->

对于 δ=4收敛性分析,取 Δx=2m,分别取 Δx=1m, Δx=0.5m和 Δx=2m,对于 m=2m,每个PD粒子与周围12个粒子以PD方式相互作用 ( Δx=1);对于 m=4m,每个PD粒子与周围48个粒子以PD方式相互作用 ( Δx=0.5);对于 m=8m,每个PD粒子与周围196个粒子以PD方式相互作用( δ). 边界条件仍按式(23)选取. 图8为不同 Δx=1值时温度梯度的误差沿横向的变化曲线. 通过图8可以看出,对 于 Δx=0.5m的情况,本文耦合方案的误差最小. Δx=1m时,误差反而大于 Δx=1m的情况,这是因为与 Δx=0.5m的情况相比, Δx=2m时两种材 料界面附近的PD粒子要与更多另外一种材料中的粒子相互作用,从而增大了误差. mm时,虽然材料界面附近的PD粒子与另外一种材料中较少的粒子相互作用,但其对应的 m值仅为2,致使其误差最大.
显示原图|下载原图ZIP|生成PPT
图8不同m值时耦合模型温度场模拟结果
-->Fig. 8Simulation results of temperature field of coupling model with different \delta values
-->

通过以上 m收敛性分析以及 m收敛性分析,可以得到结论,为了提高材料界面处的计算精度, m值不宜取太大,另一方面,为了保 证PD的计算精度, m值不能太小. 通过上述算例分析, m取4是一个恰当的选择. 在 Δx=1值恰当的前提下,提高网格密度可以提高计算精度. 在以下二维算例中,取 m=4m, δ=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)所示,本文采用第二种形式.
显示原图|下载原图ZIP|生成PPT
图9二维温度场模型模拟结果
-->Fig. 9Temperature field model simulation results of two-dimensional
-->

显示原图|下载原图ZIP|生成PPT
图9二维温度场模型模拟结果(续)
-->Fig. 9Temperature field model simulation results of two-dimensional (continued)
-->

显示原图|下载原图ZIP|生成PPT
图10含裂纹二维模型图
-->Fig. 10Two-dimensional model with a crack
-->

显示原图|下载原图ZIP|生成PPT
图11初始裂纹表示方法
-->Fig. 11Initial crack presentation method
-->

分别采用传统FEM理论、本文耦合方案和PD理论对含裂纹二维温度场进行模拟,采用式(23) (第一类边 界条件)的模拟结果如图12所示,采用式(24) (第二类边界条件)的模拟结果如图13所示. 从图12和图13可以看出,在以上两类边界条件下,裂纹上下表面的温度梯度最大,这是因为裂纹的存在阻断了热量的传播. 由于裂纹的阻隔以及左右边界的散热,热量几乎不能传到裂纹下侧.
显示原图|下载原图ZIP|生成PPT
图12第一类边界条件下含裂纹二维温度场模拟结果
-->Fig. 12Temperature field results of two-dimensional problem with crack under the first kind of boundary condition
-->

显示原图|下载原图ZIP|生成PPT
图13第二类边界条件下含裂纹二维温度场模拟结果
-->Fig. 13Temperature field results of two-dimensional problem with crack under the second kind of boundary condition
-->

对比图9与图12,图13可知,模型中的裂纹对温度场分布造成了很大的影响. 此外,本文耦合方案的结果与传统FEM理论、PD理论的结果相吻合.
需要指出,本文的耦合方案以及PD模拟结果由编制的MATLAB程序实现,程序采用显式计算,对以上两个算例均采用1000个 增量步进行计算,计算时间如表2所示
Table 2
表2
表2计算时间对比
Table 2Calculation of time contrast


新窗口打开
从表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]