# FM-DBEM Simulation of 3D Microvoid and Microcrack Graphite Models

Authors: Houdi Lu; Hongtao Wang; Haitao Wang; Lie Jin; Xinxin Wu; Yu Zhou
Journal: Science and Technology of Nuclear Installations
 (2017)

---

## Abstract

The graphite is porous medium, and the geometry and size distribution of its structural deficiencies such as microcracks and microvoids at different oxidation degrees have a great influence on the overall performance. In this paper, we adopt the FM-DBEM to study 3D models which contain spheroidal microvoids and circular microcracks. The accuracy of this method is tested by a comparison to the theoretical solution to the problem of 2D microcrack and microvoid interaction problem. Two simulations are conducted: the simulation of graphite model containing a large number of randomly distributed microcracks and microvoids and the simulation of graphite model containing microcracks and growing microvoids. The simulations investigate the effective moduli versus the two microstructures’ density and the effect of microvoid’s growth on the SIF of microcrack.


---

## Body

## 1. Introduction

Nuclear graphite is widely used as fuel blocks and reactor internals in high-temperature gas-cooled reactors (HTR) because of its high-temperature resistance and neutron moderation [1, 2]. As the structural material of HTR, the mechanical property and structural integrity of graphite structure have a great impact on HTR’s safe operation. The oxide gas and possible accidents will lead to the oxidation of graphite and influence its mesostructure and mechanical parameter. As a kind of porous media, there exist structural deficiencies such as microcracks and microvoids in graphite. During the oxidation process, the microvoids in the graphite will grow, and the microcrack tip often fractures first during the fracture process. Therefore, it is of great value to study the interaction between microcracks and microvoids in the graphite.Many researchers have developed continuum damage micromechanical models to study graphite as a kind of brittle material. Kachanov [3, 4] and Ju [5] discussed the loss of moduli in terms of evolution of the microstructure, but their research did not involve material parameters. Hu and Chandra [6] and Huang et al. [7] adopted a unit cell approach and obtained the stress intensity factors by modeling a finite number of defects. Hu et al. [8] developed an integral equation to model interactions between voids and microcracks of 2D elastic fracture mechanics. The local behavior of such defective material has been studied by many other researchers [9–12] using theoretical approaches, but the theoretical solutions are obtained for simple cases, such as a void among several regularly spaced cracks or a crack among several regularly spaced voids embedded in an infinite plate subjected to remote loading. Another characteristic of graphite is its anisotropy. The anisotropy of the coefficient of thermal expansion in crystallites and the generation of internal stresses beyond the elastic limit of graphite crystals cause the cracking in graphite [13, 14]. Yan et al. [15] developed a 2D micromechanistic model which considered intergranular cracks and intragranular cracks to predict the response of graphite under temperature changes and irradiation, and this model was appropriate for generalized plane strain condition. At present, much less literature has been focused on 3D theoretical solutions and the 3D simulation of the interaction between microvoids and microcracks.The finite element method (FEM) is also widely used to solve the microvoid and microcrack problems in some practical structures. In order to ensure the accuracy of simulation, a large number of elements are used, which largely reduces the computational efficiency. Even under such a condition, the accuracy is not that satisfactory because the assumed displacement functions do not satisfy the dominance conditions at the free boundary. In order to solve such problems, some modified methods [16, 17] are applied, but these methods need to model the vicinity of microvoids and microcracks with a mass of elements; besides, special elements are also defined to tackle the problem [18, 19].The boundary element method (BEM) is an effective method in solving microvoid and microcrack problems, since only the boundary of the analyzed domain needs to be discretized [20]. The advantages of BEM are that it can not only reduce the initial preparation data and the degrees of freedom, but also simplify the microcrack meshing process; thus, much fewer elements are needed than in FEM. When dealing with fracture problems, the dual boundary element method (DBEM) is applied to overcome the degeneration difficulty when the two surfaces of the same crack coincide [21]. In the research on the DBEM, Portela et al. [22, 23] first used this method for the analysis of 2D crack problem, and then this method was extended to 3D crack analysis by Mi and Aliabadi [24], Cisilino and Aliabadi [25], and Wilde and Aliabadi [26]. Meanwhile, many researchers [27–30] have adopted the fast multipole method (FMM) (Rokhlin [31]) in their BEM because of FMM’s ability to reduce memory requirement and calculation scale. Another fast algorithm widely adopted for crack problems is hierarchical matrices. The hierarchical matrices method is an algebraic method, and its approximation is based only on the knowledge of individual matrix entries. Research on hierarchical matrices DBEM for 3D general elastic crack problems has been reported by Benedetti et al. [32–34]. The hierarchical matrices method has better computational speed than FMM, but FMM is more capable for large scale problem than the former because FMM’s cost of storing is smaller. So far, 3D microcrack problem has been studied by several researchers. Lai and Rodin [35] presented a fast BEM for analyzing 3D linear elastic solids containing many cracks. Nishimura et al. [36] discussed a three-dimensional fast multipole boundary integral equation method for crack problems for Laplace’s equation. Wang and Yao [37] adopted fast multipole-dual boundary element method (FM-DBEM) to analyze the 3D crack problems. Lu [38] introduced a modeling method based on Python and ABAQUS/CAE to get microvoid and microcrack models. However, the interaction of microvoid and microcrack still needs intensive research.In this paper, we adopt the FM-DBEM to explore 3D graphite models which contain spheroidal microvoids and circular microcracks and study their interaction. The 3D models’ effective moduli are evaluated in several conditions. A FM-DBEM solver written by C++ code is used and its accuracy is tested by a comparison to the theoretical solution to the 2D problem of the interaction between microcrack and microvoid. This method shows high accuracy and efficiency, and it can be applied to simulate complex structures and components containing microvoids and microcracks.


## 2. Models and Methods

### 2.1. Dual Boundary Element Method for 3D Microcrack and Microvoid Interaction

