椭球粒子的本征应变边界积分方程与局部Eshelby矩阵

马　杭，方静波

（1.上海大学理学院，上海　200444；2.上海大学上海市应用数学和力学研究所，上海　200072）

椭球粒子的本征应变边界积分方程与局部Eshelby矩阵

马杭1，方静波2

（1.上海大学理学院，上海200444；2.上海大学上海市应用数学和力学研究所，上海200072）

针对粒子增强材料的大规模数值模拟问题，将局部Eshelby矩阵的概念引入到本征应变边界积分方程计算模型中，以解决粒子间的相互作用问题.局部Eshelby矩阵可以看作Eshelby张量和等效夹杂物的概念在数值方面的一种拓广.以全空间边界元子域法为参照，利用计算模型对无限域中的若干椭球粒子进行了三维应力分析.数值算例不仅验证了模型的正确性和方法的可行性，也表现出较高的计算效率，说明该计算模型和方法具有对粒子增强材料进行大规模数值分析的能力.

本征应变；Eshelby张量；等效夹杂物；近场粒子；局部Eshelby矩阵

自从20世纪50年代Eshelby［1］提出本征应变和等效夹杂物的概念以来，嵌入夹杂物的固体弹性状态成为固体力学领域的研究热点之一［2-5］.本征应变解可对应于非协调热应变、相变应变、塑性应变、等效夹杂物中的虚拟应变、量子点线［6］以及残余应力问题中的固有应变［7］等，用以表达十分广泛的物理力学问题.通过等效夹杂物的替换，本征应变解又可对应于各种异质体［5］、孔洞甚至裂纹问题［8］.因此，研究本征应变解具有重要的意义.

解析方法可以提供粒子内外应力应变分布的精细图像，为进一步研究作基础.然而，解析方法通常仅适于全空间或半空间具有简单几何形状（如椭球、圆柱和球体）问题的求解.有限元［9］以及基于边界元的数值方法［10-12］是分析复杂几何形状和材料问题的工具.对于含有大量粒子的固体，如果粒子具有随机的尺寸分布、空间分布以及各种不同的形状，出于界面的描述和数值离散的要求，数值方法求解的主要困难在于求解的规模太大.解决上述困难的途径之一是引入快速多极展开的特殊技术［13-14］，但这将大大增加算法的复杂性.因此，含有大量粒子固体的数值分析仍然是具有挑战性的课题.

最近，Ma等［15-17］将Eshelby本征应变和等效夹杂物的概念引入边界积分方程，提出了新的计算模型及相应的迭代算法，并对含有大量粒子的固体进行了二维和三维应力及等效性能的分析.对于密集分布的粒子，粒子间的相互作用将影响迭代算法的收敛性.为了解决这一问题，本工作结合Eshelby张量和等效夹杂物的概念，利用全空间近场粒子群的积分方程离散形式，提出了局部Eshelby矩阵的概念，改进了原有的算法，并以全空间边界元子域法［18］为参照，利用本征应变边界积分方程计算模型对无限域中的若干椭球粒子进行了三维应力分析，验证了模型的正确性和方法的可行性.结果表明，该算法具有较高的计算效率.

1　计算模型

1.1本征应变边界积分方程

在计算模型中，设基体与粒子均为各向同性，二者在界面上理想结合，即在界面上位移连续、面力平衡.求解域Ω表示基体，其外边界为Γ，粒子ΩI的边界为ΓI（ΓI=ΩI∩Ω）.位移和应力的本征应变边界积分方程［15-17］为

式中，

为域积分自由项，Ωε（其边界为Γε）为ΩI中的ε小区，xl=xl（q）-xl（p）为场点q和源点p距离的投影；和分别为位移、面力和应力的Kelvin基本解和为相应的基本解导出函数；NI为域Ω中的粒子总数）为本征应变.事实上，边界积分方程（1）和（2）仅仅描述了均质弹性体的位移场和应力场，而要描述含有非均质粒子的固体，必须借助Eshelby张量的概念进行等效夹杂物的替换.

1.2Eshelby张量和等效夹杂物

