3D analysis of pore effect on composite elasticity by means of the finite-element method

We developed a 3D buffer-layer finite-element method model to investigate the porosity effect on macroscopic elasticity. Using the 3D model, the effect of pores on bulk effective elastic properties was systematically analyzed by changing the degree of porosity, the aspect ratio of the ellipsoidal pore, and the elasticity of the material. The results in 3D space were compared with the previous results in 2D space. Derivatives of normalized elastic stiffness constants with respect to needle-type porosity were integers, if the Poisson ratio of a matrix material was zero; those derivatives of normalized stiffness elastic constants for C33, C44, C11, and C66 converged to −1, −2, −3, and −4, respectively, at the corresponding condition. We have developed a criterion of R < ∼1∕3, where the mutual interaction between pores became negligible for macroscopic composite elasticity. These derivatives were nearly constant at less than 5% porosity in the case of a spherical pore, suggesting that the interaction between neighboring pores was insignificant if the representative size of the pore was less than one-third of the mean distance between neighboring pores. The relations we obtained in this work were successfully applied to invert the bulk modulus and rigidity of Cmcm-CaIrO3 as a case study; the performance of the inverting scheme was confirmed through this assessment. Thus, our scheme is applicable to predict the macroscopic elasticity of porous object as well.


INTRODUCTION
Porous systems are ubiquitous.Sandstone and pumice are porous rocks.Karst terrain is representative of porous structures on the earth (e.g., Hiltunen et al., 2007), and most meteorites and asteroids also have significant porosity (Britt et al., 2002).Aside from materials associated with earth and planetary sciences, some biological tissues, such as those composing trunks and bones, are porous.Furthermore, there are many porous industrial goods, such as rubber sponges, metal foams, and semisintered ceramics.Therefore, the macroscopic physical properties of porous objects have been studied intensively not only in earth sciences but also in other disciplines, such as material engineering.
Composite elasticity of porous objects is of fundamental interest because of the nontrivial interactions among pores.Although a single spherical pore in an isotropic elastic body can be analyzed explicitly by exact analytical solutions, interactions among multiple pores cannot be solved analytically.Therefore, the porosity effect on composite elasticity has been investigated by various bound methods, such as the Voigt-Reuss bound (Hill, 1952) and the Hashin-Shtrikman bound (Hashin and Shtrikman, 1963).However, it is well known that these bounds cannot constrain the composite elasticity of porous materials in a reasonably narrow range; the lower bounds are insignificant owing to the infinitely small elasticity of pores (e.g., Watt, 1976).
The differential effective medium (DEM) method, a kind of self-consistent approach, has been applied to the problem of porous materials (e.g., Budiansky, 1965;Hill, 1965;Wu, 1966;Walpole, 1969;Guéguen et al., 1997).The validity of the DEM method has been checked using experimental data for specimens of sintered hard material with porosity.However, it has been difficult to prepare porous specimens with sufficient quality for evaluating the performance of these theoretical approaches.Thus, the finite-element method (FEM) was introduced to provide unambiguous references (e.g., Berryman, 1995;Mavko et al., 1998;Garboczi and Berryman, 2001).
We examined the composite elasticity of porous objects by a 3D FEM.An important advantage of the FEM approach is that it can deal with realistic configurations without difficulty; there have been extensive studies on "digital rock," i.e., analog models of natural rock synthesized in a computer by using the FEM and other simulation techniques (e.g., Torquato, 2002;Chen et al., 2007;Nouy and Clement, 2010;Cho et al., 2013).Furthermore, we can extend the present FEM analysis to anisotropic cases in terms of not only physical properties of the constituent material but also the geometric configuration of the pores in 3D space.
The present work is a 3D extension of the previous 2D buffered layer model analysis (Yoneda and Sohag, 2011).The uniqueness of the buffered layer model is its regular distribution of pores, which is in distinct contrast to the preceding works based on relatively irregular geometry (e.g., Cho et al., 2013;Matthies, 2015).Needless to say, its regularity is definitely more advantageous to evaluate the porosity effect systematically.
We have succeeded to obtain various interesting results through 3D FEM analysis.Furthermore, the results were successfully applied to derive intrinsic elasticity of Cmcm-CaIrO 3 from its nominal elasticity measured on a porous sample.Cmcm-CaIrO 3 is a wellknown analog of the postperovskite (PPv) phase in MgSiO 3 (Murakami et al., 2004;Oganov and Ono, 2004).It is much more elastically anisotropic than Pbnm-CaIrO 3 of the bridgmanite analog; it can also be a good analog of crustal anisotropic minerals such as quartz, feldspar, and serpentine (e.g., Vasin et al., 2013;Watanabe et al., 2014).
In this work, we limited the present FEM analysis to only a completely unsaturated dry pore system, although the degree of liquid saturation in the pores is an important factor for the macroscopic elasticity of porous rock (e.g., Takei, 2002;Vasin et al., 2013).Although it is noted that FEM is inappropriate to simulate nonstationary fluid flow inside pores, the present result is beneficial for examining porous systems with fluid because it can constrain the macroscopic elastic property of the porous host rock accommodating water or oil.
Finally, we follow the Voigt notation for elasticity (stress, σ; strain, ε; stiffness elastic constant, C; and elastic compliance constant, S) throughout this study (e.g., Nye, 1985); it is noted that the Voigt notation is a conventional way to express the exact elastic constants of a fourth-rank tensor.

