Image-based simulation of dispersed composites and porous rocks

. The article presents the results of numerical simulation of the effective elastic properties and stress concentration of the dispersed composite and porous rocks as well as the elastoplastic diagrams of tensile stress of the B 4 C / 2024Al composite using asymptotic averaging method. 2D and 3D simulations of uniaxial tension was carried out using the image-based models of the real structure obtained by X-ray tomography. We also used model structures of the B 4 C / 2024Al composite, consisting of an aluminum matrix with boron carbide inclusions in the form of ellipsoids with a random orientation. It turned out that in the elastic region 3D real structure and 3D ellipsoidal structure lead to the same e ﬀ ective moduli. On contrary, the stress concentrations are much higher in the real structure model than in the model of ellipsoidal inclusions. Also, two types of the models behave completely di ﬀ erently in the plastic region. Comparison of 2D and 3D calculations with experimental data shows that 2D models can be also used for determining the e ﬀ ective elastic moduli of B 4 C / 2024Al composite, but the use of 3D models is preferable for porous geomaterials and dispersed composites.


Introduction
Asymptotic averaging is a mathematically rigorous approach for determining the effective properties of heterogeneous and porous composites.Initially, the method was developed for determining the elastic moduli of periodic structures [1][2][3].However, this method can be used for non-periodic structures [4,5] and for determining the other properties, for example, Biot's tensor parameter of water-saturated porous rocks [6][7][8].In this case, it is important to define a representative volume element (RVE) and formulate boundary conditions at the boundary of the RVE.The first definitions of the RVE and effective elastic moduli as well as the method for their calculation were given by R. Hill [9,10].In our work, we used the generalized definition of a representative volume element formulated in [5].RVE is the volume for which the dependence of the effective property under consideration (for example, Young's modulus or Biot's coefficient) on the typical microstructure parameter (for example, inclusion volume ratio or porosity) lies on the same curve with acceptable accuracy.
In the first part of this paper, the solution of the local problems is performed for the B 4 C/2024Al composite in both elastic and plastic regions.The calculations were carried out for an image-based composite structure obtained using X-ray tomography.Further in this paper we call such a structure as the real one.For comparison, we also used the generated structure with ellipsoidal inclusions, further called as the model structure.Both material microstructures were used for subsequent finite element (FE) mesh generation.The effective stress-strain diagrams obtained using both structures were compared with the experimental results.The calculations took into account the progressive destruction of the matrix in the plastic region.
The second part of the article is devoted to the determination of the effective elastic moduli and Biot parameter of porous rocks.The development of X-ray tomography method and software for digitizing tomography images made it possible to obtain high-quality 3D models of the real structure of porous materials, which can be used for numerical experiments with the help of finite element codes.The results of 2D and 3D calculations were also compared with experimental data.

Study of stress-strain state in B 4 C/2024Al composite in computational uniaxial tension experiments
Aluminum alloys are widely used as a matrix for metal matrix composites because of its good machinability, wear resistance and mechanical properties.Boron carbide (B 4 C) is promising reinforcing inclusion in composites with an aluminum matrix due to its attractive properties such as ultra-high hardness, low density and good thermal and chemical stability [11].
Boron carbide inclusions improve the mechanical properties of the composite, such as elastic moduli, tensile strength, compressive strength, hardness and yield strength [12].As a result, aluminum matrix composites reinforced with boron carbide (B 4 C) particles are popular in automotive and aerospace industry, and especially in nuclear industry [11][12][13].The results of the experimental determination of the mechanical properties of metal matrix composites with reinforcing inclusions are widely presented in the literature.In addition to experimental studies, it is advisable to develop numerical modeling of such composites based on X-ray images.In our work, we studied the elastic/plastic properties of the dispersed B 4 C/2024Al composite depending on the volume ratio of boron carbide based on X-ray images of the internal structure (figure 1). Figure 1 1(c) shows a 2D model also built from X-ray images.As one can be seen, the inclusions have an irregular angular shape.In order to study the influence of the inclusion shape on the effective properties of the B 4 C/2024Al composite, we also used models in which boron carbide inclusions were modeled by ellipsoids with a random orientation.

