马腾飞，钱松荣，石宏顺，原群盛

(贵州大学 机械工程学院，贵阳　550025)



引用格式：马腾飞，钱松荣，石宏顺，等.基于Peridynamics理论的材料可靠性分析数值方法的优化[J].重庆理工大学学报(自然科学版)，2016(1):32-36.

Citation format：MA Teng-fei，QIAN Song-rong，SHI Hong-shun，et al.Optimization Method for Reliability Analysis of Concrete for the Theory of Peridynamics[J].Journal of Chongqing University of Technology(Natural Science)，2016(1):32-36.

基于Peridynamics理论的材料可靠性分析数值方法的优化

马腾飞，钱松荣，石宏顺，原群盛

(贵州大学 机械工程学院，贵阳550025)

摘要：Peridynamics理论(近场动力学，PD理论)是研究材料破坏变形领域的一门新理论。该理论通过建立相应的数学模型，用积分的方式描述材料的断裂变形。对PD理论及其传统数值方法作了简要的介绍，基于前人针对PD理论提出的积分方法，引入了另一种全新的固定高斯点积分方法。通过对该方法的应用，PD理论数值计算中的效率和精度都有一定改善。

关键词：近场动力学；高斯点积分；优化；精度

Peridynamics理论(PD理论)是2000年美国Sandia国家实验室的Silling提出的[1-3]。物体中任意两个物质点之间由于万有引力的作用能够产生相互作用力。如图1所示，该理论基于以上思想，假设有一物质点xi，在xi半径为δ的影响域R内，物质点xj对xi有相互作用力f，这个力称为xi点所受到的PD力。通过建立相应的数学模型，把i点在影响域R内所受到的PD力积分起来，再根据牛顿第二定律，给出PD理论的基本运动方程：

(1)

(2)

其中：c是材料的弹性模量，

(3)

s是材料的伸长率，

(4)

则

(5)

其中：s0机械材料的断裂极限；K是材料的体积模量。可以看出： 该方法用积分的方式描述机械材料的运动变形，避免了连续性假设和空间微分方程在不连续状态下求解的奇异性，在解决固体力学断裂变形中的不连续问题时具有一定的优势。由于积分在近场动力学的数值实现中有很重要的作用，因此找到合适的积分方法，使得该理论在计算变形破坏有一个合适的误差精度和一个可接受的计算效率显得尤为重要。

图1　xi点的近场区域

1PD理论常用的积分方法简介

1.1离散积分法

由于PD理论积分方法在数值实现中的重要性，Silling提出了把材料离散的积分方法(cubic-cellintergration,CCI)[3-5]。该方法把材料离散成边长为Δx等大的立方体晶格，如图2所示。影响域R内对i节点PD力的积分，就转化为R内各节点对i节点PD力的求和。此方法由于边界点的计算问题而存在很大误差。如图2所示，i节点的实际影响域为以i为圆心以δ半径为圆形区域，但此方法在计算时，只计算节点在R半径内的立方体晶格，用实心点表示。如图2中的晶格2，阴影区域外的一部分也被计算在内，而晶格1本该计算在内的阴影区域内的一部分却被遗漏掉。

图2　CCI方法的晶格的计数

1.2自适应积分法

针对Silling CCI方法计算精度低的问题，2010年，B.Kilic提出了一种可以误差控制的自适应(AI)积分方法，修正了相关点的计算[6-7]，但要对三维坐标系中每个边界相关点做几何组态划分，计算每一个节点的相关体积，导致实现过程很复杂。在该方法中，k既是i的节点又是j的节点，i和j节点影响域与在k节点上的交线分别是AB和CD，那么i和j点有不同的积分边界，在AI方法中，此处k就有2组不同的梯形点，如图3所示。另外，PD力积分需要积分点的最近位移信息，然而由运动方程决定的大量梯形积分点的位移需要每一步单独计算一次，这都使得该方法计算效率较低。

图3　各节点积分区域在k节点上不同的边界

为了避免不同积分极限带来的麻烦，且在可以接受的效率内保证精度，本研究引入固定高斯点计算方法。

2固定高斯点方法

高斯(Gauss)积分法是数值分析积分中常用的一种近似求积方法[8-9]。Gauss积分公式为

(6)

表1　位置系数权重取值

正确找到高斯点位置后，在坐标中进一步划分相关节点，使得一个高斯点都位于一个子节点的中心。如图4所示，“×”代表高斯点，图4(a)中将一个晶格划分为4个子晶格，图4(b)中将一个晶格划分为9个子晶格。

图4　j节点的高斯子节点划分

在三维坐标中子晶格的边长为

(7)

则子晶格体积为

(8)

将高斯点子节点应用到式(1)得

(9)

(10)

于是，式(9)就变为

f(ug-ui,xg-xi)Vg+bi

(11)

可以看出：本构力的积分公式就转换成了求与一系列高斯点子节点作用力的和。

3面力密度

假设有一均匀变形的物体R，在R内取一点i，让通过i的单位矢量n的法平面将R分成两部分R+和R-。令：

(12)

则i的n方向上的点的面力密度τ(i,n)为

(13)

4计算实例

现假设有一个均匀的各项同性的微弹性脆性材料的弹性模量为1.0×105N/mm2，有效积分区域为δ，泊松比为0.25，极限拉伸系数[10]s0=0.001 2。这个微弹材料受到一个Y方向上的应变矩阵：