PROCEDURES FOR 3D FINITE-ELEMENT METHOD ANALYSIS
The procedures for 3D FEM analysis are a natural extension of those for 2D analysis by Yoneda and Sohag (2011).Because the fundamental concepts and procedures have already been described, we are restricted to a brief description of the procedures of the present 3D study.
Figure 1 shows an octant part of the 5 × 5 × 5 buffered layer model used in the present 3D analysis; the figures in Appendix B show some extreme cases with mesh density.Although the outermost edge length is set as 1 m in the present work for simplicity, the absolute size of the geometry does not matter for the present purpose.This means that the ratio between the edge length and the diameter of the pore is essential.Considering the symmetry along the central interfaces, FEM analysis can be conducted only on an octant part.In this work, we conducted FEM analysis using COM-SOL®, which is commercial FEM software.
Figure 2 shows the ellipsoidal shape of pores discussed in this study.The three principal axes of the ellipsoid were set as a, b, and c in the x-, y-, and z-directions, respectively.To decrease the complexity, we assume a ¼ b, and we adjust að¼ bÞ and c.Here, the inverse of the aspect ratio α 0 ð¼ a∕cÞ is introduced to describe the present results; note that the definition of α 0 is inverse of the aspect ratio that is usually used in poroelasticity.
Here, we summarize the essence of the present model geometries: (1) Pores are embedded in an elastically isotropic matrix.(2) Pores are all identical, and (3) their centers distribute regularly at grid points in a model geometry.(4) Symmetric axes of pores align always in the z-axis direction.
Because of the equivalence between the xand y-directions, the porous bodies analyzed in this study have tetragonal symmetry.In the spherical pores cases, it has cubic symmetry because of the Figure 2. Examples of the shape of a pore as a function of aspect ratio α 0 .In this drawing, the size of the pores is normalized to have the same volume.

