橢球粒子的本征應(yīng)變邊界積分方程與局部Eshelby矩陣

馬　杭，方靜波

（1.上海大學(xué)理學(xué)院，上?！?00444；2.上海大學(xué)上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所，上?！?00072）

橢球粒子的本征應(yīng)變邊界積分方程與局部Eshelby矩陣

馬杭1，方靜波2

（1.上海大學(xué)理學(xué)院，上海200444；2.上海大學(xué)上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所，上海200072）

針對粒子增強材料的大規(guī)模數(shù)值模擬問題，將局部Eshelby矩陣的概念引入到本征應(yīng)變邊界積分方程計算模型中，以解決粒子間的相互作用問題.局部Eshelby矩陣可以看作Eshelby張量和等效夾雜物的概念在數(shù)值方面的一種拓廣.以全空間邊界元子域法為參照，利用計算模型對無限域中的若干橢球粒子進行了三維應(yīng)力分析.數(shù)值算例不僅驗證了模型的正確性和方法的可行性，也表現(xiàn)出較高的計算效率，說明該計算模型和方法具有對粒子增強材料進行大規(guī)模數(shù)值分析的能力.

本征應(yīng)變；Eshelby張量；等效夾雜物；近場粒子；局部Eshelby矩陣

自從20世紀(jì)50年代Eshelby［1］提出本征應(yīng)變和等效夾雜物的概念以來，嵌入夾雜物的固體彈性狀態(tài)成為固體力學(xué)領(lǐng)域的研究熱點之一［2-5］.本征應(yīng)變解可對應(yīng)于非協(xié)調(diào)熱應(yīng)變、相變應(yīng)變、塑性應(yīng)變、等效夾雜物中的虛擬應(yīng)變、量子點線［6］以及殘余應(yīng)力問題中的固有應(yīng)變［7］等，用以表達十分廣泛的物理力學(xué)問題.通過等效夾雜物的替換，本征應(yīng)變解又可對應(yīng)于各種異質(zhì)體［5］、孔洞甚至裂紋問題［8］.因此，研究本征應(yīng)變解具有重要的意義.

解析方法可以提供粒子內(nèi)外應(yīng)力應(yīng)變分布的精細圖像，為進一步研究作基礎(chǔ).然而，解析方法通常僅適于全空間或半空間具有簡單幾何形狀（如橢球、圓柱和球體）問題的求解.有限元［9］以及基于邊界元的數(shù)值方法［10-12］是分析復(fù)雜幾何形狀和材料問題的工具.對于含有大量粒子的固體，如果粒子具有隨機的尺寸分布、空間分布以及各種不同的形狀，出于界面的描述和數(shù)值離散的要求，數(shù)值方法求解的主要困難在于求解的規(guī)模太大.解決上述困難的途徑之一是引入快速多極展開的特殊技術(shù)［13-14］，但這將大大增加算法的復(fù)雜性.因此，含有大量粒子固體的數(shù)值分析仍然是具有挑戰(zhàn)性的課題.

最近，Ma等［15-17］將Eshelby本征應(yīng)變和等效夾雜物的概念引入邊界積分方程，提出了新的計算模型及相應(yīng)的迭代算法，并對含有大量粒子的固體進行了二維和三維應(yīng)力及等效性能的分析.對于密集分布的粒子，粒子間的相互作用將影響迭代算法的收斂性.為了解決這一問題，本工作結(jié)合Eshelby張量和等效夾雜物的概念，利用全空間近場粒子群的積分方程離散形式，提出了局部Eshelby矩陣的概念，改進了原有的算法，并以全空間邊界元子域法［18］為參照，利用本征應(yīng)變邊界積分方程計算模型對無限域中的若干橢球粒子進行了三維應(yīng)力分析，驗證了模型的正確性和方法的可行性.結(jié)果表明，該算法具有較高的計算效率.

1　計算模型

1.1本征應(yīng)變邊界積分方程

在計算模型中，設(shè)基體與粒子均為各向同性，二者在界面上理想結(jié)合，即在界面上位移連續(xù)、面力平衡.求解域?表示基體，其外邊界為Γ，粒子?I的邊界為ΓI（ΓI=?I∩?）.位移和應(yīng)力的本征應(yīng)變邊界積分方程［15-17］為

式中，

為域積分自由項，?ε（其邊界為Γε）為?I中的ε小區(qū)，xl=xl（q）-xl（p）為場點q和源點p距離的投影；和分別為位移、面力和應(yīng)力的Kelvin基本解和為相應(yīng)的基本解導(dǎo)出函數(shù)；NI為域?中的粒子總數(shù)）為本征應(yīng)變.事實上，邊界積分方程（1）和（2）僅僅描述了均質(zhì)彈性體的位移場和應(yīng)力場，而要描述含有非均質(zhì)粒子的固體，必須借助Eshelby張量的概念進行等效夾雜物的替換.

