近场动力学与有限元方法耦合求解热传导问题1)

引言

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

创新活动方式，激发代表履职激情。常委会建立了全县“双联”工作云平台，及时发布工作动态，宣传推广“双联”工作典型经验，便于人大代表相互学习、相互交流、共同提高，为广大代表依法履职、扎实做好“双联”工作奠定了坚实基础。

虽然PD在模拟裂纹扩展时有很大的优势.但作为非局部理论，PD计算效率远低于有限元法(fi nite element method,FEM).为了发挥PD处理非连续问题以及FEM计算效率高的优势，许多学者在PD与FEM结合方面做了很多研究[16-25].Macek等[12]将模型划分为PD区域与FEM区域，对两个区域重叠区域中的PD物质点与FEM结点采用位移协调的方法，PD物质点的位移可通过FEM结点位移进行插值后确定，但此过程不可逆，故边界条件只能在FEM区域实施.文献[21-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]可以写为

 

式中，Pc和Pd代表节点的热流.

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

 

将式(2)写为离散形式

 

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

 

式中，|ξ|=|x′− x|,κ代表微热传导系数.为了与长程力更好地匹配，κ不取常数，而是取一个与距离有关的量[32].

 

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

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

  

图1 界面附近的PD粒子Fig.1 PD particles close to interface

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

 

式中，l1代表粒子x′与x的距离|ξ|在材料1中的部分(图中实线部分)，其对应的微热传导系数为κ1；l2代表粒子x′与x的距离|ξ|在材料2中的部分(图中虚线部分)，其对应的微热传导系数为κ2.从图2不难看出 |ξ|=l1+l2，κ1 和 κ2 分别由对应材料的热导率k1,k2代入式(5)得到.

勘探开发一体化是精细勘探有别于常规勘探思路的一个重要方面，井位部署中必须以此为原则，统一考虑探井和开发井的位置，为下一步开发做准备。通过反复论证，将探井、开发井和水井纳入同一个井网系统。根据低渗油田的开发经验，参考临区开发方案，本区井网方式采用五点法不规则面积注水井网。探井和开发准备井采用放大的五点法，即800×566m井距，两口相邻探井、两口相邻的开发准备井之间井距为800m，探井和开发准备井之间井距为566m。根据上述部署方案，试验区内共部署探井9口，开发准备井6口。初步预测成熟区内下步开发总井数为50口，其中油井35口，注水井15口，利用老井7口，新钻井43口。

  

图2 一维情况下耦合区域节点图Fig.2 Nodes in coupling region for one-dimensional case

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

2.2.5 民族文化习惯 少数民族地区农民的思想相对保守，对家乡特定的生态环境、生产、生活方式有严重的依赖性。因此，一般来说民族文化习惯是劳动力外出务工的阻力因素，对外出务工意愿有负向影响。从表4看出，受访者大约有19.8%完全保持家乡的生活习惯，有很强的民族文化观念，非常留恋乡土，外出务工的意愿一般非常弱。34.4%的受访者保持部分民族生活习惯，13.5%的受访者没有保持民族生活习惯，这两部分劳动力通常更适应城市的生活和文化，累计比重为47.9%，外出务工意愿强烈。

 

式中

2.2.5 经预处理的空白脂质体溶液 精密量取“2.1”项下空白脂质体 100 μL，加入800 μL甲醇破乳，涡旋5 min混匀，加入50 mg/mL ZnCl2甲醇溶液100 μL，涡旋5 min混匀，室温静置15 min，12 000 r/min离心5 min，取上清液，即得。

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

 

式中，代表粒子 x 所占的空间，Vx′代表粒子x′所占的空间，对于一维情况，Vx=AΔ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所示的一维模型，总体热传导矩阵可写为

 

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

国家为自贸区的建立出台了一系列政策，也为自贸区知识产权保护制度的构建提供了便利。国务院于2013年在《中共中央关于全面深化改革若干重大问题的决定》中指出，国家应通过更为成熟的手段，加强和运用健全的知识产权保护制度，完善技术创新的激励机制，并努力探索建立健全专门知识产权法院；国家对知识产权的保护日益重视，在知识产权核心法律体系中，法定赔偿额的提高以及相关的顶层制度设计与出台均反映了中央对知识产权的重视程度。重庆自贸区作为我国西部自贸区的试验基地，其地位显得尤为重要，因为先行试验基地的成功与知识产权保护有着密切的关联，同时知识产权的保护也为以后的自由贸易区提供着重要借鉴。

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

  

