Efficient Solution of 3D Solids with Large Numbers of Fluid-Filled Pores Using Eigenstrain BIEs with Iteration Procedure

Donghong Heand Hang Ma

Abstract: To deal with the problems encountered in the large scale numerical simulation of three dimensional (3D) elastic solids with fluid-filled pores, a novel computational model with the corresponding iterative solution procedure is developed, by introducing Eshelby’s idea of eigenstrain and equivalent inclusion into the boundary integral equations (BIE). Moreover, by partitioning all the fluid-filled pores in the computing domain into the near- and the far-field groups according to the distances to the current pore and constructing the local Eshelby matrix over the near-field group, the convergence of iterative procedure is guaranteed so that the problem can be solved effectively and efficiently in the numerical simulation of solids with large numbers of fluid-filled pores.The feasibility and correctness of the proposed computational model are verified in the numerical examples in comparison with the results of the analytical solution in the case of a single spherical fluid-filled pore under uniform pressure in full space and with the results of the subdomain BIE in a number of other cases. The overall mechanical properties of solids are simulated using a representative volume element (RVE) with a single or multiple fluid-filled pores, up to one thousand in number, with the proposed computational model, showing the feasibility and high efficiency of the model. The effect of random distribution of fluid-filled on overall properties is also discussed. Through some examples, it is observed that the effective elastic properties of solids with a large number of fluid-filled pores in random distributions could be studied to some extent by those of solids with regular distributions.

Keywords: Fluid-filled pores, boundary integral equation, eigenstrain, near-filed group,local Eshelby matrix, equivalent inclusion, mechanical property.

1 Introduction

The determination of the elastic states and effective elastic properties in elastic solids containing many inclusions is of considerable interest in the research of solid mechanics and engineering. In previous investigations, researchers paid more attention to the problem of composites with solid inclusions. Some composites such as wet soils, porous rocks, fault gouges as well as biological tissues, however, are often filled with fluids in their pores which can be described as porous skeletons. Due to the effect of the property of fluid and the initial fluid pressure in pores, the mechanical behaviors of such materials containing fluid inclusions are different with those with solid inclusions. These materials with fluid-filled pores are not only existing widely in nature but also applied broadly to the geophysics, hydrogeology, petroleum science, petrophysics, etc. Especially, during exploiting energy in the deep earth, the prediction of the elastic properties of deep-seated rocks, which usually contain a large number of fluid-filled pores, is of great significance for energy exploitation.

Therefore, there has been a growing interest in studying the various properties of solids with fluid-filled pores in recent decades. Robert et al. [Robert, Conoir and Franklin(2006)] considered the propagation of elastic waves through two-dimensional lattices of water-filled inclusions in an aluminum matrix. Markov et al. [Markov, Kazatchenko,Mousatov et al. (2010)] studied the permeability of the fluid-filled inclusions in porous media. Markov et al. [Markov, Kazatchenko, Mousatov et al. (2012)] analyzed the dielectric properties of fluid-saturated bone. A variety of different techniques also has been developed and used by investigators to interpret and characterize the effective mechanical properties for such materials. The effective moduli and related problems of solids with many cracks and cavities were addressed by Kachanov et al. [Kachanov(1993); Kachanov, Tsukrov and Shafiro (1994)]. Subsequently, Kachanov et al.[Kachanov, Tsukrov and Shafiro (1995)] analyzed the fluid pressure polarization and effective elastic response to materials with fluid-saturated cracks and cavities. In particular, Shafiro et al. [Shafiro and Kachanov (1997)] also extend this method to materials with fluid-filled pores of various shapes. In addition, Levin and Markov [Levin and Markov (2005)] presented the elastic characteristics estimation for rocks containing isolated saturated pores without fluid overflow. James [James (2013)] applied the poroelasticity theory to granular composites having orthotropic porous grains and fluidfilled pores. Style et al. [Style, Boltyanskiy, Allen et al. (2015)] showed experimentally and theoretically that Eshelby’s theory breaks down for small liquid inclusions in a soft solid. Song et al. [Song, Hua and Rudnicki (2016)] studied the shear properties of heterogeneous fluid-filled porous media with spherical inclusions. Ma et al. [Ma and Yang (2018)] applied the micromechanics theory to determine the local elastic fields of the porous media with fluid-filled pores.

