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.
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
Figure 1 shows the model of an elastic structure containing edge, embedded voids, and cracks. and represent the domain and outer boundary of the solid. represent the inner boundary of the voids. and represent the two coinciding surfaces of the cracks.
The dual boundary integral equations for 3D elastic fracture mechanics without body force are shown aswhere and represent the source and field point; and are the boundary displacement and traction vectors; is a free term related to the shape of boundary at point ; is the relative displacement between and . The kernel functions and are shown aswhere represents the distance between and ; and represent the shear modulus and Poisson’s ratio. and are obtained from and :where and represent the outward normal vectors at points and ; 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 is obtained by the DBEM analysis. The SIFs at geometry point are given:where represents the distance between and ; and are the normal and tangent components of at under the local coordinate system defined at (see Figure 2).
2.2. Fast Multipole-Dual Boundary Element Method
Considering a reference point close to , if a threshold of is satisfied, the integrals in (1) and (2) can be expanded around aswhere, and are called multipole moments presented by Wang and Yao [37]; is called solid spherical harmonic function [39].
Equation (3) can be expanded asConsidering another reference point close to , if a threshold is satisfied, the integrals in (1), (2), and (3) can be expanded around in a local expansion form:where, and 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. 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 Figure 3. Based on previous circular and elliptical microcrack studies in BEM [41], this paper inherits most of the symbols and characters from previous research.
The formula for the constitutive relation of elastic problem without body force iswhere is the elastic strain, is the elastic coefficients, and the elastic strain field is defined bywhere is the elastic displacement.
The boundary conditions are given by mechanical boundary conditions:where is the surface traction and is the component of the unit outward normal vector of the surface.
The side length of the cubic is . For the spheroidal microvoid in this solid, if its radius is , then its volume is defined as . If the number of spheroidal microvoids is , the spheroidal microvoid density parameter in the solid is defined by . Similarly, the radius of the circular microcrack is and the number of elliptical microcracks is . Following Budiansky and O’connell [42] and Kachanov [43], the microcrack density parameter is defined as . Young’s modulus and Poisson’s ratio 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 modulus . The average strain and stress in the solid are defined byWhen different boundary conditions are applied, the material constants and effective moduli can be obtained.
Applying the boundary condition as inthen effective Young’s modulus in direction can be obtained byApplying the boundary condition as inwe can obtain
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 Figures 4 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 , the crack length is , , and the vertical distance from the crack tip to the hole’s margin is . 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.
(a) The schematic diagram of the 2D model
(b) The geometry and mesh of 3D model
Figure 5 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 is small. The crack tip to the left of the hole enlarges fast with the decrease of distance . The tip’s stress intensity factor shall become 2-3 times larger than when there is no hole and but only the single crack.
(a) The schematic diagram of the 2D model
(b) The geometry and mesh of 3D model
Moreover, 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 , the crack’s length maintains , and their distance is , which is equivalent to in the above model.
Figure 7 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 and shall change. When we keep the distance between centers of the crack and the hole to and change the crack’s inclined angle from 0° to 90°, the value shall decrease to 0 with the increase in the angle and there is no apparent difference between the two values. Besides, they both increase first and then decrease. When the loading direction is vertical to the crack’s plane, the 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.
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 Figure 8. By changing the sizes of microvoids and microcracks, the damage model’s effective elastic modulus decreases almost linearly as and increase. Figures 9 and 10 show the relations between the normalized effective moduli , , and , with fixed.
From this simulation, we can see that, under different microvoid densities, the two effective moduli and decrease linearly with the increase in the microcrack’s density.
Figures 11 and 12 show and ’s combined influence on the microvoid and microcrack model. The two independent variables on the horizontal axes are and , while two normalized effective moduli and are on the vertical axis. Both figures show that, due to the different definitions of the microcrack and the microvoid, the decreasing rate of and 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.
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 Figure 13 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 and the microcrack’s density are shown in Figure 13.
(a) ;
(b) ;
(c) ;
(d) ;
In 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 . 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.
(a) ;
(b) ;
(c) ;
(d) ;
In this series of models, we partly select every microcrack tip node’s stress intensity factor randomly. Figures 15 and 16 show that the value of changes with the increase in the microvoid’s size. represents the maximal value of in a certain simulation.
From the above models, it is clear that under single axis tension, the microcrack tips’ 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’ values shall remain or slowly increase.
If we average all microcrack tip nodes’ values in the models with different microvoid sizes, the 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.
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 moduli and both decrease with the increase in the two microstructures’ density. The decreasing rate of and 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.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research is supported by Tsinghua University Initiative Scientific Research Program (Grant no. 20131089332). Additional funds were also provided by the National Natural Science Foundation of China (Grant no. 11002080) and the China Scholarship Council.