Figure1 shows the model of an elastic structure containing edge, embedded voids, and cracks. V and S0 represent the domain and outer boundary of the solid. S1 represent the inner boundary of the voids. SC+ and SC- represent the two coinciding surfaces of the cracks.Figure 1
An elastic solid containing edge and embedded cracks.The dual boundary integral equations for 3D elastic fracture mechanics without body force are shown as(1)cijxujx+∫S0Tij∗x,yujy-Uij∗x,ytjydSy=-∫SC+Tij∗x,yΔujydSyx∈S0,S1(2)12ui+x+ui-x+∫S0Tij∗x,yujy-Uij∗x,ytjydSy=-∫SC+Tij∗x,yΔujydSy,x∈SC-(3)12ti+x-ti-x+∫S0Sij∗x,yujy-Dij∗x,ytjydSy=-∫SC+Sij∗x,yΔujydSy,x∈SC+,where x and y represent the source and field point; ui and ti are the boundary displacement and traction vectors; cij(x) is a free term related to the shape of boundary at point x; Δui=ui+-ui- is the relative displacement between SC+ and SC-. The kernel functions Uij∗(x,y) and Tij∗(x,y) are shown as(4)Uij∗x,y=116πG1-νr3-4νδij+rirjr2Tij∗x,y=18π1-νr1-2νnirj-njrir-∂r∂n1-2νδij+3rirjr2,where r represents the distance between x and y; G and ν represent the shear modulus and Poisson’s ratio. Dij∗x,y and Sij∗x,y are obtained from Uij∗(x,y) and Tij∗(x,y):(5)Dij∗x,y=nkxλδki∂Ulj∗x,y∂xl+G∂Ukj∗x,y∂xi+∂Uij∗x,y∂xk=18π1-vr21-2vδijrknkxr+rinjx-rjnixr+3rirjrkr3nkxSij∗x,y=nkxλδki∂Tlj∗x,y∂xl+G∂Tkj∗x,y∂xi+∂Tij∗x,y∂xk=-G4π1-vr33∂r∂n5rirjrknkx-1-2vrjnix+1-4vnjynix-3vrirjnkynkx+rkrjniynkx+δijrknkx∂r∂n+rinjx∂r∂n-1-2vniynjx+nkynkxδij+3rirknjynkx,where nkx and nky represent the outward normal vectors at points x and y; λ is the Lame constant.Collocation and eight-node quadratic are used to discretize the dual boundary integral equations. The concept of finite part integral is used to deal with the strongly singular and hypersingular integrals, and it requires that the traction and displacement derivatives should be Hölder continuous. Besides, discontinuous elements are adopted for crack modeling and edge-discontinuous elements on surfaces are adopted for approaching the corner or intersecting the crack surface. Continuous elements are adopted on all other surfaces.The two-point crack opening displacement formulas are used to calculate the stress intensity factors. The relative displacement of two surfaces of the crack at the collocation pointsΔu is obtained by the DBEM analysis. The SIFs at geometry point Q are given:(6)KIQ=4KIP1-KIP23KIP=G21-vcos⁡φπ2lPQΔunPKIIQ=4KIIP1-KIIP23KIIP=G21-vcos⁡φπ2lPQΔut2PKIIIQ=4KIIIP1-KIIIP23KIIIP=G2cos⁡φπ2lPQΔut1P,where lPQ represents the distance between P and Q; Δun and Δut are the normal and tangent components of Δu at P under the local coordinate system defined at Q (see Figure 2).Figure 2
Generation of quarter point element.


### 2.2. Fast Multipole-Dual Boundary Element Method