Although numerous researches have been done for solids with fluid-filled pores, however,most of them are analytical solutions under some specific assumptions, which are only suitable for analysis of solids with few fluid-filled pores. As for the analysis of solids with multiple solid inclusions and cracks, numerical methods including the finite element methods (FEM), a well-established numerical method, and boundary element methods(BEM) have been used [Shen and Yi (2001); Qin and Mai (2002); Dong and Lee (2005);Kakavas and Kontoni (2006)]. An advantage of the BEM exists obviously over the FEM,in that the discretization is required only on surfaces such as boundaries, interfaces and cracks in the BEM instead of discretizing the whole domain in the FEM for elastic problems containing cracks or inclusions. However, for large-scale problems of inhomogeneity using the BEM, special techniques of fast multipole expansions have to be employed [Greengard and Rokhlin (1987); Liu, Nishimura, Tanahashi et al. (2005)],which leads to complexity of the algorithm. The limitation in the use of the FEM and the BEM, however, would be the efficiency in dealing with the large scale problems, coming primarily from the huge quantity of meshes in the FEM and the full system matrix in the BEM when the number of inclusions or cracks increases since the interface unknowns needs to be solved simultaneously. Recently, by introducing Eshelby’s idea of the eigenstrain and equivalent inclusion into the BEM, Ma et al. proposed the eigenstrain formulation of the BIE for numerical analysis of solids with multiple particles [Ma, Yan and Qin (2009); Ma, Fang and Qin (2011)] and with multiple cracks [Ma, Guo,Dhanasekar et al. (2013)] and with both multiple particles and cracks [Tang, Ma and Yan(2016)]. In this way, the difficulty encountered in the large scale numerical simulation of solids with solid inclusions and cracks was solved effectively. However, for studies of solids with large numbers of fluid-filled pores, it is still a challenge task, considering the cases of random distributions in position, shape and property of fluid-filled pores as well as the interactions among them, especially for the three dimensional solids. Zhang et al.[Zhang, Lv and Zheng (2010); Wang and Henann (2016)] applied the finite element method to solids with liquid inclusions. The limitation of these methods is that the solution scale would be large since both the matrix and every fluid-filled pore should be discretized, which is not suitable for solving the large scale numerical simulation. In addition, a boundary element method was developed by Huang et al. [Huang, Zheng and Yao (2011)] for solving the problems of 2D solids with fluid-filled pores. This method, as practice confirms, is inefficient in solving solids with large numbers fluid-filled pores,especially for the three dimensional solids. Most contributions in the literatures are concerned with 2D solids and with few fluid-filled pores. The simulation of 3D solids with large numbers of fluid-filled pores has been seldom reported.

In the present work, a novel computational model of the eigenstrain BIE with iteration procedure is developed further and extended to the numerical simulation of the threedimensional solids with large numbers of fluid-filled pores. Furthermore, with the aid of the discrete form of eigenstrain BIE in full space, the local Eshelby matrix are suggested and constructed on the group of near field fluid-filled pores defined by distances to the current fluid-filled pore. In this way the strong interaction among fluid-filled pores can be taken into full consideration in the algorithm so that the convergence of iteration can be guaranteed in the numerical simulation of solids with large numbers fluid-filled pores using the eigenstrain formulation of the BIE. As the unknowns appear only on the boundary of the solution domain, the solution scale of the solids with multiple fluid-filled pores with the present model can be greatly reduced, a significant feature because such a traditionally time-consuming problem with multiple fluid-filled pores can be solved efficiently compared with the existing numerical models of the FEM or the BEM.

The paper is organized as follows: The computational model for solids with fluid-filled pores is introduced in details in Section 2, including the eigenstrain boundary integral equations, the Eshelby tensor and the replacement of equivalent inclusions. The Eshelby matrix is defined and constructed after grouping of the fluid-filled pores. The solution procedures of iteration are also described. The numerical examples are presented in Section 3, in which computed results are compared with the analytical solution in the case of a single spherical fluid-filled pore in full space and with those of the subdomain BIE in other cases. The overall mechanical properties of solids are computed using a representative volume element (RVE) with a single or multiple fluid-filled pores distributed either regularly or randomly with the proposed computational model. The efficiency is compared for the two computational models, the proposed model and the subdomain BIE. The conclusions are summarized in the last section.

2 Computational models

2.1 Formulation of eigenstrain boundary integral equations