1.2Eshelby張量和等效夾雜物

根據(jù)Eshelby［1］的工作，無應(yīng)力條件下全空間內(nèi)單個粒子的拘束應(yīng)變與本征應(yīng)變）

可以通過Eshelby張量Sijkl聯(lián)系起來，即

式中，Eshelby張量依賴于?I的幾何形狀.對于簡單的形狀，Eshelby張量可表達為顯式，而在一般情況下則表達為積分的形式［17］.對于表達為積分形式的Eshelby張量，通過數(shù)值方法可求得

式中，μ為基體的剪切模量，v為Poisson比.應(yīng)指出的是，式（5）適用于本征應(yīng)變在域?I內(nèi)均勻分布的情況.對于多粒子問題，本征應(yīng)變均勻分布的適用條件是粒子的體積足夠小或者粒子之間保持適當(dāng)?shù)木嚯x.因此，本工作采用了本征應(yīng)變均勻分布的假定.定義楊氏模量比βI=EI/EM，其中下標(biāo)I和M分別表示第I個粒子和基體.根據(jù)Hooke定律，若用等效夾雜物替換外加應(yīng)變εij作用下的粒子而不改變該處的應(yīng)力狀態(tài)，則應(yīng)滿足的關(guān)系式為

式中，

結(jié)合式（4）和（6），可利用外加應(yīng)變εij計算單個粒子的本征應(yīng)變，由此實現(xiàn)等效夾雜物的替換，則含有粒子的固體成為了形式上的均質(zhì)體，可由本征應(yīng)變邊界積分方程（1）和（2）來描述.對于單個粒子時的情況，式（6）離散后的矩陣形式為

1.3局部Eshelby矩陣

如前所述，在進行等效夾雜物替換時，需要利用外加應(yīng)變來確定粒子的本征應(yīng)變.然而，對于多個粒子的情況，粒子之間存在相互作用，其作用的大小與粒子間的距離有關(guān).除了粒子本身的狀態(tài)，粒子上的外加應(yīng)變還受到周邊其他粒子的干擾作用，從而影響算法的迭代收斂性［15-17］.如圖1所示，對于當(dāng)前粒子I，定義虛線圓內(nèi)的粒子為近場粒子，其數(shù)量為NL，虛線圓外的粒子則為遠場粒子.對于當(dāng)前粒子來說，近場粒子屬于近距離粒子，其相互作用較為強烈，而遠場粒子對于當(dāng)前粒子的影響則相對較弱.

圖1　多粒子群組的定義Fig.1 Group definitions for multiple particles

對于全空間，如果只考慮近場粒子的影響而暫時忽略遠場粒子，結(jié)合本構(gòu)關(guān)系以及積分類型的轉(zhuǎn)換關(guān)系［17］，當(dāng)前粒子I的應(yīng)力可用應(yīng)力邊界積分方程（2）表達，即為

式中，Eijkl為彈性張量.由于全空間不存在外邊界，式（2）中的2個邊界積分在式（9）中并不出現(xiàn).將Eshelby張量的概念推廣，并利用式（5）和（9），則當(dāng)前粒子的拘束應(yīng)變以及Eshelby張量可分別表達為

式中，當(dāng)上標(biāo)I/=J時，Eshelby張量表示近場粒子中不同粒子的本征應(yīng)變與拘束應(yīng)變之間的關(guān)系.與式（6）的替換關(guān)系類似，根據(jù)Hooke定律，近場粒子的等效夾雜物替換所應(yīng)滿足的關(guān)系式為

結(jié)合式（10）和（11），式（12）離散后的矩陣形式為

由此即可利用近場粒子的外加應(yīng)變向量來計算其本征應(yīng)變向量.與文獻［8］中對裂紋的分析類似，本工作將式（13）中的矩陣稱為Eshelby矩陣，該矩陣可以看作Eshelby張量和等效夾雜物的概念在數(shù)值方面的一種拓廣.應(yīng)當(dāng)指出的是：式（9）～（12）是解析式，反映了全空間近場粒子之間應(yīng)力應(yīng)變的特定關(guān)系，是對近場粒子狀態(tài)的準(zhǔn)確描述；而式（13）則是離散后的數(shù)值形式，是對近場粒子的狀態(tài)和關(guān)系的近似描述.