图3 二维情况下耦合区域节点图Fig.3 Nodes in coupling region for two-dimensional case

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

 

式中，ρ代表材料密度，cv为材料比热容，H表示以δ为半径的粒子x的邻域，fh代表热流密度函数，Sb为内热源.t和Θ分别表示时间和温度，代表粒子x′所占的体积.

 

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

 

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

 

其中，kdc和 kdd为 FEM 单元热传导矩阵中的元素[33].

这不何东下了班一进家门，老妈郑玉英就从客厅里飞过来，举一牛皮纸口袋在他眼前晃着说：“都准备好了，户口本身份证还有相片都在这儿了。”

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

将式(10)写为分块格式

 

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

 

式(15)为缩聚后的有限元方程，kc0的阶数与koo的相同，其中

 

式中，N1代表节点总数，eabs,i为节点i的绝对误差，eabs,i=Θ∗−Θ，Θ∗为本文耦合方法计算得到的温度值，Θ为解析解结果.

在控制农业污染方面，要采取工程措施和非工程措施相结合的办法，积极发展喷灌、滴灌、微灌等节水灌溉设施，大力推广符合农艺要求的化肥深施、浅湿灌溉，尽量减少带有农药、化肥的农田污水进入河道；在控制工业污染方面，对排放污水的工厂企业要落实治污措施，污水达标方能排入河道，对排污问题严重的工厂企业要限期整改，确实无能力治理的必须及时关闭；各地必须充分考虑水资源的承载能力，做到工厂企业和污水处理设施同时规划，同时设计，同时建设；在农村居民生活污水方面，要加快收集和处理设施建设，特别是对集中居住区及畜禽养殖场的生产、生活废水加快截污、处理，杜绝直接排放入河。

 

采用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所示.

以粉煤灰“一步酸溶法”提取氧化铝后尾渣—白泥为原料，通过碱液浸取—水热合成工序制备13X分子筛。研究了导向剂用量、老化温度、铝源对13X分子筛的合成的影响，当导向剂用量为原液体积的3%、老化温度为60 ℃、铝源为偏铝酸钠时，合成的13X分子筛具有较好的热稳定性，该分子筛对Pb2+的饱和吸附量达到35.82 mg/g。因此，利用白泥为原料成功制备分子筛，显著降低了分子筛的生产成本，推动了13X分子筛的生产与工业化应用，在污水处理方面有较好的应用前景。

其中Vd为粒子d所占有的体积.N代表PD粒子c作用域内的粒子总数，代表PD方式的作用

  

图4 一维热传导模型Fig.4 One-dimensional heat conduction model

第一类边界条件

 

第二类边界条件

 

其中，q代表热流密度.

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

  

图5 一维温度场模型模拟结果Fig.5 Simulation results for one-dimensional temperature fiel model

 

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

  

Cases Norm of absolute errors Norm of relative errors dirichlet boundary condition 7.84×10−13 4.51×10−14 Neumann boundary 2.59×10−12 1.03×10−13

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

 

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

理论课程：坚持课堂教学是能力培养主阵地，理论教学是能力形成的重要渠道，关注知识传授与能力培养相结合，寓能力培养于知识传授之中，知识传授服务于学生开发智力、培养应用等能力的目的。传授知识，培养学生的理论思维能力及其它能力。

相对误差的表达式为

 

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

2.2 二维情况下热传导算例

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

时至今日，在全国多地仍有佩戴香囊的习俗。人们用彩色丝线在彩绸上绣制出各种图案纹饰，缝制成形状各异、大小不等的小绣囊，内装多种芳香开窍的中草药、干花或香料，以求驱虫、避邪、安眠、保平安、传达思念等。

  

图6 FEM/PD区域划分及子结构划分示意图Fig.6 FEM/PD region division and substructure division diagram

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

 

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

如图7所示,随着处理Pb2+溶液体积由130 mL增加至150 mL，Pb2+吸附量由35.46 mg/g增至35.82 mg/g；此时，吸附量并未随体积呈现线性增长，吸附量趋于平衡，达到最大为35.82 mg/g；相比而言，溶液中Pb2+离子去除率则呈现骤降趋势，由83.12%降至72.16%，也说明分子筛对Pb2+离子的吸附量达到饱和。

2.3.5 SaO2 纳入4篇文献，各研究间存在异质性（P＜0.000 1，I2=90%），采用随机效应模型进行Meta‐分析，见图5。结果显示治疗前后试验组SaO2增加值显著大于对照组，差异有统计学意义［MD=5.32，95%CI（3.14～7.50），P＜0.000 1］。

  