L16
Yoneda and Sohag regular distribution of pores.The terminologies of "tetragonal" and "cubic" symmetry for investigated bodies are based on an analogy of those symmetries defined in crystallography.
After setting the model geometry as described above, uniform displacements were applied to the outer surfaces of the model geometry to generate a stress-strain relation inside the model.In this work, the matrix was first assumed to be elastically isotropic with a Young modulus Y ¼ 100 GPa; the Poisson ratio ν was treated as a variable between 0 and 0.45.
Let us define u, v, and w as displacements in the x-, y-, and zdirections, respectively.In the 3D analysis, we can classify six types of forced displacement on the outermost surfaces as follows: (1) There is uniform forced displacement of Δu ¼ −1 × 10 −6 m in the x-direction on the x-surface, whereas the model is free in the yand z-directions, or there is no constraint on the yand z-surfaces.Cases (2) and (3) are similar to (1), but the directions of compression are in the yand z-directions, respectively.(4) There is a coupled uniform forced displacement of Δv ¼ 1 × 10 −6 m in the y-direction on the z-surface, and that of Δw ¼ 1 × 10 −6 in the z-direction on the y-surface.This is a pure shear deformation observed from the x-direction.Cases ( 5) and ( 6) are similar to (4), but corresponding pure shear deformations are observed from the yand z-directions, respectively.
We easily obtain the average strains from the resulting displacements Δu, Δv, and Δw along the observation interfaces: (1) where the overbar for Δu, Δv, and Δw specifies their average, and x 0 , y 0 , and z 0 are the positions of the observation interfaces in the x-, y-, and z-coordinates, respectively.We then formulate the simultaneous equations between the averaged stresses over the observed interface and those averaged strains: We can make three sets of simultaneous equations according to the uniaxial compressions in the x-, y-, and z-directions.In general, we have nine redundant equations for four unknown parameters, C 11 ð¼ C 22 Þ, C 12 , C 13 ð¼ C 23 Þ, and C 33 .These four parameters were solved using the least-squares method.For two independent shear elastic constants C 44 ð¼ C 55 Þ and C 66 , we have the following equations: (3) Here, we introduce normalized moduli C Ã ð¼ C∕C 0 Þ, where C 0 is the modulus of the matrix material itself or that of the porous system at zero porosity.According to the experience of the previous 2D analysis (Yoneda and Sohag, 2011), we would mainly observe the initial slope of the normalized elastic constant against porosity ϕ, defined as The usefulness of C Ã and D were well confirmed in the previous 2D analysis (Yoneda and Sohag, 2011).

RESULTS AND DISCUSSION
We start with the case of a spherical pore in a matrix with a Poisson ratio of ν ¼ 0.3 as the simplest representative case.Figure 3 shows the variations in the longitudinal and shear elastic constants (C P and C S ) as functions of porosity obtained using the FEM and DEM methods, respectively.Although the DEM calculation was conducted to 100% porosity, the FEM calculation was terminated at 50% porosity to avoid mutual contact between neighboring pores at 0.5236, which is the volume ratio between a cube and its inscribed sphere.We recognize that the porous object still has substantial stiffness to the limit of contact between neighboring spherical pores.This contrasts with the 2D analysis, in which the object loses stiffness when approaching the limit of circular pore contact.This is an obvious dimensional effect between the 2D and 3D analyses.
In the case of a spherical pore, the macroscopic object maintains cubic symmetry, where The general trend of the FEM results is similar to that estimated by the DEM method.The DEM results are consistent with the most compliant elastic constant in the present FEM analysis.This is consistent with the preliminary results obtained in the 2D analysis (Figure 14 in Yoneda and Sohag, 2011).
Figure 3 also shows that the directional fluctuations of C P and C S are approximately 20% and 33%, respectively, at approximately 50% porosity.The directional anisotropy is smaller than those in 2D analysis, where we recognized more than 25% and 75% directional fluctuation, respectively, at approximately 50% porosity (Figure 5 in Yoneda and Sohag, 2011).
Further, we recognize that the macroscopic elasticity of the porous object is nearly isotropic and has a linear change in the low-porosity range (< 5%) for the longitudinal and shear stiffness constants.This observation is analogous to the case in 2D analysis.This finding suggests an important concept that the pore effect can be treated as a "linear system" in the porosity range lower than approximately 5%.In other words, we expect "additivity" among pores with different shapes, sizes, and orientations as long as the distribution of pores is homogeneous and the porosity is low enough.Let us define the ratio R ¼ l∕d (l, a representative size of a pore, and d, the distance between the nearest pores); R is introduced as a scale for evaluating mutual interaction among pores.We conclude that for R < ∼1∕3 (approximately 0.37, or the cubic root of 0.05), where the mutual interaction between pores is negligible for macroscopic composite elasticity.If the shape and distribution of pores are anisotropic and/or heterogeneous, we should use a maximum l and a minimum d when evaluating R.
Figure 4 compares D values derived using the FEM and DEM methods in the spherical pore case.We can see that the present FEM results are consistent with those of the DEM method, as suggested in previous works (e.g., Matthies, 2015).