對于具體的問題，如果指定了虛線圓（見圖1）的半徑，一般情況下固體中的每個粒子都會與一個獨特的近場粒子群相對應(yīng).為了計算的方便，將式（13）改寫為式中，SI是對Eshelby矩陣求逆和縮并得到的，且與當(dāng)前粒子的本征應(yīng)變向量相對應(yīng)；對應(yīng)于當(dāng)前粒子的近場粒子的外加應(yīng)變向量，與式（13）中相同；下標(biāo)I遍歷固體中的所有粒子，其總數(shù)為NI.

這樣，基于本工作提出的本征應(yīng)變邊界積分方程與局部Eshelby矩陣的算法（簡稱Eigen），固體中每個粒子本征應(yīng)變的計算均由3部分構(gòu)成：第1部分為外加載荷作用下產(chǎn)生的各個粒子上的應(yīng)力或外加應(yīng)變；第2部分為近場粒子群的作用，其相互影響較為強烈，利用式（14）進行計算；第3部分為遠場粒子群對當(dāng)前粒子的作用，其影響相對較弱，可直接利用積分方程中的域積分來計算.事實上，其他粒子對當(dāng)前粒子的作用大小和強弱是由無量綱距離來決定的，由此無量綱距離可作為粒子群中的近場和遠場粒子定義的依據(jù).顯然，近場粒子群的最小數(shù)量為1，如果NL=1，本工作算法就退化為文獻［15-17］中提出的算法，這時式（14）簡化為式（8）.

2　全空間邊界元子域法

本工作的主要目的是驗證所提出計算模型的正確性和方法的可行性，因此選擇邊界元子域法（boundary element method，BEM）作為對照的基準(zhǔn).文獻［18］中對全平面中的2個粒子進行了二維分析，但由于全空間沒有外邊界，其加載條件具有一定的特殊性，因此有必要對全空間邊界元子域法（見圖2）作簡要的介紹.

圖2　全空間邊界元子域法的示意圖Fig.2 Schematics of the sub-domain BEM in full space

如圖2（a）所示，考慮遠場載荷作用下全空間中的任一粒子?I，其邊界記作ΓI.現(xiàn)在保持界面狀態(tài)不變，將圖2（a）分解成內(nèi)域（圖2（b））和外域（圖2（c））兩種情況：圖2（b）中粒子的邊界上存在位移uI與面力τI，但無載荷作用；而圖2（c）中孔洞的邊界ΓI上，除了位移uI與面力—τI以外，還存在遠場載荷.注意內(nèi)域和外域之間滿足界面條件，即位移連續(xù)、面力平衡，因此二者面力的符號相反.由于全空間不存在外邊界，為了使邊界積分方程能夠表達遠場載荷的作用，須將圖2（c）進一步分解為圖2（d）和圖2（e），其中圖2（d）所示的虛擬界面（虛線）上的位移與面力-容易計算，圖2（e）所示的孔洞則適于全空間邊界積分方程的描述.

按照圖2所示的方法分解，對于含有NI個粒子的全空間，每個粒子都需要單獨進行描述，即需要NI個邊界積分方程，而分解后全空間中的所有孔洞（圖2（e）中只畫出了1個孔洞）則用另一個整體的邊界積分方程來描述.通過界面條件將這NI+1個方程聯(lián)系起來形成積分方程組，離散后求解，即所謂全空間邊界元子域法.邊界積分方程為

注意，式（15）中不含域積分.對于每個粒子，邊界積分方程離散后的矩陣形式為

式中，UI和TI分別為位移與面力基本解形成的系數(shù)矩陣.如果矩陣UI可逆，ΓI上的面力向量τI可表達為

對于全空間中的NI個孔洞，邊界積分方程離散后的矩陣形式為

式中，UIJ和TIJ分別為矩陣U和T中的子矩陣元素，其中下標(biāo)I和J分別表示源點p和場點q位于ΓI和ΓJ上.當(dāng)I=J時，子矩陣元素位于式（19）中“=”左邊矩陣的主對角元上，描述的是同一個孔洞.事實上，根據(jù)積分核的性質(zhì)，下列關(guān)系式成立：

式（21）成立的前提是假設(shè)粒子與基體具有相同的Poisson比，其中βI=EI/EM為粒子與基體的楊氏模量比.利用式（19）求得界面位移u后，再利用式（17）依次求出每個粒子界面上的面力τI.

3　算例

本工作的主要目地是檢驗計算模型，即Eshelby矩陣的正確性.假定粒子上的本征應(yīng)變均勻分布，且只對全空間含有少量橢球粒子的情況進行數(shù)值計算，則在計算中只選擇一個近場粒子群即可，其近場粒子群中粒子的數(shù)量與粒子的總數(shù)相同，即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個粒子的應(yīng)力分布