图7 不同δ值时耦合模型温度场模拟结果Fig.7 Simulation results of temperature fiel of coupling model with different δ values

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

刘勰采用的“虚静”这个概念，又不完全与道家的这种思想相同。刘勰将这种澄怀的状态用于文学的创作之中，要求在成文之前保持一种沉着、平静的心理状态，如一湾平静的湖水般的虚空、静谧，冷静观察，积累沉淀、感受一段时间之后再创作，最后一气呵成，而不能太依赖于主观的能动性，从而达到一种自然而然的“道”。

  

图8 不同m值时耦合模型温度场模拟结果Fig.8 Simulation results of temperature fiel of coupling model with different m values

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

 

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

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

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

 

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

  

图9二维温度场模型模拟结果Fig.9 Temperature fiel model simulation results of two-dimensional

  

图9 二维温度场模型模拟结果(续)Fig.9 Temperature fiel model simulation results of two-dimensional(continued)

  

图10 含裂纹二维模型图Fig.10 Two-dimensional model with a crack

  

图11 初始裂纹表示方法Fig.11 Initial crack presentation method

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

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

  

图12 第一类边界条件下含裂纹二维温度场模拟结果Fig.12 Temperature fiel results of two-dimensional problem with crack under the firs kind of boundary condition

  

图13 第二类边界条件下含裂纹二维温度场模拟结果Fig.13 Temperature fiel results of two-dimensional problem with crack under the second kind of boundary condition

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

 

表2 计算时间对比Table 2 Calculation of time contrast

  

Cases Coupling method PD dirichlet boundary condition 113.53s 563.30s Neumann boundary 104.21s 516.10s

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

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

处理器∶i5-4590 CPU@3.30 GHz

内存∶16.00 GB

操作系统∶Windows 10 Enterprise 64 bit

3 结论

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

参考文献

1 马天宝，任会兰，李健等.爆炸与冲击问题的大规模高精度计算.力学学报,2016,48(3)∶599-608(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))

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(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))

5 Silling SA,Lehoucq RB.Peridynamic theory of solid mechanics.Advances in Applied Mechanics,2010,44∶73-168

6 Silling SA.Origin and effect of nonlocality in a composite.Journal of Mechanics of Materials and Structures,2014,9(2)∶245-258

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

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 the-ory.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

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 fibe reinforced concrete structures.Computers&Structures,2015,161∶43-54

16 Macek RW,Silling SA.Peridynamics via finit element analysis.Finite Elements in Analysis and Design,2007,43(15)∶1169-1178

17 Kilic B,Madenci E.Coupling of peridynamic theory and the finit 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 finit 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

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

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

24 Shojaei A,Mudric T,Zaccariotto M,et al.A coupled meshless finit point/Peridynamic method for 2D dynamic fracture analysis.International Journal of Mechanical Sciences,2016,119∶419-431

25 Bie YH,Cui XY,Li ZC.A coupling approach of state-based peridynamics with node-based smoothed finit element method.Computer Methods in Applied Mechanics and Engineering,2018,331∶675-700

26 Bobaru F,Duangpanya M.The peridynamic formulation for transient heat conduction.International Journal of Heat and Mass Transfer,2010,53(19)∶4047-4059

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

28 Oterkus S,Madenci E,Agwai A.Peridynamic thermal diffusion.Journal of Computational Physics,2014,265∶71-96

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

30 王飞,马玉娥,郭妍宁.近场动力学中内核参数对非均匀材料热传导数值解的影响研究.西北工业大学学报,2017,35(2)∶203-207(Wang Fei,Ma Yu’e,Guo Yanning.Effects of kernel parameters of peridynamic theory on heat conduction numerical solution for nonhomogeneous material.Journal of Northwestern Polytechnical University,2017,35(2)∶203-207(in Chinese))

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

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

33 孔祥谦.热应力有限单元法分析.上海：上海交通大学出版社,1999∶136-139(Kong Xiangqian.Analysis of thermal stress finit element method.Shanghai∶Shanghai Jiao Tong University Press,1999∶136-139(in Chinese))

34 邱吉宝,向树红,张正平.计算结构动力学.合肥：中国科学技术大学出版社,2009∶352-354(Qiu Jibao,Xiang Shuhong,Zhang Zhengping.Computational Structure Dynamics.Hefei∶University of Science and Technology of China Press,2009∶352-354(in Chinese))

35 Ha YD,Bobaru F.Studies of dynamic crack propagation and crack branching with peridynamics.International Journal of Fracture,2010,162(1)∶229-244