Determination of elastic moduli and stress concentration in B 4 C/2024Al composite
The effective elastic moduli of dispersed B 4 C/2024Al composites were determined by the well-known asymptotic averaging method [3] based on 2D and 3D models of the real structure (figure 1) and 3D models with ellipsoids.The calculated values of Young's modulus were compared with experimental ones.According to the method mentioned above, it is necessary to solve a local problem in the RVE with a special boundary condition in the form of a linear function of coordinates on the external border Σ RVE [4,5]: Here, ξ are fast coordinates in the RVE, ξ q = x q /ε, ε = l/L, ε is a small geometric parameter, l is the characteristic linear size of the RVE, L is the characteristic linear size of the studied composite volume, n j and τ i are the components of normal and tangent vectors to Σ RVE , respectively.The following elastic properties were set in the calculations for the boron carbide inclusions: E = 450 GPa, ν = 0.16 and for the aluminum matrix: E = 74.7 GPa, ν = 0.32.Numerical experiments on uniaxial tension were carried out in a finite element (FE) program.Elastic moduli are calculated by the formulas: Figure 2 shows the results of calculating Young's modulus and comparison with the experimental data.As one can see in figure 2, for the B 4 C volume ratio range from 0.14 to 0.21, the dependence of the effective Young's modulus on the boron carbide volume ratio is linearly approximated with a fairly good accuracy for 3D and 2D models.Young's modulus values for a 3D real structure and a 3D model with ellipsoids differ from the experimental values by no more than 1-1.5% and 2-3%, respectively, at the same B 4 C volume ratio (figure 2).Thus, 3D models with ellipsoidal inclusions can be used for evaluating the effective elastic moduli, at least for the inclusion volume ratio up to 0.20.The values of Young's modulus for a 2D real structure are 6-8% lower than the experimental values at the same B 4 C volume ratio.Hence, it can be concluded that the use of 3D models is clearly preferable to 2D models.
We used 3D models of the real microstructure and 3D models with randomly oriented ellipsoidal inclusions to study the stress-strain state of the B 4 C/2024Al composite at the microlevel in the uniaxial tension and the stress concentration near boron carbide inclusions of various shapes [5].The equivalent Von Mises stress concentrations C Mises are considered as Here, σ 1 , σ 2 , σ 3 are the principal stresses, σ Mises is the equivalent Von Mises stress, σ MaxMises and σ AvMises are its maximum and average values, respectively.Table 1 shows the results of calculating the stress concentration on a mesh of nearly 2 million nodes [5].As can be seen from the table 1, the maximum concentration is observed in models of the real structure in boron carbide inclusions.In B 4 C inclusions, the von Mises stress concentration for the real structure reaches 29. 19, and for the model structure with ellipsoidsonly 5.65.In a less rigid aluminum matrix, the maximum von Mises stress concentration is significantly lower (6.19 for a real structure and 2.77 for a model with ellipsoids).Thus, studying the concentration of stress intensity, it is important to take into account the real structure of dispersed composites.

Elastoplastic properties of B 4 C/2024Al composite under uniaxial tension
The elastoplastic behavior of the B 4 C/2024Al composite was investigated in uniaxial tension computational and real experiments.Calculations were carried out taking into account the progressive failure in the matrix for 3D models of real (figure 3(a)) and model (with ellipsoids) (figure 3(b)) structures.Modeling the progressive damage in the material was performed as follows.In our problem, B 4 C inclusions are assumed to be non-destructive.At each load step the whole finite element in the matrix is considered to be damaged if the equivalent strain in it reaches 8%.The equivalent strain in the elements is calculated as the average value over all Gaussian points.In this case, the stiffness of each destroyed finite element is set to be zero.Then the next load step is carried out.
Figure 3(a) shows that for the model of the real structure accounting the progressive damage ensures the coincidence of the calculated and experimental data.As for the model with ellipsoidal inclusions (figure 3(b)), progressive damage does not lead to the coincidence of the calculated and experimental diagrams.On the other hand, the diagrams calculated without progressive failure technique included coincide with each other for both real and model inclusions but not with the experimental diagram.It should be mentioned that the numerical solution procedure with an embedded progressive damage algorithm becomes unstable at some stage of loading (figure 3(b)).