在遠場z向單位載荷作用下，全空間中的2個水平排列的扁橢球與其中心連線上的應(yīng)力分布如圖4所示.分別采用本征應(yīng)變邊界積分方程模型和全空間邊界元子域法，計算了2個扁橢球中心連線上的應(yīng)力，其中軟粒子的模量比為EI/EM=0.01，硬粒子的模量比為EI/EM=10.從計算結(jié)果（見圖4（b））可以看出，采用2種方法得到的結(jié)果十分吻合，說明了Eshelby矩陣的正確性與算法的可行性.另外，由圖4（b）可知，無論是軟粒子還是硬粒子，其應(yīng)力分量σ33在界面上都出現(xiàn)了跳躍.

圖4　2個水平排列的扁橢球粒子中心連線上的應(yīng)力分布Fig.4 Two horizontally placed oblate ellipsoidal particles and the stress distributions

圖5為全空間中水平排列的2個不同彈性模量的球形粒子和其中心連線上的應(yīng)力分布.圖5（a）中左側(cè)為軟粒子，右側(cè)為硬粒子.由圖5（b）可以看出，2種方法得到的結(jié)果吻合良好.

圖5　2個水平排列的不同彈性模量球形粒子中心連線上的應(yīng)力分布Fig.5 Two horizontally placed spheroidal particles with different modulus and the stress distributions

圖6為全空間中垂直排列的2個球形粒子和其中心連線上的應(yīng)力分布.圖6（b）所示的結(jié)果同樣說明2種方法的結(jié)果吻合良好.另外，由圖6（b）可以看出，在2個球形粒子垂直排列的情況下，軟粒子和硬粒子的應(yīng)力分量σ11都在界面上出現(xiàn)了跳躍.

圖6　2個垂直排列的球形粒子中心連線上的應(yīng)力分布Fig.6 Two vertically placed spheroidal particles and the stress distributions

3.23個粒子的應(yīng)力分布

在遠場z向單位載荷作用下，全空間中的3個水平排列的長橢球粒子如圖7（a）所示，分別采用本征應(yīng)變邊界積分方程模型和全空間邊界元子域法計算了右側(cè)2個長橢球粒子中心連線上的應(yīng)力，其中軟粒子的模量比為EI/EM=0.01，硬粒子的模量比為EI/EM=10.計算結(jié)果（見圖7（b））依然表明2種方法得到的結(jié)果十分吻合，說明了Eshelby矩陣的正確性與算法的可行性.同樣，由圖7（b）可知，無論是軟粒子還是硬粒子，其應(yīng)力分量σ33在界面上都出現(xiàn)了跳躍，但是2個長橢球粒子的狀態(tài)并不完全相同，另外，2個長橢球粒子內(nèi)部的應(yīng)力也是有區(qū)別的.

圖7　3個水平排列的長橢球粒子中右側(cè)2個橢球粒子中心連線上的應(yīng)力分布Fig.7 Three horizontally placed prolate ellipsoidal particles and the stress distributions

圖8為全空間中水平排列的3個球形粒子和右側(cè)2個粒子中心連線上的應(yīng)力分布.圖8（b）所示的結(jié)果也同樣說明2種方法得到的結(jié)果吻合良好.比較圖8（b）與圖7（b）的結(jié)果，可以看出粒子幾何形狀對應(yīng)力分布的影響.對于硬粒子，球形粒子內(nèi)部的應(yīng)力分量σ33高于長橢球粒子，而基體的應(yīng)力分量σ33有所降低；對于軟粒子，基體的應(yīng)力分量σ33有所上升.

圖8　3個水平排列的球形粒子中右側(cè)2個粒子中心連線上的應(yīng)力分布Fig.8 Three horizontally placed spheroidal particles and the stress distributions

3.35個粒子的應(yīng)力分布

在遠場z向單位載荷作用下，全空間中的5個在豎直平面內(nèi)排列的球形粒子如圖9（a）所示，分別采用本征應(yīng)變邊界積分方程模型和全空間邊界元子域法計算了中間粒子和右側(cè)粒子二者中心連線上的應(yīng)力，其中軟粒子的模量比為EI/EM=0.01，硬粒子的模量比為EI/EM=10.計算結(jié)果（見圖9（b））依然表明2種方法得到的結(jié)果十分吻合，說明了Eshelby矩陣的正確性與算法的可行性.同樣，由圖9（b）可知，無論是軟粒子還是硬粒子，其應(yīng)力分量σ33在界面上都出現(xiàn)了跳躍.另外，由于中間粒子和右側(cè)粒子的狀態(tài)有差別，圖9中2個粒子內(nèi)部應(yīng)力的差別比圖8中2個粒子差別更大.