3D pore effect on composite elasticity
We recognize that the 5 × 5 × 5 model shows better consistency with DEM results than does the 7 × 7 × 7 model.This situation is caused by the limitation of machine capacity making a precise mesh in the 7 × 7 × 7 model.Therefore, we use the 5 × 5 × 5 model because of the better consistency with the DEM method in this work.
According to an analogy with the previous 2D analysis, we first check the relation between the D values and aspect ratio α 0 at ν ¼ 0. We find the following from Figure 5: 1) D 33 converges to −1 as α 0 → 0. This is consistent with intuition because the pores approach tubes extending in the z-direction at the limit of α 0 → 0. 2) D 11 and D 66 converge to −3 and −4, respectively, as α 0 → 0. This is consistent with the previous 2D analysis because of the geometric similarity.diverge to −∞ as α 0 → ∞.This is consistent with intuition because the pores approach plates in the x − y space at the limit of α 0 → ∞. 5) At α 0 ¼ 1, all the D ij values are close to −2, which is consistent with previous analysis for selected moduli by means of the selfconsistent method, DEM, etc. (e.g., Guéguen et al., 1997;Cho et al., 2013;Matthies, 2015).
Although our findings (equations 2 and 3) are difficult not only to understand intuitively but also to derive deductively, they are inter-esting and useful to qualitatively predict the macroscopic elasticity of porous object.
In the previous 2D analysis of a circular pore, we found useful relations consistent with 2D FEM results, such as those for a circular pore: We tried to find similar relations for 3D FEM analysis, starting with the analogous functional forms used in 2D analysis.We assumed the following functional forms for longitudinal and shear constants: 3D pore effect on composite elasticity D 44 ðα 0 ; νÞ ¼ −r 1 ðα 0 Þð1 − r 2 ðα 0 ÞνÞ r 3 ðα 0 Þ ; D 66 ðα 0 ; νÞ ¼ −s 1 ðα 0 Þð1 − s 2 ðα 0 ÞνÞ s 3 ðα 0 Þ : (7) Figure 6 shows the results of fitting equations 6 and 7 to the FEM results.We recognize that the functions reproduced the FEM results reasonably.The differences between the FEM results and the fitting values are within approximately 1.0 for D 11 and D 33 , and approximately 0.1 for D 44 and D 66 , for the fitting range of ν from 0.0 to 0.4.It is worth mentioning that the fitting performance improves to approximately 0.2 for D 11 and D 33 , and approximately 0.1 for D 44 and D 66 , for the fitting range of ν from 0.0 to 0.35. Figure 7 shows the resulting variations in the parameters p, q, r, and s with aspect ratio α 0 for the fitting range of ν from 0.0 to 0.4.It is worth mentioning that the global features of p, q, r, and s are similar for the fitting range of ν from 0.0 to 0.35 as well.
Finally, the fittings for D 12 and D 23 are mentioned.According to analogy with the previous 2D analysis, we tested the following functional forms to find a suitable functional form inductively: Figure 8 compares the results of fitting equations 8 and 9 to the FEM results.We see that the performance of the functions is not satisfactory except for α 0 ¼ 1.0 and ν of 0.1-0.3.This is in contrast with the previous 2D study, where an equation analogous with equations 8 and 9 worked very well (see equation 20 and Figure 10 in Yoneda and Sohag, 2011).The survey for more suitable fictional forms in 3D remains as a future subject.
The data sets of the D values are given in Appendix A, and we recommend that readers carry out their own fitting of the data around specific values of ν and α 0 using equations 6-9, a spline function, polynomial function, or another function.
Before the single crystal elasticity of Cmcm-CaIrO 3 was available, ultrasonic velocity was measured on a porous aggregate of Cmcm-CaIrO 3 synthesized at 8 GPa and 1200°C in the Kawai-type high-pressure apparatus.The specimen was confirmed to be a single phase of Cmcm-CaIrO 3 by means of X-ray diffraction.Figure 9 is a sectional image of the specimen showing the pore shape (α 0 ∼ 0.3) and distribution; its porosity was estimated to be 6% by comparing its nominal density and the X-ray density of 8211 kg∕m 3 (Sugahara et al., 2008).
We measured ultrasonic velocities in the three directions of a rectangular specimen (approximately 1 mm edge length), and we found it to be nearly isotropic within 2% in P-and S-wave velocities.
Figure 7. Resulting parameters p, q, r, and s as functions of aspect ratio α 0 , after the fitting of the FEM results based on equations 6 and 7 (see Figure 6).The solid, dashed, and dotted lines correspond to subscripts of 1, 2, and 3 for p, q, r, and s, respectively.A significant anomaly for r 2 (the dashed curve in the D 44 plot) at the left end of the plot area seems to be due to numerical instability.2 with nominal values of bulk modulus and rigidity; the two experimental data agree well with each other despite the 2% porosity difference.
The experimental data are enough to constrain intrinsic K and G from the four parameters of the porous object (porosity, aspect ratio of pore, and Pand S-wave velocities) as suggested in previous research (e.g., Takei, 2002).Thus, we tried to correct the nominal values of the bulk modulus and rigidity to the intrinsic ones by using the results of the present FEM analysis.
The present target specimen was assumed to be a macroscopically isotropic elastic object, whereas the present FEM porosity effect analysis was conducted on an object with tetragonal symmetry (equivalence between the xand y-directions) resulting in D 11 ¼ D 22 , D 13 ¼ D 23 , and D 44 ¼ D 55 .We have to take that difference into account in conducting the porosity effect correction.
The correction of the porosity effect was conducted using the following method: The initial values for C ðiÞ 44 , at the ith iteration.The factor of "0.02" in equation 10 is one-third of the 6% porosity, which was equally divided into the three directions of the x-, y-, and z-axes.In the present iteration, ν ðiÞ was slightly shifted from 0.285 to 0.287 throughout.
Figure 10 shows the results of the correction for the porosity effect.We found that K ¼ 177ð6Þ GPa and G ¼ 85ð2Þ GPa, where the magnitude of error was estimated from the uncertainty in aspect ratio α 0 , porosity, and the anisotropy of the ultrasonic velocities.We recognize that the present results are consistent with the intrinsic bulk modulus and rigidity (178 and 84 GPa as shown in Table 1) within the uncertainties.Thus, the present case study confirms the reliability of the present correction procedure for the porosity effect.
Figure 9. Scanning electron microscope image of the Cmcm-CaIrO 3 specimen with 6% porosity.The black parts correspond to pores, suggesting a tubelike geometry.Although it is difficult to determine the average aspect ratio α 0 of pores from the image, we used α 0 ∼ 0.3 as the starting value of the correction for the porosity effect.