Determination of effective properties of porous rocks
Asymptotic averaging was also used to determine the effective elastic moduli and Biot's parameter of porous rocks based on 2D and 3D models of the real structure obtained using X-ray tomography.For the calculations, two types of porous rocks were considered: limestone and hyaloclastite.Calculated Young's moduli were compared with experimental data.Calculated Biot's coefficients were verified by the well-known formula, which utilizes bulk moduli of the porous rock and skeleton material.The dependences of Young's modulus and Biot's coefficient on porosity and pore shape were studied.

Theoretical aspects and solution method
Biot's parameter (pore pressure transfer tensor) is included in the formula for calculating effective stresses [15,16]: Here �p� is the fluid pressure (positive under compression).Angle brackets indicate averaging over the RVE.For isotropic rocks, Biot's parameter is a scalar coefficient that varies from 0 to 1 and depends on the rock properties (porosity, shape and orientation of pores).In [17,18], an empirical formula was proposed for calculating the isotropic pore pressure transfer coefficient α: where K s is the bulk modulus of the skeleton material, K e f f is the effective bulk modulus.
A rigorous mathematical derivation of formula ( 1) is presented in [19].The pore pressure transfer tensor (Biot's parameter) is difficult to determine experimentally.Moreover, the experimental determination has not been yet developed for anisotropic rocks.Therefore, the computational method for determining this parameter is relevant and of great practical importance.The asymptotic averaging method used to determine Biot's parameter is similar to the procedure for determining the effective elastic moduli and leads to the solution of local problems in the RVE of the porous rock [6][7][8].
To determine the effective elastic moduli of a porous medium in the RVE, it is necessary to solve a set of local boundary-value problems with special boundary conditions in the form of linear functions of coordinates on the outer RVE boundary Σ RVE and zero pore pressure p on both pores boundaries Σ int and Σ RVE [6][7][8]: Finally, elastic moduli are calculated by the formulas: If the material is isotropic, effective Young's modulus E e f f , Poisson's ratio ν e f f and bulk modulus K e f f are calculated using the formulas of elasticity theory (the superscript "eff" is omitted for ease of notation): Here λ and µ are the elastic Lame moduli.
For determining Biot's parameter, a boundary value problem is formulated in the RVE as follows.A constant fluid pressure is set on the pore boundary Σ int , while the outer RVE boundary Σ RVE is fixed in the perpendicular direction [6][7][8]: The last boundary condition replaces the periodicity condition and is rigorously correct, if the RVE has a symmetric shape.The components of Biot's tensor α i j are calculated by the formula: Certainly, this method is supposed to be implemented using a finite element code.

Examples of calculations
The calculation of the effective elastic moduli and components of Biot's tensor was carried out for two types of rocks such as dolomitic limestone and hyaloclastite.Dolomitic limestone core samples were taken in the south-eastern part of Moscow from wells in the depth interval from 15 to 38 m.Dolomitic limestone consists of calcite (92%), dolomite (6.5%), quartz and clay minerals (< 2%).Limestone has a microcrystalline structure, a homogeneous porous texture (figure 5(a)).The pores are evenly distributed, the pore size is from 0.1 to 1-2 mm.In the calculations, the material of the skeleton was considered homogeneous with elastic properties: E = 50 GPa, ν = 0.33, given according to literature data.
To study the stress-strain state of rocks at the microlevel and calculate the effective properties, 3D (figure 5(b)) and 2D (figure 5(c)) models of the real microstructure were used, obtained using a Yamato TDM-1000H-II X-ray computer microtomograph (Lomonosov Moscow State University) and digitization in the VG MAX 3.3 program in OOO Sovtest-Service (Kursk).Two-dimensional calculations were performed for a plane strain state.
Figure 6(a) shows the dependence of Young's modulus on porosity for dolomitic limestones based on the results of calculations and experiments.As one can see, the results obtained using 3D models are close to the experimental data, the difference is less than 5%.Young's modulus values calculated using 2D models differ from the results of 3D numerical calculations and experimental data by 20-30% at porosity from 7 to 28%.At low values of porosity (less than 5%), the values of Young's modulus for 2D and 3D models are close to each other.) that the results of 2D and 3D calculations by the averaging method are close to each other at low porosity (up to 5%).The values of Biot's coefficient for 2D models are 20-35% higher than for 3D models at porosity from 5 to 20%.The discrepancy becomes smaller and amounts to 10-15% with a porosity of more than 20%.As one can see (figure 6(b)), the results obtained by 3D calculations and given by the formula (1) coincide.
Hyaloclastites are volcanogenic-sedimentary rocks obtained from the southern and southwestern regions of Iceland.They were formed in the process of solidification and crystallization of magmatic melt during volcanic eruptions in underwater conditions [20].Slightly altered hyaloclastites consist of fragments of volcanic basalt glass.Hyaloclastites have got pores of a round (figure 7(a)) or angular (figure 7(b)) shape due to the peculiarities of the origin.Therefore, such geomaterials can be used to study the effect of pore shape on effective properties.The following elastic properties of the skeleton material were set in the calculations: E = 12.3 GPa, ν = 0.35.
We investigated the effect of pore shape on Young's modulus (figure 8(a)) and Biot's parameter (figure 8(b)) using 2D models of hyaloclastites.It turned out that for the same porosity in models with round pores, Young's modulus is higher (figure 8(a)), and Biot's coefficient, on the contrary, is lower (figure 8(b)) than in models with elongated angular pores.Probably, this is due to the fact that the contact area of the pores and the solid phase for round pores is smaller than for elongated pores.

Conclusion
The so-called local problems related to asymptotic averaging were numerically implemented in application to B 4 C/2024Al composite and porous geomaterials.Numerical implementation uses 3D structure images obtained by means of X-ray tomography.Commercial code VG MAX 3.3 was used for FE mesh generation and following porting to FE software.The FE calculations were also done for a model structure in which B 4 C inclusions are modeled by ellipsoids.A significant discrepancy in the computed effective stress/strain diagrams in the plastic strain range was revealed between the real structure of B 4 C/2024Al and the model structure with ellipsoids, if the progressive destruction of the matrix was taken into account.The calculated stress/strain diagram coincides with the experimental data quite well, if the image-based composite structure is exploited.
The results of calculating the effective elastic properties of the dispersed composite and porous rocks showed that 2D models can be used to determine the Young's modulus of B 4 C/2024Al composite, and for porous rocks 2D models are unsuitable and 3D models should be used instead.
The averaging method provides a general way to calculate Biot's tensor of porous rocks.Knowledge of the accurate assessment of the pore pressure transfer tensor is necessary for the correct determination of effective stresses.Adequate determination of Biot's tensor of porous water-saturated rocks by the asymptotic averaging method is possible only with the use of 3D models.A significant influence of the pore shape on the effective Young's modulus and Biot's parameter of porous rocks was revealed.

Figure 1 .
Figure 1.3D model (a), B 4 C inclusions (b) and 2D model (c) of the B 4 C/2024Al composite (a) shows a 3D two-component model of the B 4 C/2024Al composite constructed on the basis of X-ray images obtained at the Shanghai Institute of Applied Physics and processed in the VG MAX 3.3 program [14] at Sovtest-Service LLC (Kursk).The images were supplied by Prof. Qiang Zhang from Harbin Institute of Technology, China in framework of collaborative project.In figure 1(b) B 4 C inclusions are depicted, and figure

Figure 2 .
Figure 2. Dependence of Young's modulus on B 4 C volume ratio: 1) experiments, 2) 3D model of the real structure, 3) 3D model with ellipsoidal inclusions, 4) 2D model of the real structure

Figure 7 .Figure 8 .
Figure 7. Photographs of thin sections of hyaloclastite-with round (a) and angular (b) pores

Table 1 .
The equivalent Von Mises stress concentration in computational uniaxial compression experiments