Considering a reference pointO close to y, if a threshold of Oy→<Ox→ is satisfied, the integrals in (1) and (2) can be expanded around O as(7)∫S0Uij∗x,ytjydSy=18πG∑n=0∞∑m=-nnFij,n,mS¯Ox→Mj,n,mU1O+Gi,n,mS¯Ox→Mn,mU2O∫S0Tij∗x,yujydSy=18πG∑n=0∞∑m=-nnFij,n,mS¯Ox→Mj,n,mT1O+Gi,n,mS¯Ox→Mn,mT2O∫Sc+Tij∗x,yΔujydSy=18πG∑n=0∞∑m=-nnFij,n,mS¯Ox→Mj,n,mT1,ScO+Gi,n,mS¯Ox→Mn,mT2,ScO,where(8)Fij,n,mSOx→=λ+3Gλ+2GδijSn,mOx→-λ+Gλ+2GOx→j∂∂xiSn,mOx→Gi,n,mSOx→=λ+Gλ+2G∂∂xiSn,mOx→Sn,mOy→=n-m!Pnmcos⁡θeimϕ1rn+1.Mj,n,mU1O,Mn,mU2O,Mj,n,mT1O,Mn,mT2O,Mj,n,mT1,ScO, and Mn,mT2,ScO are called multipole moments presented by Wang and Yao [37]; Sn,m is called solid spherical harmonic function [39].Equation (3) can be expanded as(9)∫S0Dij∗x,ytjydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mS¯Ox→∂xl+G∂Fkj,n,mS¯Ox→∂xi+∂Fij,n,mS¯Ox→∂xkMj,n,mU1O+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mS¯Ox→∂xl+G∂Gk,n,mS¯Ox→∂xi+∂Gi,n,mS¯Ox→∂xkMn,mU2O∫S0Sij∗x,yujydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mS¯Ox→∂xl+G∂Fkj,n,mS¯Ox→∂xi+∂Fij,n,mS¯Ox→∂xkMj,n,mT1O+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mS¯Ox→∂xl+G∂Gk,n,mS¯Ox→∂xi+∂Gi,n,mS¯Ox→∂xkMn,mT2O∫Sc+Sij∗x,yΔujydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mS¯Ox→∂xl+G∂Fkj,n,mS¯Ox→∂xi+∂Fij,n,mS¯Ox→∂xkMj,n,mT1,ScO+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mS¯Ox→∂xl+G∂Gk,n,mS¯Ox→∂xi+∂Gi,n,mS¯Ox→∂xkMn,mT2,ScO.Considering another reference point O′ close to x, if a threshold O′x→<O′y→ is satisfied, the integrals in (1), (2), and (3) can be expanded around O′ in a local expansion form:(10)∫S0Uij∗x,ytjydSy=18πG∑n=0∞∑m=-nnFij,n,mRO′x→Lj,n,mU1O′+Gi,n,mRO′x→Ln,mU2O′∫S0Tij∗x,yujydSy=18πG∑n=0∞∑m=-nnFij,n,mRO′x→Lj,n,mT1O′+Gi,n,mRO′x→Ln,mT2O′∫Sc+Tij∗x,yΔujydSy=18πG∑n=0∞∑m=-nnFij,n,mRO′x→Lj,n,mT1,ScO′+Gi,n,mRO′x→Ln,mT2,ScO′∫S0Dij∗x,ytjydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mRO′x→∂xl+G∂Fkj,n,mRO′x→∂xi+∂Fij,n,mRO′x→∂xkLj,n,mU1O′+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mRO′x→∂xl+G∂Gk,n,mRO′x→∂xi+∂Gi,n,mRO′x→∂xkLn,mU2O′∫S0Sij∗x,yujydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mRO′x→∂xl+G∂Fkj,n,mRO′x→∂xi+∂Fij,n,mRO′x→∂xkLj,n,mT1O′+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mRO′x→∂xl+G∂Gk,n,mRO′x→∂xi+∂Gi,n,mRO′x→∂xkLn,mT2O′∫Sc+Sij∗x,yΔujydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mRO′x→∂xl+G∂Fkj,n,mRO′x→∂xi+∂Fij,n,mRO′x→∂xkLj,n,mT1,ScO′+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mRO′x→∂xl+G∂Gk,n,mRO′x→∂xi+∂Gi,n,mRO′x→∂xkLn,mT2,ScO′,where(11)Fij,n,mRO′x→=λ+3Gλ+2GδijRn,mO′x→-λ+Gλ+2GO′x→j∂∂xiRn,mO′x→Gi,n,mRO′x→=λ+Gλ+2G∂∂xiRn,mO′x→Rn,mOy→=1n+m!Pnmcos⁡θeimϕrn.Lj,n,mU1O′,Ln,mU2O′,Lj,n,mT1O′,Ln,mT2O′,Lj,n,mT1,ScO′, and Ln,mT2,ScO′ are called local moments. Instead of being evaluated directly, the local moments are derived by a linear mapping acting on the multipole moments in the FM framework. Rn,m is another solid spherical harmonic function.Based on the geometry information of the boundary element, an adaptive tree is constructed. By introducing an exponential expansion, three operators, namely, the multipole to local translation, exponential to exponential translation, and exponential to local translation, are carried out recursively throughout the tree in order to obtain both the multipole moments and local moments for the tree nodes of various levels. Yoshida [39] presented the formulas of the translation operators for 3D elastostatics problems.In the iterative solution of DBEM, GMRES is adopted. A left preconditioner matrix with block diagonal forms presented by Nishimura et al. [36] is used in this paper. Each diagonal block corresponds to one tree leaf, and entries of the block are evaluated by the collocation nodes contained in the leaf directly. In this paper, multipole, local, and exponential expansions are all of eighteenth order, and the 2-norm of the error in the FMM approximation presented by Greengard and Rokhlin [40] is less than 10−6. When the multipole and local expansion are of twenty-fourth order and the exponential expansion is of twenty-seventh order, the value difference of SIF [37] between this expansion order and the eighteenth order is within 0.5% of analytical value, showing the adequacy of the eighteenth order.


### 2.3. 3D Microvoid and Microcrack Graphite Model Simulated Using FM-DBEM

Considering a representative elastic solid cube with spheroidal microvoids and circular microcracks, see Figure3. Based on previous circular and elliptical microcrack studies in BEM [41], this paper inherits most of the symbols and characters from previous research.Figure 3
An isotropic cube containing randomly located spheroidal microvoids and circular microcracks.The formula for the constitutive relation of elastic problem without body force is(12)σij=Cijklεkl,where εkl is the elastic strain, Cijkl is the elastic coefficients, and the elastic strain field εkl is defined by(13)εkl=12ui,j+uj,i,where ui is the elastic displacement.The boundary conditions are given by mechanical boundary conditions:(14)ti=σijnj=ti¯ui=ui¯,where ti is the surface traction and ni is the component of the unit outward normal vector of the surface.The side length of the cubic isL. For the spheroidal microvoid in this solid, if its radius is Rv, then its volume is defined as Vv=4/3πRv3. If the number of spheroidal microvoids is Nv, the spheroidal microvoid density parameter in the solid is defined by ω1=4πRv3Nv/3L3. Similarly, the radius of the circular microcrack is Rc and the number of elliptical microcracks is Nc. Following Budiansky and O’connell [42] and Kachanov [43], the microcrack density parameter is defined as ω2=Rc3Nc/L3. Young’s modulus E0 and Poisson’s ratio ν0 of the matrix are 9.8 Gpa and 0.14, respectively.The weakened solid model shows similar behaviors as an isotropic elastic medium when the amount of microvoids and microcracks is large enough. The procedure can be used to obtain the effective Young’s modulusEeff. The average strain and stress in the solid are defined by(15)εij¯=1V∫VεijdVσij¯=1V∫VσijdV.When different boundary conditions are applied, the material constants and effective moduli can be obtained.Applying the boundary condition as in(16)Surfacex1+:u-1=u0,t-2=t-3=0u0≠0Surfacex1-:u-1=0,t-2=t-3=0Surfacex2+:t-1=t-2=t-3=0Surfacex2-:t-1=0,u-2=0,t-3=0Surfacex3+:t-1=t-2=t-3=0Surfacex3-:t-1=t-2=0,u-3=0then effective Young’s modulus in direction X1 can be obtained by(17)E1eff=σ11¯ε11¯Applying the boundary condition as in(18)Surfacex1+:u-1=u0,t-2=t-3=0u0≠0Surfacex1-:u-1=0,t-2=t-3=0Surfacex2+:t-1=0,u-2=u0,t-3=0Surfacex2-:t-1=0,u-2=0,t-3=0Surfacex3+:t-1=t-2=0,u-3=u0Surfacex3-:t-1=t-2=0,u-3=0.we can obtain(19)Keff=σ11¯+σ22¯+σ33¯3ε11¯+ε22¯+ε33¯.