In the present model, a 3D solid containing NIpores filled with fluids is considered as shown in Fig. 1, where the solid is assumed to be isotropic material denoted as Ω with the outer boundary Γ, representing the solution domain of elastic medium or the base matrix.The fluid filled in pores is assumed to be linear nonviscous and compressible. The domain of the filled-pores is denoted by ΩI(I=1,2,…, NI) with the boundary ΓI(ΓI=ΩI∩Ω)or the interface. For the solution domain Ω considered, the displacement and stress fields of the problem can be described by the eigenstrain formulations of the BIE respectively as follows [Ma, Yan and Qin (2009)]:

where

stands for the free term resulted from the domain integral of Eq. (2). p and q in Eqs. (1),(2) and (3) denote the source and field point, respectively.ui,σijand τistand for the displacements, stresses and tractions, respectively.represents the eigenstrain in ΩIafter the fluid-filled pores being replaced by the equivalent inclusions. Ωεin Eq. (3) is an infinitesimal region in ΩIwith the boundary Γεaround the source point p located inside ΩIand xl= xl( q ) -xl(p). In Eqs. (1) and (2),,andstand for the Kelvin’s fundamental solutions for displacements, tractions and stresses, respectively.,andare correspondingly the derived fundamental solutions. In the boundary integral equations above, the eigenstrains or the states of pores filled with fluids are all unknowns to be determined step by step in the solution procedures. It needs to be pointed out, however, that these eigenstrain BIEs can only describe, as a matter of fact, the displacement and stress fields of homogeneous linear elastic media. In order to obtain the solutions of solids with fluid-filled pores via the eigenstrain BIEs, the replacement of pores by equivalent inclusions has to be carried out beforehand with the aid of Eshelby tensors, which will be described in what follows.

Figure 1: Model of 3D solid with multiple fluid-filled pores

2.2 Eshelby tensor and equivalent inclusion

The following relation holds true according to Eshelby’s original work [Eshelby (1957);Eshelby (1959)].

where ν and μ are the Poisson’s ratio and shear modulus of the elastic media, respectively.δijis the Kronecker symbol. It needs to be pointed out that Eq. (5) is exact only for the uniform distribution of eigenstrain in ΩI.

It is noticed that the fluid can bear only pressure and correspondingly the compressive‘stresses’ and ‘strains’ with minus sign. As mentioned previously, the fluid in the pore is assumed to be linear nonviscous and compressible, therefore the constitutive relation for fluids confined in the pores in the present work can be expressed as:

where the deviatoric components of the applied and the constrained strains as well as the eigenstrains are defined asrespectively. Eq. (8) reflects the fact that there is no shear component either for stress or for strain in the pore filled with fluid. Eqs. (7) and (8) can be rewritten in explicit form as

where

By substituting (4) into (9), the eigenstrain in the equivalent inclusion can be obtained via the applied strains for the single fluid-filled pore using the following equations

which can be written concisely in a matrix form as

It is seen from the property of static fluid that although the scalar component of strains,, of fluids in Eq. (7) are uniformly constant and position independent within each fluid-filled pore. However, the constrained strains so as to the eigenstrains are, in general case, position-dependent within the pore because of the geometry dependence of Eshelby tensor. This is also true for the deviatoric components. As the pores with spherical and ellipsoidal shapes are considered in the present work, the uniform distribution of constant eigenstrain is true for a single fluid-filled pore in full space under the uniform far-field loading, which suggests that the constrained strainsare uniformly constant in Eq. (9)for a single fluid-filled pore.

2.3 Grouping of pores and local Eshelby matrix

For solids with multiple fluid-filled pores, it is obvious that the applied strains or the applied stress at each fluid-filled pore will be disturbed by other fluid-filled pores. That is,the interactions exist among pores, especially for those in the local zone surrounding the current pore in concern. The gratitude of the interactions depends on the distances among fluid-filled pores, coming obviously from the self-evident physical effect among the pores and being seen readily from the distance dependence of fundamental solutions in the BIE formulations. In this regard, when computing the applied strains for the current pore via applied stresses, in addition to the far-field loading, all of the other pores surrounding the current pore as well as those in the domain of concern should be taken into consideration. More precisely, the interactions give rise to two effects to be taken into consideration adequately. The first effect lies apparently in that the applied strains owing to loading are no longer constant within the pores, so as to the eigenstrains.However, the constant assumption of eigenstrain would be accurate enough provided that the distance conditions are satisfied [Tang, Ma and Yan (2016)], which has been verified previously in the computational practices [Ma, Yan and Qin (2009); Ma, Fang and Qin(2011)]. In the present work, the eigenstrain and the applied strain are computed at the geometrical center of each pore. The second effect is considered to be the principal factor that the convergence of iteration procedures will be interfered by the interactions. The gratitude of the interactions depends on the distances of pores to the current pore in concern. The shorter the distance it is, the stronger the interactions among them, and vice versa. In the present work, in order to deal with the strong interactions of fluid-filled pores with short distances in the local region on the current pore in consideration, all the pores in the computing domain are divided into two groups as shown in Fig. 2. The nearfield group for the current pore I in full space is defined as those placed within the circle of dashed line centered at the current pore, neglecting tentatively the effect of the far-field group of pores, consisting of those placed outside the circle. With such definitions, the near-field pores belong to the short distance group with relatively strong effect of interactions on the current pore while the far-field pores belong to the long distance group with relatively weak effect of interactions on the current pore.