由弹塑性力学应力应变转换公式[11]得

(14)

此处取E=1.0×105N/mm2。由弹塑性力学知识可得其应力矩阵如下：

根据式(11)将此微弹材料进行数学建模，以i点为所求面力的中心点，i的影响域为δ，将δ均分m段，令m=12，即Δx=δ/12。取式(6)中n=2，即在三维情况下把1个晶格划分为8个高斯点子晶格，那么R+区域就离散为很多含有一个高斯点xg的小晶格，相对于i点的体力积分就转换为许多xg晶格相对于中心点i的体力的求和。

由弹塑性力学知识知：

(15)

其中u1,u2,u3分别为点xg相对于i点的3个方向上的位移，即式(11)中的ug～ui，x,y,z分别为3个方向上xg点相对i点的相对位置，为式(11)中的xg～xi。

依次改变Δx的大小，编程计算得到应力张量，与理论值对比，得出误差曲线随m值的变化，如图4所示。

图5　精度与Δx值之间关系(n=2)

取式(6)中n=3，即在三维坐标中把1个晶格划分为27个高斯子晶格，仍取Δx=δ/12，编程计算得到相应的应力张量为

5结论

1) 将通过编程计算得到的应力矩阵与实际的应力矩阵对比可知：PD理论的固定高斯点方法计算所得的结果与实际数据相符(误差为4.09%)，证明该方法是可行的。

2) 通过改变m值可以看出：对应力的计算精度随Δx的减小而增大。当m=8时，计算精度开始趋向稳定；但是随着m值增大，影响域R内的相关节点增多，增加了计算量，降低了计算效率。

3) 当高斯点精细时，随着子节点的增多，计算精度得到明显提高(当n=3，m=12时，误差为1.05%)。

4) 当m值相同时，固定高斯点方法的计算精度高于Silling提出的CCI方法的精度[3]，但计算时间略长。对比AI方法的实现过程及精度[6]可知：固定高斯点方法的精度略低于AI方法,但实现过程大为简化，计算效率有明显提高。

通过以上的计算可以看出：固定高斯点积分法是PD理论运动方程数值求解方法的改进，为PD理论的数值积分的实现，为对精度和效率有不同需求的情况提供了一条切实可行的途径。

参考文献：

[1]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.

[2]SILLING S A,ASKARI E.Peridynamic modeling of impact damage[C]//Proceeding of the ASME/JSME Pressure Vessels and Piping Conference.San Diego:[s.n.],2000:197-205.

[3]SILLING S A,ASKARI E.A meshfree method based on the peridynamic model of solid mechanics[J].Computers and Structures,2005,83:1526-1535.

[4]SILLING S A. Dynamic fracture modeling with a meshfree peridynamic code.Second MIT Conference on Computional Fluid and Solid Mechanics.Massachusettes：MIT,2003:641-644.

[5]MACEK R W,SILLING S A.Peridynamics via finite element analysis[J].Finite Elements in Analysis and Design,2007,43(15):1169-1178.

[6]YU K,XIN X J,LEASE K B.A new method of adaptive integration with error control for bond-based peridynamics[C]//Proceeding of the World Congress on Engineering and Computer Science.San Francisco:[s.n.],2010:1041-1046.

[7]YU K,XIN X J，LEASE K B.A new adaptive integration method for the peridynamic theory[J].Modeling and Simulation in Materials Science and Engineering,2011,19:45003.

[8]张民选,罗贤兵.数值分析[M].南京：南京大学出版社,2013.

[9]蔡大用.数值分析与实验学习指导[M].北京：清华大学出版社，2002.

[10]沈峰，章青，黄丹，等.基于近场动力学理论的混凝土轴拉破坏过程模拟[J].计算力学学报,.2013,30(6):79-83.

[11]陈明祥.弹塑性力学[M].武汉：科学出版社，2010.

(责任编辑刘舸)

Optimization Method for Reliability Analysis of

Concrete for the Theory of Peridynamics

MA Teng-fei，QIAN Song-rong，SHI Hong-shun，YUAN Qun-sheng

(School of Mechanical Engineering, Guizhou University, Guiyang 550025，China)

Abstract:Peridynamics (PD) is an emerging new theory used in the research of material damage. This theory describes the fracture and deformation of material by establishing corresponding mathematical model and using integral method. This paper simply introduced the theory of peridynamics and traditional numerical method and a new integration method with fixed Gaussian points was introduced based on previous methods. The efficiency and accuracy of the numerical calculation of PD theory also have relative improvement through the application of the new method.

Key words:peridynamics; Gauss integral; optimization; accuracy

文章编号：1674-8425(2016)01-0032-05

中图分类号：TH115

文献标识码：A

doi:10.3969/j.issn.1674-8425(z).2016.01.006

作者简介：马腾飞(1989—)，男，河南商丘人，硕士研究生，主要从事机械材料的可靠性分析研究。

基金项目：贵州省国际合作项目(黔科合外G字[2013]7006号)；贵州省联合 (黔科合LH字[2014]7624号)；贵州省留学人员科技活动项目 (黔人项目资助合同[2014]13号)

收稿日期：2015-10-23