### 2.4. A Computational Test of Microcrack and Microvoid Interaction under Plane Strain State

The current literature on the theoretical analysis of the microcrack-microvoid interaction model is based on the 2D model [8, 11, 44, 45], and there is less research on the theoretical solution of the 3D microcrack-microvoid interaction. Therefore, the test algorithm in this section adopts a model with a large thickness to approximate the 2D plane strain state.The 2D crack and hole model and 3D meshing schematic diagram of BEM are shown in Figures4 and 6. In Figure 4, the centers of the single crack and the single hole are located on the same horizontal axis, and there is balanced vertical tensile stress in the up and down direction. The radius of the 2D circular hole is R, the crack length is 2L, R=L, and the vertical distance from the crack tip to the hole’s margin is D. In the simulated model, a large external matrix is established. The thickness of the circular hole and the crack is much larger than the sizes of the hole and the crack. The midsection in the model can be equaled to the approximate plane strain state.Figure 4
The schematic diagram of a crack aligned along the center of the hole.
(a)
The schematic diagram of the 2D model
(b)
The geometry and mesh of 3D modelFigure5 shows the comparison between the simulation result of this test model and the 2D theoretical approximate solution. The crack tips’ SIFs from FM-DBEM conform to the theoretical approximate solution, especially when distance D is small. The crack tip to the left of the hole enlarges fast with the decrease of distance D. The tip’s stress intensity factor K1 shall become 2-3 times larger than when there is no hole and but only the single crack.Figure 5
The comparison between the result from the 3D model in Figure4 and the 2D theoretical solution.Figure 6
The schematic diagram of the inclined crack and hole model.
(a)
The schematic diagram of the 2D model
(b)
The geometry and mesh of 3D modelMoreover, if we spin the crack in the above model with anα angle, a 2D inclined crack and hole model is formed. In Figure 6, the hole’s radius is R, the crack’s length maintains 2L, and their distance is D1=2.2R, which is equivalent to D=1.2R in the above model.Figure7 shows the comparison between the simulation result and the 2D theoretical approximate solution. In this inclined crack and hole model, because of the existence of crack’s inclined angle, its tip’s SIFs K1 and K2 shall change. When we keep the distance between centers of the crack and the hole to D1 and change the crack’s inclined angle α from 0° to 90°, the K1 value shall decrease to 0 with the increase in the angle and there is no apparent difference between the two K2 values. Besides, they both increase first and then decrease. When the loading direction is vertical to the crack’s plane, the K1 value is the largest and the crack is most likely to start cracking. When the loading direction is parallel to the crack’s plane, there is basically no stress concentration.Figure 7
The comparison between the result from the 3D model in Figure6 and the 2D theoretical solution.


## 2.1. Dual Boundary Element Method for 3D Microcrack and Microvoid Interaction

Figure1 shows the model of an elastic structure containing edge, embedded voids, and cracks. V and S0 represent the domain and outer boundary of the solid. S1 represent the inner boundary of the voids. SC+ and SC- represent the two coinciding surfaces of the cracks.Figure 1
An elastic solid containing edge and embedded cracks.The dual boundary integral equations for 3D elastic fracture mechanics without body force are shown as(1)cijxujx+∫S0Tij∗x,yujy-Uij∗x,ytjydSy=-∫SC+Tij∗x,yΔujydSyx∈S0,S1(2)12ui+x+ui-x+∫S0Tij∗x,yujy-Uij∗x,ytjydSy=-∫SC+Tij∗x,yΔujydSy,x∈SC-(3)12ti+x-ti-x+∫S0Sij∗x,yujy-Dij∗x,ytjydSy=-∫SC+Sij∗x,yΔujydSy,x∈SC+,where x and y represent the source and field point; ui and ti are the boundary displacement and traction vectors; cij(x) is a free term related to the shape of boundary at point x; Δui=ui+-ui- is the relative displacement between SC+ and SC-. The kernel functions Uij∗(x,y) and Tij∗(x,y) are shown as(4)Uij∗x,y=116πG1-νr3-4νδij+rirjr2Tij∗x,y=18π1-νr1-2νnirj-njrir-∂r∂n1-2νδij+3rirjr2,where r represents the distance between x and y; G and ν represent the shear modulus and Poisson’s ratio. Dij∗x,y and Sij∗x,y are obtained from Uij∗(x,y) and Tij∗(x,y):(5)Dij∗x,y=nkxλδki∂Ulj∗x,y∂xl+G∂Ukj∗x,y∂xi+∂Uij∗x,y∂xk=18π1-vr21-2vδijrknkxr+rinjx-rjnixr+3rirjrkr3nkxSij∗x,y=nkxλδki∂Tlj∗x,y∂xl+G∂Tkj∗x,y∂xi+∂Tij∗x,y∂xk=-G4π1-vr33∂r∂n5rirjrknkx-1-2vrjnix+1-4vnjynix-3vrirjnkynkx+rkrjniynkx+δijrknkx∂r∂n+rinjx∂r∂n-1-2vniynjx+nkynkxδij+3rirknjynkx,where nkx and nky represent the outward normal vectors at points x and y; λ is the Lame constant.Collocation and eight-node quadratic are used to discretize the dual boundary integral equations. The concept of finite part integral is used to deal with the strongly singular and hypersingular integrals, and it requires that the traction and displacement derivatives should be Hölder continuous. Besides, discontinuous elements are adopted for crack modeling and edge-discontinuous elements on surfaces are adopted for approaching the corner or intersecting the crack surface. Continuous elements are adopted on all other surfaces.The two-point crack opening displacement formulas are used to calculate the stress intensity factors. The relative displacement of two surfaces of the crack at the collocation pointsΔu is obtained by the DBEM analysis. The SIFs at geometry point Q are given:(6)KIQ=4KIP1-KIP23KIP=G21-vcos⁡φπ2lPQΔunPKIIQ=4KIIP1-KIIP23KIIP=G21-vcos⁡φπ2lPQΔut2PKIIIQ=4KIIIP1-KIIIP23KIIIP=G2cos⁡φπ2lPQΔut1P,where lPQ represents the distance between P and Q; Δun and Δut are the normal and tangent components of Δu at P under the local coordinate system defined at Q (see Figure 2).Figure 2
Generation of quarter point element.


