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

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-CaIrO3$ from its nominal elasticity measured on a porous sample. $Cmcm-CaIrO3$ is a well-known analog of the postperovskite (PPv) phase in $MgSiO3$ (Murakami et al., 2004; Oganov and Ono, 2004). It is much more elastically anisotropic than $Pbnm-CaIrO3$ 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.

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 COMSOL®, 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 $α′(=a/c)$ is introduced to describe the present results; note that the definition of $α′$ 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 $x$- and $y$-directions, the porous bodies analyzed in this study have tetragonal symmetry. In the spherical pores cases, it has cubic symmetry because of the 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=100GPa$⁠; 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 $z$-directions, 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−6m$ in the $x$-direction on the $x$-surface, whereas the model is free in the $y$- and $z$-directions, or there is no constraint on the $y$- and $z$-surfaces. Cases (2) and (3) are similar to (1), but the directions of compression are in the $y$- and $z$-directions, respectively. (4) There is a coupled uniform forced displacement of $Δv=1×10−6m$ 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 $y$- and $z$-directions, respectively.

The usefulness of $C*$ and $D$ were well confirmed in the previous 2D analysis (Yoneda and Sohag, 2011).

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 (⁠$CP$ and $CS$⁠) 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 $C11=C22=C33$⁠, $C12=C23=C31$⁠, and $C44=C55=C66$ are satisfied. 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 $CP$ and $CS$ 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).

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 $α′$ at $ν=0$⁠. We find the following from Figure 5:

1. $D33$ converges to $−1$ as $α′→0$⁠. This is consistent with intuition because the pores approach tubes extending in the $z$-direction at the limit of $α′→0$⁠.

2. $D11$ and $D66$ converge to $−3$ and $−4$⁠, respectively, as $α′→0$⁠. This is consistent with the previous 2D analysis because of the geometric similarity.

3. $D44$ converges to $−2$ as $α′→0$⁠. This is a new finding of the present 3D analysis. Consequently, $D33$⁠, $D44$⁠, $D11$⁠, and $D66$ converge to $−1$⁠, $−2$⁠, $−3$⁠, and $−4$⁠, respectively, at $ν=0$ and as $α′→0$⁠.

4. $D11$ and $D66$ converge to $−1$ as $α′→∞$⁠, whereas $D33$ and $D44$ diverge to $−∞$ as $α′→∞$⁠. This is consistent with intuition because the pores approach plates in the $x−y$ space at the limit of $α′→∞$⁠.

5. At $α′=1$⁠, all the $Dij$ values are close to $-2$⁠, which is consistent with previous analysis for selected moduli by means of the self-consistent 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 interesting and useful to qualitatively predict the macroscopic elasticity of porous object.

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 $D11$ and $D33$⁠, and approximately 0.1 for $D44$ and $D66$⁠, 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 $D11$ and $D33$⁠, and approximately 0.1 for $D44$ and $D66$⁠, 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 $α′$ 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.

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 $α′=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 $α′$ using equations 6–9, a spline function, polynomial function, or another function.

$Cmcm-CaIrO3$ has been frequently investigated as a representative analog of unquenchable $PPv-MgSiO3$ (e.g., Hunt et al. 2009), since Murakami et al. (2004) and Oganov and Ono (2004) discover the PPv phase in $MgSiO3$ with $Cmcm-CaIrO3$ structure. Yoneda et al. (2014) report single crystal elasticity of $Cmcm-CaIrO3$ by means of the inelastic X ray scattering (IXS) technique (Table 1).

Before the single crystal elasticity of $Cmcm-CaIrO3$ was available, ultrasonic velocity was measured on a porous aggregate of $Cmcm-CaIrO3$ 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-CaIrO3$ by means of X-ray diffraction. Figure 9 is a sectional image of the specimen showing the pore shape (⁠$α′∼0.3$⁠) and distribution; its porosity was estimated to be 6% by comparing its nominal density and the X-ray density of $8211kg/m3$ (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. Liu et al. (2011) do independent study on $Cmcm-CaIrO3$ aggregates synthesized by themselves (8% porosity); those results are summarized in Table 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 $P$- and $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 $x$- and $y$-directions) resulting in $D11=D22$⁠, $D13=D23$⁠, and $D44=D55$⁠. 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 $C11(0)(=C22(0),C33(0))$⁠, $C12(0)(=C23(0),C31(0))$⁠, and $C44(0)(=C55(0),C66(0))$ were calculated to be 251, 101, and 75 GPa from the present 6% porosity object; the superscript suffix “0” indicates the initial value of the iterative process. It is noted that $C11$⁠, $C12$⁠, and $C44$ are equivalent with $λ+2μ$⁠, $λ$⁠, and $μ(=G)$⁠, respectively, by using Lame’s elastic constants for isotropic elasticity.

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 $α′$⁠, 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.

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 $D33$⁠, $D44$⁠, $D11$⁠, and $D66$ converge to $−1$⁠, $−2$⁠, $−3$⁠, and $−4$⁠, respectively, at a Poisson ratio of $ν=0$ and an aspect ratio of $α′→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-CaIrO3$⁠. 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.

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.

The numerical data sets of $Dij$ as a function of aspect ratio of pore (⁠$α′$⁠) and Poisson ratio (⁠$ν$⁠) are given in Tables A1–A6.

Examples of mesh in the present FEM analysis are shown in Figures B1–B3.