Figure 2: Definition of the near-field group for the current pore I in full space

Suppose that the number of pores in the near-field group is N in the full space and neglect tentatively the pores in the far-field group, the stresses in the pores of the near-field group can be expressed as follows, using the stress Eq. (2) of the eigenstrain BIE:

where the following integral identity [Ma, Fang and Qin (2011)] has been employed:

After interfaces being discretized and following the derivation procedures of pore replacement by equivalent inclusions with the aid of Eshelby tensor similar to those carried out and described in Subsection 2.2, a matrix form of the discretized eigenstrain BIE for the near-field group can be derived as

where the submatrix of each diagonal termSkkin [S] reflects the relation between the eigenstrain and the applied strain for each fluid-filled pore itself, having similar structure to that ofSin Eq. (12). The explicit expression ofSplease refers to Eq. (11). In contrast,the submatrix of each off diagonal termSjk(j≠k) in [S] correlates the relations between the eigenstrain and the applied strain for different fluid-filled pores in the near-field group,reflecting the mutual effects of interactions among pores. It needs to be pointed out that every near-field group for each fluid-filled pore is generally unique in the solution domain, especially in the case of random distributions of pores. Suppose the current fluid-filled pore is numbered by I in global series, the eigenstrain vectorfor each fluid-filled pore can be computed by the applied strains vector,{ε} ={ε1,ε2,… ,εN}Tof the fluid-filled pores in its own near-field group by the following matrix expression

where TIis obtained from [S] by inversion then contraction of it, named as the local Eshelby matrix in the present work. From the derivation process of the local Eshelby matrix obtained from the BIE with the aid of Eshelby’s idea for the near-field group of pores irrespective of those in the far-field group, it can be seen that the computed eigenstrains in the current pore can be considered accurate in the numerical sense because the BIE itself is an analytical expression. The errors will come partly from the discretization and partly from the assumption of the constant eigenstrain inside fluidfilled pores so that the strong interactions can be coped with thoroughly and successfully among the fluid-filled pores in the near-field group.

Figure 3: Flow chart of the solution procedure

2.4 Solution procedure of iteration

After replacing all the fluid-filled pores by the equivalent inclusions following Eshelby’s idea, the solution domain becomes formally homogeneous with eigenstrains distributed in it, which can be solved by the BIE. Eq. (1) can be rewritten in matrix form as

where A, B represent the system matrix and the coefficient matrix formed from the domain integrals Eq. (1), respectively. b stands for the right vector formed by known boundary quantities with the corresponding integral coefficients, andxdenotesthe boundary unknowns.ε0represents the total eigenstrain vector to be solved. It is noticed thatA,Bandbare coefficients which just need to be computed once. At the initiation stage, the local Eshelby matrices are computed firstly for each pore one by one, and the eigenstrain vector ε0is assigned an initial value by neglecting tentatively the effect of the fluid-filled pores (i.e., the whole solution domain only contains the homogeneous matrix without pores). Then the boundary unknownsxcan be computed by the following iterative formula:

where k stands for the iteration count. The maximum differences of the eigenstrain components between two consecutive iterations are defined as

The convergence criterion is checked at each iteration step as follows

where E denotes the Young’s modulus of the solids. If the criterion does not meet, the boundary unknowns are renewed using Eq. (19). The flow chart of the iterative procedures is shown in Fig. 3. The crucial steps can be summarized as follows:

(a) Prepare the coefficientsA,B,bin Eq. (18) and the local Eshelby matrix for all the near-field groups.