## 2.2. Fast Multipole-Dual Boundary Element Method

Considering a reference pointO close to y, if a threshold of Oy→<Ox→ is satisfied, the integrals in (1) and (2) can be expanded around O as(7)∫S0Uij∗x,ytjydSy=18πG∑n=0∞∑m=-nnFij,n,mS¯Ox→Mj,n,mU1O+Gi,n,mS¯Ox→Mn,mU2O∫S0Tij∗x,yujydSy=18πG∑n=0∞∑m=-nnFij,n,mS¯Ox→Mj,n,mT1O+Gi,n,mS¯Ox→Mn,mT2O∫Sc+Tij∗x,yΔujydSy=18πG∑n=0∞∑m=-nnFij,n,mS¯Ox→Mj,n,mT1,ScO+Gi,n,mS¯Ox→Mn,mT2,ScO,where(8)Fij,n,mSOx→=λ+3Gλ+2GδijSn,mOx→-λ+Gλ+2GOx→j∂∂xiSn,mOx→Gi,n,mSOx→=λ+Gλ+2G∂∂xiSn,mOx→Sn,mOy→=n-m!Pnmcos⁡θeimϕ1rn+1.Mj,n,mU1O,Mn,mU2O,Mj,n,mT1O,Mn,mT2O,Mj,n,mT1,ScO, and Mn,mT2,ScO are called multipole moments presented by Wang and Yao [37]; Sn,m is called solid spherical harmonic function [39].Equation (3) can be expanded as(9)∫S0Dij∗x,ytjydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mS¯Ox→∂xl+G∂Fkj,n,mS¯Ox→∂xi+∂Fij,n,mS¯Ox→∂xkMj,n,mU1O+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mS¯Ox→∂xl+G∂Gk,n,mS¯Ox→∂xi+∂Gi,n,mS¯Ox→∂xkMn,mU2O∫S0Sij∗x,yujydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mS¯Ox→∂xl+G∂Fkj,n,mS¯Ox→∂xi+∂Fij,n,mS¯Ox→∂xkMj,n,mT1O+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mS¯Ox→∂xl+G∂Gk,n,mS¯Ox→∂xi+∂Gi,n,mS¯Ox→∂xkMn,mT2O∫Sc+Sij∗x,yΔujydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mS¯Ox→∂xl+G∂Fkj,n,mS¯Ox→∂xi+∂Fij,n,mS¯Ox→∂xkMj,n,mT1,ScO+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mS¯Ox→∂xl+G∂Gk,n,mS¯Ox→∂xi+∂Gi,n,mS¯Ox→∂xkMn,mT2,ScO.Considering another reference point O′ close to x, if a threshold O′x→<O′y→ is satisfied, the integrals in (1), (2), and (3) can be expanded around O′ in a local expansion form:(10)∫S0Uij∗x,ytjydSy=18πG∑n=0∞∑m=-nnFij,n,mRO′x→Lj,n,mU1O′+Gi,n,mRO′x→Ln,mU2O′∫S0Tij∗x,yujydSy=18πG∑n=0∞∑m=-nnFij,n,mRO′x→Lj,n,mT1O′+Gi,n,mRO′x→Ln,mT2O′∫Sc+Tij∗x,yΔujydSy=18πG∑n=0∞∑m=-nnFij,n,mRO′x→Lj,n,mT1,ScO′+Gi,n,mRO′x→Ln,mT2,ScO′∫S0Dij∗x,ytjydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mRO′x→∂xl+G∂Fkj,n,mRO′x→∂xi+∂Fij,n,mRO′x→∂xkLj,n,mU1O′+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mRO′x→∂xl+G∂Gk,n,mRO′x→∂xi+∂Gi,n,mRO′x→∂xkLn,mU2O′∫S0Sij∗x,yujydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mRO′x→∂xl+G∂Fkj,n,mRO′x→∂xi+∂Fij,n,mRO′x→∂xkLj,n,mT1O′+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mRO′x→∂xl+G∂Gk,n,mRO′x→∂xi+∂Gi,n,mRO′x→∂xkLn,mT2O′∫Sc+Sij∗x,yΔujydSy=18πG∑n=0∞∑m=-nnnkλδki∂Flj,n,mRO′x→∂xl+G∂Fkj,n,mRO′x→∂xi+∂Fij,n,mRO′x→∂xkLj,n,mT1,ScO′+18πG∑n=0∞∑m=-nnnkλδki∂Gl,n,mRO′x→∂xl+G∂Gk,n,mRO′x→∂xi+∂Gi,n,mRO′x→∂xkLn,mT2,ScO′,where(11)Fij,n,mRO′x→=λ+3Gλ+2GδijRn,mO′x→-λ+Gλ+2GO′x→j∂∂xiRn,mO′x→Gi,n,mRO′x→=λ+Gλ+2G∂∂xiRn,mO′x→Rn,mOy→=1n+m!Pnmcos⁡θeimϕrn.Lj,n,mU1O′,Ln,mU2O′,Lj,n,mT1O′,Ln,mT2O′,Lj,n,mT1,ScO′, and Ln,mT2,ScO′ are called local moments. Instead of being evaluated directly, the local moments are derived by a linear mapping acting on the multipole moments in the FM framework. Rn,m is another solid spherical harmonic function.Based on the geometry information of the boundary element, an adaptive tree is constructed. By introducing an exponential expansion, three operators, namely, the multipole to local translation, exponential to exponential translation, and exponential to local translation, are carried out recursively throughout the tree in order to obtain both the multipole moments and local moments for the tree nodes of various levels. Yoshida [39] presented the formulas of the translation operators for 3D elastostatics problems.In the iterative solution of DBEM, GMRES is adopted. A left preconditioner matrix with block diagonal forms presented by Nishimura et al. [36] is used in this paper. Each diagonal block corresponds to one tree leaf, and entries of the block are evaluated by the collocation nodes contained in the leaf directly. In this paper, multipole, local, and exponential expansions are all of eighteenth order, and the 2-norm of the error in the FMM approximation presented by Greengard and Rokhlin [40] is less than 10−6. When the multipole and local expansion are of twenty-fourth order and the exponential expansion is of twenty-seventh order, the value difference of SIF [37] between this expansion order and the eighteenth order is within 0.5% of analytical value, showing the adequacy of the eighteenth order.