CONCLUSION
Three-dimensional FEM analysis was conducted to study the porosity effect on macroscopic elastic properties.It was confirmed that the FEM results are consistent with the DEM results in the case of a low-porosity spherical pore.
We found that D 33 , D 44 , D 11 , and D 66 converge to −1, −2, −3, and −4, respectively, at a Poisson ratio of ν ¼ 0 and an aspect ratio of α 0 →0.We proposed a criterion of R < ∼1∕3, where the mutual interaction between the pores becomes negligible for macroscopic composite elasticity.
We proposed and examined functional forms to fit the present D values.Finally, the concept was applied to experimental data for a porous sintered specimen of Cmcm-CaIrO 3 .The resulting bulk modulus and rigidity is consistent with those calculated from single crystal elastic constants determined by IXS.This results convinced us to systematically evaluate the macroscopic elasticity of a porous object as a function of porosity and the aspect ratio of the pore.

ACKNOWLEDGMENTS
This work was supported by a Grant-in Aid for Scientific Research (nos. 19204044, 22224008, and 15H02128) from the Japan Society for the Promotion of Science.

APPENDIX A THE NUMERICAL DATA SET OF D IJ
The numerical data sets of D ij as a function of aspect ratio of pore (α 0 ) and Poisson ratio (ν) are given in Tables A1-A6.3D pore effect on composite elasticity