(b) Compute the outer boundary unknowns via Eq. (18) without eigenstrains.

(c) Assign an initial value for eigenstrain vector ε0.

(d) Renew the boundary unknowns.

(e) Compute the eigenstrain with applied strains in the near-field and the far-field groups.

(f) Check the convergence criterion (21).

(g) If the convergence criterion (21) is not met, return to Step (e). Otherwise, go to next Step (h).

(h) Compute the stresses in solid and the pressures of fluids in pores as well as the overall properties of solid.

In the present work, when the distances between the source and the field points, p and q,are relatively small, the following boundary type integrals are used instead of the domain integrals [Ma, Fang and Qin (2011)], which corresponds to the case for the near-field group

respectively for the domain integrals in Eqs. (1) and (2). Notice the difference between Eqs. (23) and (14) because the domain integral in (2) or (23) is strong singular. The free term O*ijklexists or not depending on if the field and source points, p and q, are coincide in ΩIor not. Whereas when the distances between the source and field points are relatively large, the one-point computing [Ma, Fang and Qin (2011)] is employed to approximate these domain integrals as follows:

where VIrepresents the volume of the pore, which corresponds to the case for the farfield group.

Figure 4: (a) The axes of ellipsoid; (b) Meshes of the outer surface; (c) Meshes of the interface of a pore in one octant

3 Numerical examples