## 2.3. 3D Microvoid and Microcrack Graphite Model Simulated Using FM-DBEM

Considering a representative elastic solid cube with spheroidal microvoids and circular microcracks, see Figure3. Based on previous circular and elliptical microcrack studies in BEM [41], this paper inherits most of the symbols and characters from previous research.Figure 3
An isotropic cube containing randomly located spheroidal microvoids and circular microcracks.The formula for the constitutive relation of elastic problem without body force is(12)σij=Cijklεkl,where εkl is the elastic strain, Cijkl is the elastic coefficients, and the elastic strain field εkl is defined by(13)εkl=12ui,j+uj,i,where ui is the elastic displacement.The boundary conditions are given by mechanical boundary conditions:(14)ti=σijnj=ti¯ui=ui¯,where ti is the surface traction and ni is the component of the unit outward normal vector of the surface.The side length of the cubic isL. For the spheroidal microvoid in this solid, if its radius is Rv, then its volume is defined as Vv=4/3πRv3. If the number of spheroidal microvoids is Nv, the spheroidal microvoid density parameter in the solid is defined by ω1=4πRv3Nv/3L3. Similarly, the radius of the circular microcrack is Rc and the number of elliptical microcracks is Nc. Following Budiansky and O’connell [42] and Kachanov [43], the microcrack density parameter is defined as ω2=Rc3Nc/L3. Young’s modulus E0 and Poisson’s ratio ν0 of the matrix are 9.8 Gpa and 0.14, respectively.The weakened solid model shows similar behaviors as an isotropic elastic medium when the amount of microvoids and microcracks is large enough. The procedure can be used to obtain the effective Young’s modulusEeff. The average strain and stress in the solid are defined by(15)εij¯=1V∫VεijdVσij¯=1V∫VσijdV.When different boundary conditions are applied, the material constants and effective moduli can be obtained.Applying the boundary condition as in(16)Surfacex1+:u-1=u0,t-2=t-3=0u0≠0Surfacex1-:u-1=0,t-2=t-3=0Surfacex2+:t-1=t-2=t-3=0Surfacex2-:t-1=0,u-2=0,t-3=0Surfacex3+:t-1=t-2=t-3=0Surfacex3-:t-1=t-2=0,u-3=0then effective Young’s modulus in direction X1 can be obtained by(17)E1eff=σ11¯ε11¯Applying the boundary condition as in(18)Surfacex1+:u-1=u0,t-2=t-3=0u0≠0Surfacex1-:u-1=0,t-2=t-3=0Surfacex2+:t-1=0,u-2=u0,t-3=0Surfacex2-:t-1=0,u-2=0,t-3=0Surfacex3+:t-1=t-2=0,u-3=u0Surfacex3-:t-1=t-2=0,u-3=0.we can obtain(19)Keff=σ11¯+σ22¯+σ33¯3ε11¯+ε22¯+ε33¯.


## 2.4. A Computational Test of Microcrack and Microvoid Interaction under Plane Strain State

The current literature on the theoretical analysis of the microcrack-microvoid interaction model is based on the 2D model [8, 11, 44, 45], and there is less research on the theoretical solution of the 3D microcrack-microvoid interaction. Therefore, the test algorithm in this section adopts a model with a large thickness to approximate the 2D plane strain state.The 2D crack and hole model and 3D meshing schematic diagram of BEM are shown in Figures4 and 6. In Figure 4, the centers of the single crack and the single hole are located on the same horizontal axis, and there is balanced vertical tensile stress in the up and down direction. The radius of the 2D circular hole is R, the crack length is 2L, R=L, and the vertical distance from the crack tip to the hole’s margin is D. In the simulated model, a large external matrix is established. The thickness of the circular hole and the crack is much larger than the sizes of the hole and the crack. The midsection in the model can be equaled to the approximate plane strain state.Figure 4
The schematic diagram of a crack aligned along the center of the hole.
(a)
The schematic diagram of the 2D model
(b)
The geometry and mesh of 3D modelFigure5 shows the comparison between the simulation result of this test model and the 2D theoretical approximate solution. The crack tips’ SIFs from FM-DBEM conform to the theoretical approximate solution, especially when distance D is small. The crack tip to the left of the hole enlarges fast with the decrease of distance D. The tip’s stress intensity factor K1 shall become 2-3 times larger than when there is no hole and but only the single crack.Figure 5
The comparison between the result from the 3D model in Figure4 and the 2D theoretical solution.Figure 6
The schematic diagram of the inclined crack and hole model.
(a)
The schematic diagram of the 2D model
(b)
The geometry and mesh of 3D modelMoreover, if we spin the crack in the above model with anα angle, a 2D inclined crack and hole model is formed. In Figure 6, the hole’s radius is R, the crack’s length maintains 2L, and their distance is D1=2.2R, which is equivalent to D=1.2R in the above model.Figure7 shows the comparison between the simulation result and the 2D theoretical approximate solution. In this inclined crack and hole model, because of the existence of crack’s inclined angle, its tip’s SIFs K1 and K2 shall change. When we keep the distance between centers of the crack and the hole to D1 and change the crack’s inclined angle α from 0° to 90°, the K1 value shall decrease to 0 with the increase in the angle and there is no apparent difference between the two K2 values. Besides, they both increase first and then decrease. When the loading direction is vertical to the crack’s plane, the K1 value is the largest and the crack is most likely to start cracking. When the loading direction is parallel to the crack’s plane, there is basically no stress concentration.Figure 7
The comparison between the result from the 3D model in Figure6 and the 2D theoretical solution.