根据Eshelby［1］的工作，无应力条件下全空间内单个粒子的拘束应变与本征应变）

可以通过Eshelby张量Sijkl联系起来，即

式中，Eshelby张量依赖于ΩI的几何形状.对于简单的形状，Eshelby张量可表达为显式，而在一般情况下则表达为积分的形式［17］.对于表达为积分形式的Eshelby张量，通过数值方法可求得

式中，µ为基体的剪切模量，v为Poisson比.应指出的是，式（5）适用于本征应变在域ΩI内均匀分布的情况.对于多粒子问题，本征应变均匀分布的适用条件是粒子的体积足够小或者粒子之间保持适当的距离.因此，本工作采用了本征应变均匀分布的假定.定义杨氏模量比βI=EI/EM，其中下标I和M分别表示第I个粒子和基体.根据Hooke定律，若用等效夹杂物替换外加应变εij作用下的粒子而不改变该处的应力状态，则应满足的关系式为

式中，

结合式（4）和（6），可利用外加应变εij计算单个粒子的本征应变，由此实现等效夹杂物的替换，则含有粒子的固体成为了形式上的均质体，可由本征应变边界积分方程（1）和（2）来描述.对于单个粒子时的情况，式（6）离散后的矩阵形式为

1.3局部Eshelby矩阵

如前所述，在进行等效夹杂物替换时，需要利用外加应变来确定粒子的本征应变.然而，对于多个粒子的情况，粒子之间存在相互作用，其作用的大小与粒子间的距离有关.除了粒子本身的状态，粒子上的外加应变还受到周边其他粒子的干扰作用，从而影响算法的迭代收敛性［15-17］.如图1所示，对于当前粒子I，定义虚线圆内的粒子为近场粒子，其数量为NL，虚线圆外的粒子则为远场粒子.对于当前粒子来说，近场粒子属于近距离粒子，其相互作用较为强烈，而远场粒子对于当前粒子的影响则相对较弱.

图1　多粒子群组的定义Fig.1 Group definitions for multiple particles

对于全空间，如果只考虑近场粒子的影响而暂时忽略远场粒子，结合本构关系以及积分类型的转换关系［17］，当前粒子I的应力可用应力边界积分方程（2）表达，即为

式中，Eijkl为弹性张量.由于全空间不存在外边界，式（2）中的2个边界积分在式（9）中并不出现.将Eshelby张量的概念推广，并利用式（5）和（9），则当前粒子的拘束应变以及Eshelby张量可分别表达为

式中，当上标I/=J时，Eshelby张量表示近场粒子中不同粒子的本征应变与拘束应变之间的关系.与式（6）的替换关系类似，根据Hooke定律，近场粒子的等效夹杂物替换所应满足的关系式为

结合式（10）和（11），式（12）离散后的矩阵形式为

由此即可利用近场粒子的外加应变向量来计算其本征应变向量.与文献［8］中对裂纹的分析类似，本工作将式（13）中的矩阵称为Eshelby矩阵，该矩阵可以看作Eshelby张量和等效夹杂物的概念在数值方面的一种拓广.应当指出的是：式（9）～（12）是解析式，反映了全空间近场粒子之间应力应变的特定关系，是对近场粒子状态的准确描述；而式（13）则是离散后的数值形式，是对近场粒子的状态和关系的近似描述.

对于具体的问题，如果指定了虚线圆（见图1）的半径，一般情况下固体中的每个粒子都会与一个独特的近场粒子群相对应.为了计算的方便，将式（13）改写为式中，SI是对Eshelby矩阵求逆和缩并得到的，且与当前粒子的本征应变向量相对应；对应于当前粒子的近场粒子的外加应变向量，与式（13）中相同；下标I遍历固体中的所有粒子，其总数为NI.