Figure 1 .
Figure 1.Octant part of a multi-buffer-layer model for N ¼ 5, or a 5 × 5 × 5 model with a cube having an edge length of 1 m.The edges of the central cube, indicated by short arrows, correspond to the observation interfaces for stress and displacement in the deformation tests.

Figure 4 .
Figure 4.The D 11 ð¼ D 22 ; D 33 Þ and D 44 ð¼ D 55 ; D 66 Þ versus thePoisson ratio ν for a spherical pore, or aspect ratio α 0 ¼ 1.The dashed lines are the results obtained using the DEM method.Symbols "○" and "▪" are used to show the results of the 5 × 5 × 5 and 7 × 7 × 7 FEM models, respectively.Note that the result obtained using the 7 × 7 × 7 model at ν ¼ 0.45 is not shown because of divergence in the FEM calculation.

Figure 8 .
Figure 8. D 12 versus the Poisson ratio ν at aspect ratio α 0 ¼ 0.1, 1.0, and 10.The solid lines are calculated using equation 8, and the dashed lines are calculated using equation 9, and the symbols "•" and "○" show the D 12 and D 23 obtained in FEM analysis.

Figure 10 .
Figure 10.Reconstructed elasticity of Cmcm-CaIrO 3 as a function of aspect ratio α 0 after the correction for the porosity effect.The solid and dashed curves are based on the present experimental results and those by Liu et al. (2011), respectively.The thick dotted rectangles show the magnitude of uncertainty around the central values shown by solid circles; it is caused by uncertainty of velocity, porosity, and the aspect ratio of the pore.The solid squares show the intrinsic (i.e., no-porosity) values determined by IXS measurement.(a) Bulk modulus versus aspect ratio.The thicker curve is C 11 − ð4∕3ÞC 44 , and the thinner curve is ð1∕3ÞðC 11 þ 2C 12 Þ.(b) Rigidity versus aspect ratio.The thicker curve is 0.5ðC 11 − C 12 Þ, and the thinner one is C 44 .
IN THE PRESENT FINITE-ELEMENT METHOD ANALYSISExamples of mesh in the present FEM analysis are shown in FiguresB1-B3.

Figure B- 1 .
Figure B-1.Mesh for the case of large spherical pores.Porosity is 0.45.The total element number is 77,156.

Figure B- 2 .
Figure B-2.Mesh for the case of needle-type pores (α 0 ¼ 0.1).Porosity is 0.001.The total element number is 62,411.The left-side image is the object before meshing.

Table A -
1.The D 11 α 0 ; ν values obtained in the present 3D FEM analysis.Table A-2.The D 33 α 0 ; ν values obtained in the present 3D FEM analysis.Table A-4.The D 66 α 0 ; ν values obtained in the present 3D FEM analysis.

Table A -
5. The D 12 α 0 ; ν values obtained in the present 3D FEM analysis.Note that D 12 cannot be calculated at ν 0.

Table A -
6.The D 23 α 0 ; ν values obtained in the present 3D FEM analysis.Note that D 23 cannot be calculated at ν 0.