## 3. Results and Discussion

### 3.1. The Simulation of Graphite Model Containing a Large Number of Randomly Distributed Microcracks and Microvoids

Related research on graphite’s microstructure has shown that there are obvious microvoids and microcracks in nuclear graphite’s microstructure, especially the oxidized graphite material.In this section, we study the model where there are a large number of microvoids and microcracks in the elastic matrix. The influence of microvoids’ and microcracks’ distribution, size, density, and so forth on the performance of the overall damage model is studied. A graphite model containing 50 randomly distributed microcracks and 50 randomly distributed microvoids is built, as in Figure8. By changing the sizes of microvoids and microcracks, the damage model’s effective elastic modulus decreases almost linearly as ω1 and ω2 increase. Figures 9 and 10 show the relations between the normalized effective moduli Eeff/E0, Keff/K0, and ω2, with ω1 fixed.Figure 8
A global translucent view of element distribution in a cube containing randomly distributed microcracks and microvoids.Figure 9
E
e
f
f
/
E
0 versus the microcrack’s and microvoid’s density.Figure 10
K
e
f
f
/
K
0 versus the microcrack’s and microvoid’s density.From this simulation, we can see that, under different microvoid densities, the two effective moduliEeff and Keff decrease linearly with the increase in the microcrack’s density.Figures11 and 12 show ω1 and ω2’s combined influence on the microvoid and microcrack model. The two independent variables on the horizontal axes are ω1 and ω2, while two normalized effective moduli Eeff/E0 and Keff/K0 are on the vertical axis. Both figures show that, due to the different definitions of the microcrack and the microvoid, the decreasing rate of Eeff and Keff with the microcrack density increase is about twice as fast as that with microvoid density increase. In other words, under the same density increase, the increase in microcrack density shall make the whole model’s moduli decrease faster.Figure 11
E
e
f
f
/
E
0versus the microcrack’s and microvoid’s density.Figure 12
K
e
f
f
/
K
0 versus the microcrack’s and microvoid’s density.


### 3.2. The Simulation of Graphite Model Containing Microcracks and Growing Microvoids

On the basis of the simulated analysis of the graphite model containing randomly distributed microcracks and microvoids, the simulation of graphite model with the growing microvoids is conducted in order to study the effect of the microvoid’s growth on the existing microcrack during the graphite’s oxidization process.The model in Figure13 also contains 50 microcracks and 50 microvoids with random locations and directions. The microcrack keeps its geometrical parameters during the microvoid’s size change. The microvoid’s density ω1 and the microcrack’s density ω2 are shown in Figure 13.Figure 13
A model with a fixed microcrack density and growing microvoids.
(a)
ω
1
=
0.03;ω2=0.03
(b)
ω
1
=
0.03; ω2=0.075
(c)
ω
1
=
0.03; ω2=0.105
(d)
ω
1
=
0.03; ω2=0.15In order to study the effect of growing microvoids on the microcracks’ opening displacement and the tips’ SIFs under the same force boundary conditions, we simulate models with differentω2. All microcracks’ opening displacement is drawn in Figure 14. With the microvoid’s growth, most microcracks’ opening distances are gradually increasing, which means that with the increase in the microvoid’s size the microcrack bears more local stress concentration and is easier to be damaged, especially those microcracks with relatively coincident normal direction and loading direction. Their opening distances show a more apparent increasing tendency.Figure 14
The microcracks’ opening displacement.
(a)
ω
1
=
0.03; ω2=0.03
(b)
ω
1
=
0.03; ω2=0.075
(c)
ω
1
=
0.03; ω2=0.105
(d)
ω
1
=
0.03; ω2=0.15In this series of models, we partly select every microcrack tip node’s stress intensity factorK1 randomly. Figures 15 and 16 show that the value of K1 changes with the increase in the microvoid’s size. ω1.0 represents the maximal value of ω1 in a certain simulation.Figure 15
K
1 versus relative microvoid density when ω2=0.02.Figure 16
K
1 versus relative microvoid density when ω2=0.06.From the above models, it is clear that under single axis tension, the microcrack tips’K1 values are located in the range between 0 and 1 basically, which is also a result of the randomness of microcrack directions and microvoid distributions. With the increase in the microvoid’s size, the changing tendency of microcrack tip nodes’ SIFs also shows a large randomness: some microcrack tips’ SIFs shall increase in a relatively large scale because the distance between the microcrack and the microvoid decreases and the microcrack plane is vertical to the loading direction; some microcrack tips’ SIFs shall decrease or change randomly; many relatively normal microcrack tips’ K1 values shall remain or slowly increase.If we average all microcrack tip nodes’K1 values in the models with different microvoid sizes, the K1 values versus the increase in microvoids’ sizes are shown in Figure 17. This figure shows that the microcrack tips’ SIFs slowly increase with the microvoid growing.Figure 17
AverageK1 versus relative microvoid density.


## 3.1. The Simulation of Graphite Model Containing a Large Number of Randomly Distributed Microcracks and Microvoids