In this section, some numerical examples are presented to demonstrate the accuracy and efficiency of the eigenstrain formulation of BIE and the corresponding iterative solution procedure for the 3D solids with fluid-filled pores. As very few analytical solutions exist,all the results are obtained at least by two procedures for the purpose of comparison, the proposed eigenstrain BIE computational model in the present paper and the subdomain BIE proposed by Huang et al. [Huang, Zheng and Yao (2011)], where the solids with fluid-filled pores were analyzed using 2D model. For the pores with spherical shape, the radius of the pore is denoted as r0, while for the pores with elliptical shape, the definition of the axes of the ellipsoid is shown in Fig. 4(a), where a, b and c are three half radii of the ellipsoidal pore. In the following analysis, a=b is used and the aspect ratio is defined as the ratio of c and a (i.e. the shape of the elliptical pore is chosen as either oblate(a=b>c) or prolate (a=b

Figure 5: A thick-wall spherical shell containing compressible fluid under uniform outer pressure

3.1 Verification of the proposed method

In order to verify the correctness and accuracy of the proposed computational model, a simple example of a spherical fluid-filled pore in full space is presented first. The analytical solution of it can be obtained from a thick-wall spherical shell filled with compressible fluid subjected to a uniform pressure, p0, on its outer boundary as shown in Fig. 5. The inner pressure in the pore, p, is induced by p0. As mentioned at the beginning of Section 2, the pore with linear compressible fluid is assumed in the present model. The analytical solution of the radical stress and displacement can be written as follows for the shell

where a and b denote the inner and outer radius of the spherical shell, respectively. E and ν represent the elastic modulus and Poisson’s ratio of the shell. Similar to the deriving process in the literature [Huang, Zheng and Yao (2011)], the inner pressure, p, as shown in Fig. 5, can be obtained as

where K is the compressibility of the fluid. When the outer radius, b, approaches to infinite, Eq. (28) becomes

which is the analytical solution of the inner pressure with a spherical fluid-filled pore in full space, under a uniform pressure.

The problem of a single spherical fluid-filled pore in full space is computed by the two numerical procedures, the subdomain BIE and the eigenstrain BIE. The numerical results of the inner pressures, the radial displacements and the radial stresses of solid media at point(x1=2r0, x2=0, x3=0) with different compressibility of the fluid are listed in Tabs. 1, 2 and 3 in comparison with the analytical solutions. The minus values in the table indicate compression. It can be seen that the results of both the two numerical procedures are in good agreement with those of the analytical ones. Moreover, the results of eigenstrain BIE perform even better than that of the subdomain BIE, showing the correctness and feasibility of the proposed computational model. From Tab. 1 through 3, it also can be observed that whether for the inner pressures or the radial displacements and stresses of solid media, the results will decrease with the increasing of the compressibility of the fluid.

In the case of elliptical fluid-filled pores and all the following examples, as there is no analytical solution found, the comparison is carried out between the two numerical procedures, the subdomain BIE and the eigenstrain BIE. The inner pressures in a single elliptical fluid-filled pore as a function of aspect ratio are compared in Figs. 6(a) and 6(b)under different far-filed loading modes and different compressibility of fluid, respectively.In the computation, the fluid-filled pore is loaded with unit load in far-field. It can be seen from Fig. 6 that the results of both the two numerical procedures are in good agreement with each other. Moreover, it also indicates that the inner pressures are greatly affected by the shape of the pores and compressibility of fluid as well as the external loading modes.

Table 1: Comparison of inner pressures in a single spherical fluid-filled pore

Table 2: Comparison of radial displacements at the point (x1=2r0, x2=0, x3=0) in solid media

Figure 6: Inner pressures in a single elliptical pore as a function of aspect ratio, c/a,under different load modes (a) and different compressibility of fluid (b)

In addition, the stresses at solid media adjacent to the single elliptical fluid-filled pore under far-field single compression in x3direction is also compared and presented in Fig. 7. Further,the stresses at solid media adjacent to two and four spherical fluid-filled pores under farfield uniform triaxial compression are compared and presented in Figs. 8(a) and 8(b),respectively. It can be seen from Figs. 7 and 8 that the results of both the two numerical procedures, the eigenstrain BIE and the subdomain BIE, are in good agreement with each other. It also can be seen from Fig. 8 that the assumption of constant eigenstrain in pores is appropriate and feasible for at the present computational conditions, if the non-dimensional distances between pores are not too small in the case of multiple pores.

Figure 7: Dimensionless stresses at solid media adjacent to the pore under single compression in x3

Figure 8: Dimensionless stresses at solid media adjacent to the two spherical fluid-filled pores (a) and to the four spherical fluid-filled pores (b) under far-field uniform triaxial compression.

Figure 9: The RVE with a single fluid-filled pore (a) and with triply periodically spaced fluid-filled pores (b)

3.2 The overall properties of the RVE

In this subsection, the effective elastic properties of a cube RVE, the 3D solid with fluidfilled pores as shown in Fig. 9, are solved, containing either a single or triply periodically spaced spherical fluid-filled pores, respectively. The size of the RVE is set as H=W=20. The outer boundary of the RVE and each of the interfaces are discretized by using 164 and 290 nodes, respectively, as shown in Figs. 4(b) and 4(c). In the computation, the RVE is loaded with unit load in either single compression or uniform triaxial compression. The volume fraction is defined as the ratio of the total volume of pores and the volume of the RVE. In the present work, the range of volume fraction varies from 0 to 0.2.

3.2.1 Overall properties of RVE with a single fluid-filled pore

In order to further verify the correctness and feasibility of the proposed algorithm in present work, the effective elastic properties of RVE with a single fluid-filled pore (NI=1) is solved using the proposed computational model and compared with the subdomain BIE. The overall bulk modulus,, of RVE with an ellipsoidal pore is shown in Fig. 10 as a function of volume fraction under uniform triaxial unit compression. The numerical results show that the overall bulk modulus will increase with the increase of volume fraction, and that the lower the aspect ratio of the ellipsoid is, the lower the overall bulk modulus will be.

The various non-dimensional overall properties, including elastic modulus,/ E, Poisson’s ratio,, shear modulus,/E, bulk modulus,, of RVEs with a single fluid-filled pore as a function of aspect ratio, c/a, are shown in Fig. 11 under uniaxial compression or uniform triaxial compression. The anisotropic phenomenon can be observed in the case of elliptical fluid-filled pores. From the intersection point of the dashed lines in Figs. 11(b) and 11(c), it can be observed that if c/a=1 then the fluid-filled pore is spherical, and the overall properties become isotropic, consistent with the physical facts.

Figure 10: Overall bulk modulus of RVE with a single fluid-filled pore as a function of volume fraction

Figure 11: Overall properties of RVE with a single fluid-filled pore as a function of aspect ratio: Overall elastic modulus (a), overall Poisson’s ratio (b), overall shear modulus (c) and overall bulk modulus (d)

It can be seen also from Figs. 10 and 11 that the results of both the two numerical procedures, the eigenstrain and the subdomain, are in good agreement with each other,showing further the correctness and feasibility of the proposed computational model in the present work.

3.2.2 Overall properties of RVE with multiple fluid-filled pores

In this subsection, the situation of regular distribution of multiple pores is considered.The RVE with triply periodically spaced multiple fluid-filled pores is shown schematically in Fig. 9(b). As the degree of freedom increases very fast with the increase of the number of fluid-filled pores, the computing program using the subdomain BIE will not work on the desktop computer (Intel i7-4770 CPU, 3.40 GHz). Therefore, the comparison between the two numerical procedures, the eigenstrain BIE and subdomain BIE, cannot be performed for the RVE containing more fluid-filled pores. However, there is no such limitation for the computational model proposed in the present work.

Figure 12: Overall properties of RVE with multiple fluid-filled pores as a function of pore number, NI: Elastic modulus (a) and Poisson’s ratio (b)

The results for the overall properties of the RVE with triply periodically spaced spheroidal fluid-filled pores as a function of the total fluid-filled pore number, NI, are shown in Fig. 12, using the eigenstrain BIE, while the volume fraction of pores is kept constant. It can be seen that the values of computed overall elastic modulus and Poisson’s ratio become stable gradually when the pore number grows large enough, that is, greater approximately than 216. Therefore, 216 fluid-filled pores in the RVE are employed in the following numerical examples in the present work.

Figure 13: Overall bulk modulus of RVE with multiple fluid-filled pores as a function of fluid compressibility, K

Figure 14: Overall elastic modulus (a) and Poisson’s ratio (b) of RVE with multiple fluid-filled pores as a function of volume fraction

The overall properties of the RVE with multiple spherical pores under uniaxial compression are presented in Figs. 14(a) and 14(b), respectively, as a function of volume fraction. It is seen from Fig. 14() that the overall elastic modulus decrease monotonically as expected with the increase of fluid volume fraction. The higher the compressibility of fluid is, the lower the overall elastic modulus will be. However, the overall Poisson’s ratio behaves more complicated with volume fraction as shown in Fig. 14(b).

The overall Poisson’s ratio and shear modulus of the RVE with multiple pores are presented in Figs. 15(a) and 15(b), respectively, as a function of aspect ratio under uniaxial compression at different fluid compressibility. The isotropic phenomenon is also observed in the case of multiple spherical fluid-filled pores (i.e., a=b, c/a=1) from the intersection point of the two dashed lines in Fig. 15. The values of Poisson’s ratio again behave in a more complex fashion with the increase of the aspect ratio, c/a.

Figure 15: (a) Overall Poisson’s ratio and (b) shear modulus of RVE with multiple fluidfilled pores as a function of aspect ratio, c/a

Figure 16: Schematics of the RVE with multiple fluid-filled pores in random distribution by position, size and angle (a) and the definition of tilting angle, θ, of elliptical pore (b)

3.3. Effect of random distribution on overall properties

It is more often in reality that the fluid-filled pores are randomly distributed in solids by positions, shapes, sizes and properties, etc. In the present work, the pseudo random function in computer program written in Fortran language is employed to generate random distributions. The random distributions of pores by position and size as well as angle is shown schematically in Fig. 16(a) in the RVE with multiple fluid-filled pores, and the tilting angle, θ,of elliptical pore is shown in Fig. 16(b), where the angle, β, is kept constant. The definitions of the near-field group are similar to those correspondingly in regular distributions except that the number of pores may be different in the near-field group in random distributions.

The non-dimensional overall properties, including elastic modulus,/ E , Poisson’s ratio,, shear modulus,/E, of RVEs with multiple fluid-filled pores as a function of volume fraction, are shown in Figs. 17(a) and 17(b) for both the regular and random distributions,respectively, under uniaxial compression, and the random distribution is consisted of the random variation of positions of elliptical pores combined with the random variation of tilting angles, θ, of ellipsoid as shown in Fig. 16(b). It can be seen from Fig. 17 that both the overall elastic modulus and shear modulus decrease with the increase of volume fraction as expected, and the differences of the various overall properties are negligibly small between the regular and random distributions.

The bulk modulus of the RVE with two kind of randomly distributed fluid-filled pores are also computed and compared with those of solids with regularly distributed fluid-filled pores.The first kind of distribution is composed of the random variation of compressibility of pores in the range K=10-9-100combined with the random variation of positions of elliptical pores.The second kind includes, in addition to those of compressibility and position being the same with that of the first kind, the random variation of tilting angles, θ, of ellipsoid as shown in Fig. 16(b) as well as the sizes of pores with a fixed aspect ratio, c/a, of 0.7. The computed results, the overall bulk modulus of the RVE with multiple pores as a function of volume fraction under uniformly triaxial compression, are summarized in Fig. 18, showing that the overall bulk modulus of the two random distributions correspond to those of the regular distributions with the compressibility K=4.95×10-2and K=2.25×10-1, respectively,which implies that when the number of randomly distributed fluid-filled pores becomes large enough in the solid media of finite size, the effect of fortuitous cancellation can be observed for the fluid-filled pores in random distributions. The effective elastic properties of solids with the fluid-filled pores in random distributions could be studied to some extent by those of solids with regular distributions. However, as the correlation between the two seems of intricacy, further researches need to be carried out carefully and systematically for this issue in future work.

Figure 17: Comparison of overall elastic modulus (a) and overall Poisson’s ratio and shear modulus (b) of RVE with multiple pores between random and regular distributions as a function of volume fraction

Figure 18: Comparison of overall bulk modulus of RVE with multiple pores between random and regular distributions as a function of volume fraction

3.4 Comparison of efficiencies

As mentioned previously, when the number of fluid-filled pores becomes larger, the program of the subdomain BIE method will not work in the present desktop computer(Intel i7-4770CPU, 3.40 GHz). Therefore, using a few fluid-filled pores, the degree of freedom and CPU time of two computational models, the eigenstrain BIE and subdomain BIE, are compared in Tab. 4 as well as in Fig. 19 under uniform triaxial compression,showing that the efficiency of the proposed eigenstrain computational model is much higher than that of the subdomain procedure when the total number of pores grows larger.

Table 4: Comparison of the degrees of freedom and CPU times for the two algorithms

Figure 19: Comparison of CPU times for the two computing procedures, the eigenstrain and the subdomain, as a function of total pore numbers, NI

This is because the size of the system matrix of the subdomain procedure will increase with the total pore number, while the size of the system matrix of the proposed eigenstrain BIE model remain unchanged as the unknowns appear only on the outer boundary of solution domain, while the interface unknowns do not appear in the system equations. It is known that the implementation of the conventional BEM results in a full system matrix. The cost to solve such a matrix will grow exponentially with the increase of its size so that special techniques of fast multipole expansions [Greengard and Rokhlin(1987); Liu, Nishimura, Tanahashi et al. (2005)] have to be employed for large scale problems. In the proposed computational model, as a matter of fact, the solution of local Eshelby matrices does pay cost as the number of local Eshelby matrices, which needs to be solved also, will increase with the growing of pore number. However, from the computational point of view, the workload of solving a problem with many small matrices, corresponding to the case of the proposed model, would be much efficient than that of solving a single huge matrix, which is just the case of the subdomain procedure.More the number of pores involved, much greater the difference of efficiencies between the two approaches will be, as shown in Tab. 4 and Fig. 19, especially for the 3D problem.In the present work, it is observed from the computational practice that the convergence can usually be achieved by 2-6 iterations. The maximum iteration time is no more than 9 so that the convergence can be guaranteed, indicating that the strong interactions among the pores in the near-field group are overcome by introducing the corresponding local Eshelby matrix.

4 Conclusions

In this paper, a novel computational model of eigenstrain BIE with the corresponding iterative solution procedure is developed for analysis of 3D solids with fluid-filled pores in great numbers. All the fluid-filled pores in the computing domain are divided into the near-field groups and far-field groups according to the distances to the current pore. The concept of local Eshelby matrix, defined on pores in the near-field group having strong effects of interaction on the current pore, has been introduced into the computational model to resolve the problem of interactions among fluid-filled pores. Moreover, the construction of local Eshelby matrix is considered as a key step of the present method to guarantee the convergence of iteration. The feasibility and correctness of the proposed computational model are verified in comparison with the results of the analytical solution in the case of a single spherical fluid-filled pore in full space and of the subdomain BIE in all other cases. The results for the overall mechanical properties are simulated and presented using a cube RVE with single or multiple fluid-filled pores, up to one thousand in number, with the proposed computational model. As the unknowns appear only on the boundary of the solution domain, the solution scale of solids with multiple fluid-filled pores with the present model can be remained fairly small, a significant feature because such a traditionally time-consuming problem with multiple fluid-filled pores can be solved efficiently compared with the existing numerical models of the FEM or the BEM.In addition, the effect of random distribution of fluid-filled on overall properties is also discussed in the present work.

Acknowledgements:The research work has received funding from the National Natural Science Foundation of China (Grant Nos. 11672173, 11272195).