这样，基于本工作提出的本征应变边界积分方程与局部Eshelby矩阵的算法（简称Eigen），固体中每个粒子本征应变的计算均由3部分构成：第1部分为外加载荷作用下产生的各个粒子上的应力或外加应变；第2部分为近场粒子群的作用，其相互影响较为强烈，利用式（14）进行计算；第3部分为远场粒子群对当前粒子的作用，其影响相对较弱，可直接利用积分方程中的域积分来计算.事实上，其他粒子对当前粒子的作用大小和强弱是由无量纲距离来决定的，由此无量纲距离可作为粒子群中的近场和远场粒子定义的依据.显然，近场粒子群的最小数量为1，如果NL=1，本工作算法就退化为文献［15-17］中提出的算法，这时式（14）简化为式（8）.

2　全空间边界元子域法

本工作的主要目的是验证所提出计算模型的正确性和方法的可行性，因此选择边界元子域法（boundary element method，BEM）作为对照的基准.文献［18］中对全平面中的2个粒子进行了二维分析，但由于全空间没有外边界，其加载条件具有一定的特殊性，因此有必要对全空间边界元子域法（见图2）作简要的介绍.

图2　全空间边界元子域法的示意图Fig.2 Schematics of the sub-domain BEM in full space

如图2（a）所示，考虑远场载荷作用下全空间中的任一粒子ΩI，其边界记作ΓI.现在保持界面状态不变，将图2（a）分解成内域（图2（b））和外域（图2（c））两种情况：图2（b）中粒子的边界上存在位移uI与面力τI，但无载荷作用；而图2（c）中孔洞的边界ΓI上，除了位移uI与面力—τI以外，还存在远场载荷.注意内域和外域之间满足界面条件，即位移连续、面力平衡，因此二者面力的符号相反.由于全空间不存在外边界，为了使边界积分方程能够表达远场载荷的作用，须将图2（c）进一步分解为图2（d）和图2（e），其中图2（d）所示的虚拟界面（虚线）上的位移与面力-容易计算，图2（e）所示的孔洞则适于全空间边界积分方程的描述.

按照图2所示的方法分解，对于含有NI个粒子的全空间，每个粒子都需要单独进行描述，即需要NI个边界积分方程，而分解后全空间中的所有孔洞（图2（e）中只画出了1个孔洞）则用另一个整体的边界积分方程来描述.通过界面条件将这NI+1个方程联系起来形成积分方程组，离散后求解，即所谓全空间边界元子域法.边界积分方程为

注意，式（15）中不含域积分.对于每个粒子，边界积分方程离散后的矩阵形式为

式中，UI和TI分别为位移与面力基本解形成的系数矩阵.如果矩阵UI可逆，ΓI上的面力向量τI可表达为

对于全空间中的NI个孔洞，边界积分方程离散后的矩阵形式为

式中，UIJ和TIJ分别为矩阵U和T中的子矩阵元素，其中下标I和J分别表示源点p和场点q位于ΓI和ΓJ上.当I=J时，子矩阵元素位于式（19）中“=”左边矩阵的主对角元上，描述的是同一个孔洞.事实上，根据积分核的性质，下列关系式成立：

式（21）成立的前提是假设粒子与基体具有相同的Poisson比，其中βI=EI/EM为粒子与基体的杨氏模量比.利用式（19）求得界面位移u后，再利用式（17）依次求出每个粒子界面上的面力τI.

3　算例

本工作的主要目地是检验计算模型，即Eshelby矩阵的正确性.假定粒子上的本征应变均匀分布，且只对全空间含有少量椭球粒子的情况进行数值计算，则在计算中只选择一个近场粒子群即可，其近场粒子群中粒子的数量与粒子的总数相同，即NL=NI.因此，本工作中的算例并不需要进行迭代计算.如果采用文献［15-17］中给出的算法，即NL=1，则必须进行迭代计算.图3为椭球的轴和第一卦限界面的单元剖分，其中图3（b）所示的界面剖分用于计算粒子的Eshelby张量与Eshelby矩阵，或作为边界元子域法的离散单元.

图3　椭球的轴和第一卦限界面的单元剖分Fig.3 Definition of the axes of ellipsoid and the boundary elements in one octant

3.12个粒子的应力分布