Related research on graphite’s microstructure has shown that there are obvious microvoids and microcracks in nuclear graphite’s microstructure, especially the oxidized graphite material.In this section, we study the model where there are a large number of microvoids and microcracks in the elastic matrix. The influence of microvoids’ and microcracks’ distribution, size, density, and so forth on the performance of the overall damage model is studied. A graphite model containing 50 randomly distributed microcracks and 50 randomly distributed microvoids is built, as in Figure8. By changing the sizes of microvoids and microcracks, the damage model’s effective elastic modulus decreases almost linearly as ω1 and ω2 increase. Figures 9 and 10 show the relations between the normalized effective moduli Eeff/E0, Keff/K0, and ω2, with ω1 fixed.Figure 8
A global translucent view of element distribution in a cube containing randomly distributed microcracks and microvoids.Figure 9
E
e
f
f
/
E
0 versus the microcrack’s and microvoid’s density.Figure 10
K
e
f
f
/
K
0 versus the microcrack’s and microvoid’s density.From this simulation, we can see that, under different microvoid densities, the two effective moduliEeff and Keff decrease linearly with the increase in the microcrack’s density.Figures11 and 12 show ω1 and ω2’s combined influence on the microvoid and microcrack model. The two independent variables on the horizontal axes are ω1 and ω2, while two normalized effective moduli Eeff/E0 and Keff/K0 are on the vertical axis. Both figures show that, due to the different definitions of the microcrack and the microvoid, the decreasing rate of Eeff and Keff with the microcrack density increase is about twice as fast as that with microvoid density increase. In other words, under the same density increase, the increase in microcrack density shall make the whole model’s moduli decrease faster.Figure 11
E
e
f
f
/
E
0versus the microcrack’s and microvoid’s density.Figure 12
K
e
f
f
/
K
0 versus the microcrack’s and microvoid’s density.


## 3.2. The Simulation of Graphite Model Containing Microcracks and Growing Microvoids

On the basis of the simulated analysis of the graphite model containing randomly distributed microcracks and microvoids, the simulation of graphite model with the growing microvoids is conducted in order to study the effect of the microvoid’s growth on the existing microcrack during the graphite’s oxidization process.The model in Figure13 also contains 50 microcracks and 50 microvoids with random locations and directions. The microcrack keeps its geometrical parameters during the microvoid’s size change. The microvoid’s density ω1 and the microcrack’s density ω2 are shown in Figure 13.Figure 13
A model with a fixed microcrack density and growing microvoids.
(a)
ω
1
=
0.03;ω2=0.03
(b)
ω
1
=
0.03; ω2=0.075
(c)
ω
1
=
0.03; ω2=0.105
(d)
ω
1
=
0.03; ω2=0.15In order to study the effect of growing microvoids on the microcracks’ opening displacement and the tips’ SIFs under the same force boundary conditions, we simulate models with differentω2. All microcracks’ opening displacement is drawn in Figure 14. With the microvoid’s growth, most microcracks’ opening distances are gradually increasing, which means that with the increase in the microvoid’s size the microcrack bears more local stress concentration and is easier to be damaged, especially those microcracks with relatively coincident normal direction and loading direction. Their opening distances show a more apparent increasing tendency.Figure 14
The microcracks’ opening displacement.
(a)
ω
1
=
0.03; ω2=0.03
(b)
ω
1
=
0.03; ω2=0.075
(c)
ω
1
=
0.03; ω2=0.105
(d)
ω
1
=
0.03; ω2=0.15In this series of models, we partly select every microcrack tip node’s stress intensity factorK1 randomly. Figures 15 and 16 show that the value of K1 changes with the increase in the microvoid’s size. ω1.0 represents the maximal value of ω1 in a certain simulation.Figure 15
K
1 versus relative microvoid density when ω2=0.02.Figure 16
K
1 versus relative microvoid density when ω2=0.06.From the above models, it is clear that under single axis tension, the microcrack tips’K1 values are located in the range between 0 and 1 basically, which is also a result of the randomness of microcrack directions and microvoid distributions. With the increase in the microvoid’s size, the changing tendency of microcrack tip nodes’ SIFs also shows a large randomness: some microcrack tips’ SIFs shall increase in a relatively large scale because the distance between the microcrack and the microvoid decreases and the microcrack plane is vertical to the loading direction; some microcrack tips’ SIFs shall decrease or change randomly; many relatively normal microcrack tips’ K1 values shall remain or slowly increase.If we average all microcrack tip nodes’K1 values in the models with different microvoid sizes, the K1 values versus the increase in microvoids’ sizes are shown in Figure 17. This figure shows that the microcrack tips’ SIFs slowly increase with the microvoid growing.Figure 17
AverageK1 versus relative microvoid density.


## 4. Conclusions

This paper adopted the 3D FM-DBEM to simulate the graphite model containing a large quantity of microvoids and microcracks. By the simulation under the approximate 3D equivalent plane strain state and the comparison to several existing theoretical approximate solutions, we have ensured the effectiveness and accurateness of the solver and proved that a result with high precision can be guaranteed by using this FM-DBEM to simulate the structural model containing microvoids and microcracks.There is an obvious strengthening effect at the microcrack tip nearest the microvoid’s edge, and the SIFs of the microcrack tips increase rapidly with the decrease in the distance between the microcrack and the microvoid, whereas the microcrack tips far from the microvoids’ edge are influenced slightly and the SIFs remain with the original values.Through the 3D simulation of the graphite model containing a large number of microvoids and microcracks, we find that the whole model’s two effective moduliEeff and Keff both decrease with the increase in the two microstructures’ density. The decreasing rate of Eeff and Keff with microcrack density increasing is about twice as fast as that with microvoid density. Meanwhile, with the growth of the microvoid, the microcrack shall bear more local stress concentration under the same load condition; therefore, there shall be larger crack opening displacement and the graphite is more likely to be damaged.



---

*Source: 1071709-2017-06-06.xml*