圖9　5個球形粒子中的中間粒子與右側(cè)粒子中心連線上的應(yīng)力分布Fig.9 Horizontal computing line for five spheroidal particles and the stress distributions

圖10（b）為5個球形粒子中間粒子的和上方粒子二者中心連線上的應(yīng)力分布.圖10（b）的結(jié)果同樣說明2種方法的結(jié)果吻合良好.比較圖10與圖6，可以看出粒子排布對應(yīng)力分布的影響.

圖10　5個球形粒子中的中間粒子與上側(cè)粒子中心連線上的應(yīng)力分布Fig.10 Vertical computing line for five spheroidal particles and the stress distributions

眾所周知，含有粒子的介質(zhì)屬于異質(zhì)體問題，界面上的應(yīng)力分量有些是連續(xù)的，而有些會出現(xiàn)跳躍，例如軟粒子（EI/EM=0.01）的應(yīng)力分量σ33.采用邊界元子域法時，界面相當(dāng)于子域的邊界，因此界面應(yīng)力的計算需要分別在兩側(cè)的邊界上同時進行.應(yīng)當(dāng)指出的是，界面上的應(yīng)力分量出現(xiàn)跳躍恰恰反映了界面兩側(cè)材料具有不同彈性模量的影響.采用本征應(yīng)變邊界積分方程模型時，由于等效夾雜物的替換，計算區(qū)域成為表觀均質(zhì)體，因此界面應(yīng)力的計算可以直接在界面上進行.計算結(jié)果表明，對于應(yīng)力分量的數(shù)值出現(xiàn)跳躍的情況，按本工作中的算法得到的界面應(yīng)力分量的數(shù)值取兩側(cè)數(shù)值的平均值，如圖4～9中的σ33，圖6和10中的σ11，σ22等，這一結(jié)果與前期工作得到的結(jié)果［17］相同.

3.4計算效率

本工作提出的本征應(yīng)變邊界積分方程模型的計算效率與全空間邊界元子域法相比有較大的差別.2種算法的平均CPU時間如表1所示.由表1可以看出，2種算法的CPU時間相差達2個數(shù)量級，其中本征應(yīng)變邊界積分方程模型具有較高的計算效率.另外，由表1中2種算法CPU時間的比值可以看出，隨著粒子數(shù)量的增加，2種算法計算效率的差別將增大.這是因為全空間邊界元子域法的計算時間主要取決于總體系數(shù)矩陣的求解，其矩陣規(guī)模隨著粒子的數(shù)量而增大，因此CPU時間按幾何級數(shù)的方式增長；而在本征應(yīng)變邊界積分方程的算法中，當(dāng)指定了近場粒子群的大小之后，其CPU時間大體上按算術(shù)級數(shù)的方式增長.

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

4　結(jié)束語

針對具有密集分布粒子的粒子增強材料的大規(guī)模數(shù)值模擬問題，本工作利用全空間近場粒子群的積分方程離散形式，提出了局部Eshelby矩陣的概念，這一概念可以看作Eshelby張量和等效夾雜物的概念在數(shù)值方面的一種拓廣.通過將局部Eshelby矩陣引入到本征應(yīng)變邊界積分方程計算模型中，改進了原有的算法，解決了近場粒子間的相互作用問題，并且該問題將影響算法的迭代收斂性.以邊界元子域法為參照，利用本征應(yīng)變邊界積分方程計算模型對無限域中的若干橢球粒子和球形粒子進行了三維應(yīng)力分析，驗證了模型的正確性、方法的可行性.結(jié)果表明，Eigen模型具有較高的計算效率，具有對粒子增強材料大規(guī)模數(shù)值分析的能力.

作為初步工作，粒子中的本征應(yīng)變被假定為常數(shù)，這是本工作中計算模型的局限所在.進一步的工作需要選取適當(dāng)?shù)牟逯捣绞?，取消本征?yīng)變的常數(shù)假定，從而將計算模型推廣到一般幾何形狀粒子的分析中.

［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，F(xiàn)ang 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.本文彩色版可登陸本刊網(wǎng)站查詢：http：//www.journal.shu.edu.cn

Eigenstrain boundary integral equation with local Eshelby matrix for ellipsoidal particles

MA Hang1，F(xiàn)ANG 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

國家自然科學(xué)基金資助項目（11272195，11332005）

馬杭（1951—），男，教授，博士生導(dǎo)師，研究方向為計算固體力學(xué).E-mail：hangma@staff.shu.edu.cn