在远场z向单位载荷作用下，全空间中的2个水平排列的扁椭球与其中心连线上的应力分布如图4所示.分别采用本征应变边界积分方程模型和全空间边界元子域法，计算了2个扁椭球中心连线上的应力，其中软粒子的模量比为EI/EM=0.01，硬粒子的模量比为EI/EM=10.从计算结果（见图4（b））可以看出，采用2种方法得到的结果十分吻合，说明了Eshelby矩阵的正确性与算法的可行性.另外，由图4（b）可知，无论是软粒子还是硬粒子，其应力分量σ33在界面上都出现了跳跃.

图4　2个水平排列的扁椭球粒子中心连线上的应力分布Fig.4 Two horizontally placed oblate ellipsoidal particles and the stress distributions

图5为全空间中水平排列的2个不同弹性模量的球形粒子和其中心连线上的应力分布.图5（a）中左侧为软粒子，右侧为硬粒子.由图5（b）可以看出，2种方法得到的结果吻合良好.

图5　2个水平排列的不同弹性模量球形粒子中心连线上的应力分布Fig.5 Two horizontally placed spheroidal particles with different modulus and the stress distributions

图6为全空间中垂直排列的2个球形粒子和其中心连线上的应力分布.图6（b）所示的结果同样说明2种方法的结果吻合良好.另外，由图6（b）可以看出，在2个球形粒子垂直排列的情况下，软粒子和硬粒子的应力分量σ11都在界面上出现了跳跃.

图6　2个垂直排列的球形粒子中心连线上的应力分布Fig.6 Two vertically placed spheroidal particles and the stress distributions

3.23个粒子的应力分布

在远场z向单位载荷作用下，全空间中的3个水平排列的长椭球粒子如图7（a）所示，分别采用本征应变边界积分方程模型和全空间边界元子域法计算了右侧2个长椭球粒子中心连线上的应力，其中软粒子的模量比为EI/EM=0.01，硬粒子的模量比为EI/EM=10.计算结果（见图7（b））依然表明2种方法得到的结果十分吻合，说明了Eshelby矩阵的正确性与算法的可行性.同样，由图7（b）可知，无论是软粒子还是硬粒子，其应力分量σ33在界面上都出现了跳跃，但是2个长椭球粒子的状态并不完全相同，另外，2个长椭球粒子内部的应力也是有区别的.

图7　3个水平排列的长椭球粒子中右侧2个椭球粒子中心连线上的应力分布Fig.7 Three horizontally placed prolate ellipsoidal particles and the stress distributions

图8为全空间中水平排列的3个球形粒子和右侧2个粒子中心连线上的应力分布.图8（b）所示的结果也同样说明2种方法得到的结果吻合良好.比较图8（b）与图7（b）的结果，可以看出粒子几何形状对应力分布的影响.对于硬粒子，球形粒子内部的应力分量σ33高于长椭球粒子，而基体的应力分量σ33有所降低；对于软粒子，基体的应力分量σ33有所上升.

图8　3个水平排列的球形粒子中右侧2个粒子中心连线上的应力分布Fig.8 Three horizontally placed spheroidal particles and the stress distributions

3.35个粒子的应力分布

在远场z向单位载荷作用下，全空间中的5个在竖直平面内排列的球形粒子如图9（a）所示，分别采用本征应变边界积分方程模型和全空间边界元子域法计算了中间粒子和右侧粒子二者中心连线上的应力，其中软粒子的模量比为EI/EM=0.01，硬粒子的模量比为EI/EM=10.计算结果（见图9（b））依然表明2种方法得到的结果十分吻合，说明了Eshelby矩阵的正确性与算法的可行性.同样，由图9（b）可知，无论是软粒子还是硬粒子，其应力分量σ33在界面上都出现了跳跃.另外，由于中间粒子和右侧粒子的状态有差别，图9中2个粒子内部应力的差别比图8中2个粒子差别更大.

图9　5个球形粒子中的中间粒子与右侧粒子中心连线上的应力分布Fig.9 Horizontal computing line for five spheroidal particles and the stress distributions

图10（b）为5个球形粒子中间粒子的和上方粒子二者中心连线上的应力分布.图10（b）的结果同样说明2种方法的结果吻合良好.比较图10与图6，可以看出粒子排布对应力分布的影响.

图10　5个球形粒子中的中间粒子与上侧粒子中心连线上的应力分布Fig.10 Vertical computing line for five spheroidal particles and the stress distributions

众所周知，含有粒子的介质属于异质体问题，界面上的应力分量有些是连续的，而有些会出现跳跃，例如软粒子（EI/EM=0.01）的应力分量σ33.采用边界元子域法时，界面相当于子域的边界，因此界面应力的计算需要分别在两侧的边界上同时进行.应当指出的是，界面上的应力分量出现跳跃恰恰反映了界面两侧材料具有不同弹性模量的影响.采用本征应变边界积分方程模型时，由于等效夹杂物的替换，计算区域成为表观均质体，因此界面应力的计算可以直接在界面上进行.计算结果表明，对于应力分量的数值出现跳跃的情况，按本工作中的算法得到的界面应力分量的数值取两侧数值的平均值，如图4～9中的σ33，图6和10中的σ11，σ22等，这一结果与前期工作得到的结果［17］相同.

3.4计算效率

本工作提出的本征应变边界积分方程模型的计算效率与全空间边界元子域法相比有较大的差别.2种算法的平均CPU时间如表1所示.由表1可以看出，2种算法的CPU时间相差达2个数量级，其中本征应变边界积分方程模型具有较高的计算效率.另外，由表1中2种算法CPU时间的比值可以看出，随着粒子数量的增加，2种算法计算效率的差别将增大.这是因为全空间边界元子域法的计算时间主要取决于总体系数矩阵的求解，其矩阵规模随着粒子的数量而增大，因此CPU时间按几何级数的方式增长；而在本征应变边界积分方程的算法中，当指定了近场粒子群的大小之后，其CPU时间大体上按算术级数的方式增长.

表1　BEM和Eigen算法的平均CPU时间Table 1 Mean CPU times of BEM and Eigen

4　结束语

针对具有密集分布粒子的粒子增强材料的大规模数值模拟问题，本工作利用全空间近场粒子群的积分方程离散形式，提出了局部Eshelby矩阵的概念，这一概念可以看作Eshelby张量和等效夹杂物的概念在数值方面的一种拓广.通过将局部Eshelby矩阵引入到本征应变边界积分方程计算模型中，改进了原有的算法，解决了近场粒子间的相互作用问题，并且该问题将影响算法的迭代收敛性.以边界元子域法为参照，利用本征应变边界积分方程计算模型对无限域中的若干椭球粒子和球形粒子进行了三维应力分析，验证了模型的正确性、方法的可行性.结果表明，Eigen模型具有较高的计算效率，具有对粒子增强材料大规模数值分析的能力.

作为初步工作，粒子中的本征应变被假定为常数，这是本工作中计算模型的局限所在.进一步的工作需要选取适当的插值方式，取消本征应变的常数假定，从而将计算模型推广到一般几何形状粒子的分析中.

［1］Eshelby J D.The determination of the elastic field of an ellipsoidal inclusion and related problems［J］.Proceedings of the Royal Society of London Series A：Mathematical and Physical Sciences，1957，A241（1226）：376-396.

［2］Mura T，Shodja H M，Hirose Y.Inclusion problems（part 2）［J］.Applied Mechanics Review，1996，49（10）：S118-S127.

［3］Kiris A，Inan E.Eshelby tensors for a spherical inclusion in microelongated elastic fields［J］. International Journal of Engineering Science，2005，43：49-58.

［4］Mercier S，Jacques N，Molinari A.Validation of an interaction law for the Eshelby inclusion problem in elasto-viscoplasticity［J］.International Journal of Solids and Structures，2005，42：1923-1941.

［5］Shen L X，Yi S.An effective inclusion model for effective moluli of heterogeneous materials with ellipsoidal inhomogeneities［J］.International Journal of Solids and Structures，2001，38：5789-5805.

［6］Pan E.Elastic or piezoelastic fields around a quantum dot：fully coupled or semicoupled model？［J］Journal of Applied Physics，2002，91（6）：3785-3796.

［7］Ma H，Deng H L.Nondestructive determination of welding residual stresses by boundary element method［J］.Advances in Engineering Software，1998，29：89-95.

［8］Ma H，Guo Z，Dhanasekar M，et al.Efficient solution of multiple cracks in great number using eigen COD boundary integral equations with iteration procedure［J］.Engineering Analysis with Boundary Elements，2013，37（3）：487-500.

［9］Kakavas P A，Kontoni D N.Numerical investigation of the stress field of particulate reinforced polymeric composites subjected to tension［J］.International Journal for Numerical Methods in Engineering，2006，65：1145-1164.

［10］Lee J，Choi S，Mal A.Stress analysis of an unbounded elastic solid with orthotropic inclusions and voids using a new integral equation technique［J］.International Journal of Solids and Structures，2001，38：2789-2802.

［11］Dong C Y，Cheung Y K，Lo S H.A regularized domain integral formulation for inclusion problems of various shapes by equivalent inclusion method［J］.Computer Methods in Applied Mechanics and Engineering，2002，191（31）：3411-3421.

［12］Nakasone Y，Nishiyama H，Nojiri T.Numerical equivalent inclusion method：a new computational method for analyzing stress fields in and around inclusions of various shapes［J］. Materials Science and Engineering，2000，A285：229-238.

［13］Greengard L F，Rokhlin V.A fast algorithm for particle simulations［J］.Journal of Computational Physics，1987，73：325-48.

［14］Liu Y J，Nishimura N，Tanahashi T，et al.A fast boundary element method for the analysis of fiber-reinforced composites based on a rigid-inclusion model［J］.ASME Journal of Applied Mechanics，2005，72：115-128.

［15］Ma H，Yan C，Qin Q H.Eigenstrain formulation of boundary integral equations for modeling particle-reinforced composites［J］.Engineering Analysis with Boundary Elements，2009，33（3）：410-419.

［16］Ma H，Xia L W，Qin Q H.Computational model for short-fiber composites with eigen-strain formulation of boundary integral equations［J］.Applied Mathematics and Mechanics，2008，29（6）：757-767.

［17］Ma H，Fang J B，Qin Q H.Simulation of ellipsoidal particle-reinforced materials with eigenstrain formulation of 3D BIE［J］.Advances in Engineering Software，2011，42（10）：750-759.

［18］Chen Y Z.Boundary integral equation method for two dissimilar elastic inclusions in an infinite plate［J］.Engineering Analysis with Boundary Elements，2012，36（1）：137-146.本文彩色版可登陆本刊网站查询：http：//www.journal.shu.edu.cn

Eigenstrain boundary integral equation with local Eshelby matrix for ellipsoidal particles

MA Hang1，FANG Jing-bo2
（1.College of Sciences，Shanghai University，Shanghai 200444，China；2.Shanghai Institute of Applied Mathematics and Mechanics，Shanghai University，Shanghai 200072，China）

Aiming at large scale numerical simulation of particle reinforced materials，a concept of local Eshelby matrix is introduced into a computational model of the eigenstrain boundary integral equation to solve the problem of interactions among particles.The local Eshelby matrix can be considered as an extension of Eshelby tensor and an equivalent inclusion in a numerical form.Taking the sub-domain boundary element method as the control，three-dimensional stress analyses are carried out for some ellipsoidal particles in infinite media with the proposed computational model.Numerical examples verify correctness，feasibility and high efficiency of the present model with the corresponding solution procedure，showing potential of solving large scale numerical simulations for particle reinforced materials..

eigenstrain；Eshelby tensor；equivalent inclusion；near field particle；local Eshelby matrix

O 241

A

1007-2861（2015）03-0344-12

10.3969/j.issn.1007-2861.2014.01.039

2013-12-24

国家自然科学基金资助项目（11272195，11332005）

马杭（1951—），男，教授，博士生导师，研究方向为计算固体力学.E-mail：hangma@staff.shu.edu.cn