On the respect to the hashin-shtrikman bounds of some analytical methods applying to porous media for estimating elastic moduli

district, Hanoi, Vietnam Article history: Received 30/03/2021, Revised 12/04/2021, Accepted 19/04/2021 Abstract In this work, some popular analytic formulas such as Maxwell (MA), Mori-Tanaka approximation (MTA), and a recent method, named the Polarization approximation (PA) will be applied to estimate the elastic moduli for some porous media. These approximations are simple and robust but can be lack reliability in many cases. The Hashin-Shtrikman (H-S) bounds do not supply an exact value but a range that has been admitted by researchers in material science. Meanwhile, the effective properties by unit cell method using the finite element method (FEM) are considered accurate. Different shapes of void inclusions in two or three dimensions are employed to investigate. Results generated by H-S bounds and FEM will be utilized as references. The comparison suggests that the method constructed from the minimum energy principle PA can give a better estimation in some cases. The discussion gives out some remarks which are helpful for the evaluation of effective elastic moduli. Keywords: Maxwell approximation; polarization approximation; Mori-Tanaka approximation; effective elastic moduli; porous medium. https://doi.org/10.31814/stce.nuce2021-15(2)-02 © 2021 National University of Civil Engineering 1. Introduction Most realistic materials, natural or man-made, such as rock, concrete, 3D printing materials con- tain several phases including pores inside their micro-structures. Modern technologies allow the description of unit cell materials in detail which facilitates extremely convincing results of the ef- fective properties by using computational homogenization [1–8]. In practical engineering, an "in- stant"estimation, which does not depend too much on resources, is expected. As the distribution of material is random, it is supposed that the possible effective moduli vary in a range. This promoted analytical methods [9–13], which have developed formulas to construct the upper and lower bounds for this effective coefficient. Unfortunately, these formulas may give a large range of effective values, especially in the case of the high contrast between the properties of the matrix and the inclusion. The effective medium approximations (EMA) [14–17] have developed to avoids this drawback, such as ∗Corresponding author. E-mail address: nhunth@nuce.edu.vn (Nguyen, N.) 14 Nguyen, N., et al. / Journal of Science and Technology in Civil Engineering the self-consistent, differential, correlation approximation [18–22], the Maxwell approximation (MA) [23, 24], the Mori-Tanaka approximations (MTA) [25], and the recent Polarization approximation (PA) [26]. These approximations are applicable for only limited types of inclusions. To overcome this drawback, the equivalent-inclusion approach using artificial neural network has been proposed in [27]. Several works studying PAs have clarified the advantages of this method applying for composite materials. This work will study the application of MA,MTA, and PA to compute effective elastic mod- uli of some porous microstructures. In the next section, we will review the Maxwell, Mori-Tanaka, and Polarization approximations. After that, numerical examples will be presented to compare the results of MA, MTA, and PA with Hashin-Shtrikman bounds (H-S bounds) and the finite element method (FEM). Finally, some discussion will be presented in the last section. 2. Briefly review of MA, MTA and PA predicting the effective elastic moduli In this section, we briefly review some analytical approximations which have been used in a wide range of composite materials to estimate the elastic moduli. Considering an isotropic multicomponent material in d-dimensional space (d = 2, 3) consisting of n isotropic components. The matrix phase has the volume fraction vIα and the α-inclusion has the volume fraction vIα. The bulk modulus and shear modulus of the matrix are KM and µM, respectively. Those of the α inclusion phases are KIα and µIα. 2.1. Maxwell approximation Maxwell Approximations, also called as Maxwell-Garnett or Clausius Mossotti approximations [23, 24], for predicting effective elastic moduli of 2-phase materials are written as: Keff = ( vI KI + (d − 1)K∗M + vM KM + K∗M )−1 − K∗M (1) where K∗M = KM 2(d − 1)µM d (2) and µeff = ( vI µI + µ∗M + vM µM + µ∗M )−1 − µ∗M (3) where µ∗M = µM d2KM + 2(d + 1)(d − 2)µM 2dKM + 4dµM (4) 2.2. Mori–Tanaka approximation The MTA, derived as an approximate solution to the field equations for the composite to compute the elastic moduli CMTA, has the expression: CMTA = vMCM + n∑ α=2 vIαCIα : D0α vMI + n∑ α=2 vIαD0α (5) 15 Nguyen, N., et al. / Journal of Science and Technology in Civil Engineering where D0α = [ I + Pα : C−1M : (CIα − CM ]−1 (6) In (5)–(6), CIα,CM are elastic moduli of the α-inclusion and the matrix; I is quadratic unit tenso. The Eshelby tensor P in the 2D case is the symmetric depolarization tensor of the ellipsoids from the α-inclusion phase, determined according to [12]: P1111 = KM KM + µM a22 + 2a1.a2 (a1 + a2)2 + µM KM . a2 a1 + a2  (7) P2222 = KM KM + µM a21 + 2a1.a2 (a1 + a2)2 + µM KM . a1 a1 + a2  (8) P1122 = KM KM + µM  a22 (a1 + a2)2 − µM KM . a2 a1 + a2  (9) P2211 = KM KM + µM  a22 (a1 + a2)2 − µM KM . a1 a1 + a2  (10) P1212 = KM KM + µM  a21 + a22 2(a1 + a2)2 + µM 2KM  (11) where α1, α2 are the semi axes of the ellipse. For the 3-D case, the formula of Eshelby tensor is more complicated, we refer to [28] for more details. From (5), the bulk modulus K and the elastic shear modulus µ formula of Mori-Tanaka approxi- mation can be written as: KMTA = vMKM + n∑ α=2 vIαKIαDKα vM + n∑ α=2 vIαDKα (12) and µMTA = vMµM + n∑ α=2 vIαµIαDµα vM + n∑ α=2 vIαDµα (13) DKα,Dµα are functions depending on the inclusion-shape,DKα,Dµα with α-ellipsoid inclusion phases, are specified: DKα = αµ (P1111 − P1122 − P2211 + P2222) + 2 Pˆ (14) Dµα = αK (P1111 + P1122 + P2211 + P2222) + 2 2Pˆ + 1 2 ( 2αµP1212 + 1 ) (15) Pˆ = 2αMαK (P1111P2222 − P2211P1122) + ( αK + αµ ) (P1111 + P2222) + ( αK − αµ ) + 2 (16) αK = KI KM − 1, αµ = µI µM − 1 (17) We list in Table 1 the function DKα,Dµα for several types of inclusion [? ]: 16 Nguyen, N., et al. / Journal of Science and Technology in Civil Engineering Table 1. DKα, Dµα of several types of inclusion Spherical Needle Platelet DK = KM + 43µM KI + 43µM Dµ = µM + µ ∗ µI + µ∗ µ∗ = µM 9kM + 8µM 6kM + 12µM DK = KM + µM + 13µI KI + µM + 13µI Dµ = 1 5  4µMµM + µI + 2µM + γMµI + γI + 2KI + 4 3µM KI + 13µI  γM = µM 3KM + µM 3KM + 7µM DK = KM + 43µM KI + 43µM Dµ = µM + µ ∗ µI + µ∗ µ∗ = µI 9kI + 8µI 6kI + 12µI 2.3. Polarization approximation The effective elastic moduli Ceff (Keff, µeff) of the isotropic composite maybe defined via the minimum energy principle: ε0 : Ceff : ε0 = inf〈ε〉=ε0 I(ε), I(ε) = ∫ V ε : C : εdV (18) for all macroscopic constant strain tensor ε0 where ε is expressed through the displacement field u, written as ε = 1 2 ( ∇u + (∇u)T ) ,C is fourth rank material stiffness tensor and 〈.〉 denotes the average over the volume V or via the minimum complementary energy principle: σ0 : C−1eff : σ 0 = inf 〈σ〉=σ0 ∫ V σ : C−1 : σdV (19) for all macroscopic constant stress tensor σ0 where the trial stress field σ should satisfy ∇ · σ = 0. Avoiding the complicated problem of (18), I(ε) is reformulated using polarization, then, minimiz- ing only principal part of the formula yields the following trial strain field: εi j = ε 0 i j + 3K0 + µ0 µ0 (3K0 + 4µ0) n∑ α=1 pαklψ,i jkl − 1 2µ0 n∑ α=1 ( pαmiϕ α , jm + p α mjϕ α ,im ) (20) where (K0, µ0) are elastic moduli of the reference material, pαi j is the component ij of the polarization field of the second order tensor pα,ϕα,ψα are harmonic and biharmonic potentials, see [26, 29, 30] for more details. There are several ways to determine the free reference parameters K0, µ0. In this paper, the PA uses dilute solution reference, as most of other EMAs use this solution as the starting point. By using (20) as the optimal polarization trial fields, the PA for the macroscopic elasticity of a general isotropic n-component material has the particular form: KPA =  n∑ α=1 vα Kα + K∗ −1 − K∗ (21) µPA =  n∑ α=1 vα µα + µ∗ −1 − µ∗ (22) 17 Nguyen, N., et al. / Journal of Science and Technology in Civil Engineering where K∗ and µ∗ are solutions of the following equations: n∑ α=2 vα (Kα − KM) ( KM + K∗ Kα + K∗ − DKα ) = 0 (23) n∑ α=2 vα (µα − µM) ( µM + µ∗ µα + µ∗ − Dµα ) = 0 (24) Note that using a suitable trial stress tensor to solve problem (18) came to the same results as (21) and (22) interestingly [29]. 3. Numerical examples In this section, we will examine some micro-structures using MA, PA, and MTA in 2D and 3D cases. Several shapes of inclusions will be considered: circle, ellipse (2D) and platelet, spherical, needle (3D). 3.1. 2D porous examples We consider a sample of the size 1 × √3 mm in two cases of porous medium: (i) void circular inclusions (I1) and (ii) void ellipse inclusions (I2). The axis ratio of ellipse inclusion a/b equals 1/2. The distribution of inclusions is shown in Fig. 2 in which the nearest distance between the center of inclusions is 0.5 mm. The bulk modulus KM and the shear modulus µM of the matrix are 1 kN/mm2 and 0.4 kN/mm2, respectively.   * 2 * 0, n M M K K Kv K K DK KD D DD D § ·  ¨ ¸© ¹ ¦ (24)   * 2 * 0. n M Mv DD D PD D D P PP P P P § ·  ¨ ¸© ¹ ¦ (25) Note that using a suitable trial stress tensor to solve problem (18) came to the same results as (21) and (22) interestingly [26]. 3. Numerical examples In this section, we will examine some micro-structures using MA, PA, and MTA in 2D and 3D cases. Several shapes of inclusions will be considered: circle, ellipse (2D) and platelet, spherical, needle (3D). 3.1. 2D porous examples We consider a sample of the size 1x 3 mm in two cases of porous medium: (i) void circular inclusions (I1) and (ii) void ellipse inclusions (I2). The axis ratio of ellipse inclusion a/b equals 1/2. The distribution of inclusions is shown in Figs. 2a, b in which the nearest distance between the center of inclusions is 0.5 mm. The bulk modulus MK and the shear modulus MP of the matrix are 1 kN/mm2 and 0.4 kN/mm2, respectively. (a) I1 (b) I2 Figure 1. Unit cells with void-circular inclusions I1 and void-ellipse inclusions I2 The bulk modulus and the shear modulus estimated by MA, PA, MTA are shown in Figs. 2 and 3 in the comparison with FEM and H-S bounds. We can see that: (i) with void circular inclusion, the results estimated by MA, PA, MTA coincide. Simultaneously, these moduli show a good agreement with the result from the unit-cell method (FEM); (ii) with void ellipse circular inclusions, the MA results coincide with the upper H-S bounds (HSU), which defer lightly from the results by MTA and PA. We note that in FEM implementation, 160000 regular tri-elements have been utilized. The material coefficients of the void phase are extremely small (1E-6). This is easy to (a) I1   * 2 * 0, n M M K K Kv K K DK KD D DD D § ·  ¨ ¸© ¹ ¦ (24)   * 2 * 0. n M Mv DD D PD D D P PP P P P § ·  ¨ ¸© ¹ ¦ (25) Note that using a uitable tri l stress tensor o lve pr bl m (18) came to th same results as (21) and (22) interesti gly [26]. 3. Numerical xamples In this section, we will examine some micro-structures using MA, PA, and MTA in 2D and 3D cases. Several shapes of inclusio s w ll be considered: ci cle, ellipse (2D) and platelet, spherical, needle (3D). 3.1. 2D porous exampl s We consider a sample of th size 1x 3 mm in two cases of porous medium: (i) void circular inc usio s (I1) and (ii) void ellipse inclusio s (I2). The axis ratio of ellipse inclusio a/b equals 1/2. The distribut on of inclusio s is shown in Figs. 2a, b in which the near st distance betw en th center of inclusio s i 0.5 mm The bulk modulus MK and the shear modulus MP of the ma rix are 1 kN/mm2 and 0.4 kN/mm2, respectively. (a) I1 (b) I2 Figure 1. Unit cells with void-circular nclusio s I1 a d void-ellipse inclusio s I2 The bulk modulus and the shear modulus estimat d by MA, PA, MT are shown in Figs. 2 and 3 in the comparison w th FEM and H-S bounds. We can see that: (i) with void circular inclusio , the results timat d by MA, PA, MTA coincide. Simultaneously, the e moduli show a g od a reement with the result f om he uni -cell method (FEM); (ii) with void ellipse circular inclusio s, the MA results coincide with the upper H-S bounds (HSU), which defer lightly from the results by MTA and PA. We note that in FEM implementation, 160000 regular tri-elements have been utilized. The material coefficients of the void phase are extrem ly s all (1E-6). This is easy to (b) I2 Figure 1. Unit cells with void-circular inclusions I1 and void-ellipse inclusions I2 The bulk modulus and the shear modulus estimated byMA, PA,MTA are shown in Figs. 2 and 3 in the comparison with FEM and H-S bounds. We can see that: (i) with voi circular inclusion, the results estimated by MA, PA, MTA coincide. Simultaneously, these moduli show a good agreement with the result from the unit-cell method (FEM); (ii) with void ellipse circular inclusions, the MA results coincide with the upper H-S bounds (HSU), which defer lightly from the results by MTA and PA. We note that in FEM implementation, 160000 regular tri-elements have been utilized. The material coefficients of the void phase are extremely small (1E-6). This is easy to learn when comparing the formula of these estimations. However, the agreement reduces remarkably between the results of analytic methods and FEM when the volum fraction of the void phase incr ases. 8 Nguyen, N., et al. / Journal of Science and Technology in Civil Engineering learn when comparing the formula of these estimations. However, the agreement reduces remarkably between the results of analytic methods and FEM when the volume fraction of the void phase increases. (a) Bulk modulus (b) Shear modulus Figure 2. Comparison of elastic moduli estimated by PA, MTA, MA in the case of 2D void-circular inclusions I1. (a) Buck modulus (b) Shear modulus Figure 3. Comparison of elastic moduli estimated by PA, MTA, MA in the case of 2D void-ellipse inclusions I2. 3.2. 2D 3-component examples This section employs some 3-component porous media to compare results generated by PA, MTA, and the H-S bounds. First, a square sample in 2D of the size 1 x 1 mm2 as shown in Fig. 4 is employed to investigate. (a) Bulk modulus learn when comparing the formula of these estimations. However, the agreement reduces remarkably between the results of analytic methods and FEM when the volume fraction of the void phase increases. ( ) l l s (b) Shear modulus i . i f l stic oduli esti ated by PA, MTA, MA in the case of 2D id-circular inclusions I1. (a) c l s (b) Shear modulus igure 3. o parison of elastic oduli esti ated by PA, MTA, MA in the case of 2D void-ellipse inclusions I2. 3.2. 2 3-co ponent exa ples This section e ploys so e 3-co ponent porous media to compare results generated by PA, TA, and the H-S bounds. First, a square sample in 2D of the size 1 x 1 mm2 as shown in Fig. 4 is employed to investigate. (b) Shear modulus Figure 2. Comparison of elastic moduli estimated by PA, MTA, MA in the case of 2D void-circular inclusions I1 hen co paring the formula of these estimations. Howev r, the agreement s re arkably betwe n the results of analytic methods and FEM when the l e fraction of the void phase increases. (a) Bulk modulus (b) Shear modulus Figure 2. Comparison of elastic moduli estimated by PA, MTA, in the case of 2D void-circular inclusions I1. (a) Buck modulus (b) Shear modulus Figure 3. Comparison of elastic moduli estimated by PA, MTA, MA in the case of 2D void-ellipse inclusions I2. 3.2. 2D 3-component examples This section employs some 3-component porous media to compare results generated by PA, MTA, and the H-S bounds. First, a square sample in 2D of the size 1 x 1 mm2 as shown in Fig. 4 is employed to investigate. (a) Bulk modulus learn when comparing the formula of these estimations. However, the agreement reduces remarkably between the results of analytic methods and FEM when the volume fraction of the void phase increases. (a) Bulk modulus (b) Shear modulus Figure 2. Comparison of elastic moduli esti ated by PA, MT , MA in th case of 2D void-circular inclusions I1. (a) Buck modulus (b) Shear modulus Figure 3. Comparison of elastic moduli estimated by PA, MT , MA in th case of 2D v id-ell pse inclusions I2. 3.2. 2D 3-c mponent examples This section employs some 3-c mponent porous media to compare results generated by PA, MTA, and th H-S bounds. First, a squ re sample in 2D of th size 1 x 1 mm2 as shown in Fig. 4 is employed to investigate. (b) Shear modulus Figure 3. Comparison of elastic moduli estimated by PA, MTA, A in the case of 2D void-ellipse inclusions I2 3.2. 2D 3-component examples This sect on empl ys s e 3-component por us media c mpare results generated by PA, MTA, a d the H-S bounds. First, a square sample in 2D of the size 1 × 1m 2 as shown in Fig. 4 is employed to investigate. The elastic properties of components are: (i) th matrix (KM, µM) = (1, 0.4) kN/mm2, (ii) void ellipse inclusions (KI1, µI1) = (0, 0) kN/mm2, (iii) circular inclu ions (KI2, µI2) = (20, 12) kN /mm2. The ellipse inclusions (in dark) have the volume fraction of 10.1% while that of circular inclusions (in white) varies. The ratio between radius of ellipse a/b = 1/5. In FEM implementation, we used 180000 triangular elements, and the material properties of the voids phase are nearly zeros as they are in 2D 2-component examples. 19 Nguyen, N., et al. / Journal of Science and Technology in Civil Engineering Figure 4. A 2D 3-component unit cell with ellipse and circle inclusions The elastic properties of components are: (i) the matrix ( MK , MP ) = (1, 0.4) kN/mm2, (ii) void ellipse inclusions ( I1K , I1P ) = (0, 0) kN/mm2, (iii) circular inclusions ( I2K , I2P ) = (20, 12) kN /mm2. The ellipse inclusions (in dark) have the volume fraction of 10.1% while that of circular inclusions (in white) varies. The ratio between radius of ellipse a/b = 1/5. In FEM implementation, we used 180000 triangular elements, and the material properties of the voids phase are nearly zeros as they are in 2D 2-component examples. Fig. 5 plots the bulk modulus (a) and shear modulus (b). In this case, there is no significant discrepancy between MTA and PA estimation while the discrepancy between those and FEM increases in proportion to the volume fraction of inclusions. (a) Buck modulus (b) Shear modulus Figure 5. Comparison of the results estimated by PA, MTA, FEM of a 2D 3- component unit cell. The square sample consists of three phases: the matrix ( MK , MP ) = (40, 20) kN/mm2, ellipse inclusions ( I1K , I1P ) = (0, 0) kN/mm2, circular inclusions ( I2K , I2P ) = (10, 0.4) kN/mm2. We consider another 2D 3-component sample in which the properties of the matrix are large than those of the inclusions. The properties of the matrix ( MK , MP ), Figure 4. A 2D 3-component unit cell with ellipse and circle inclusions Fig. 5 plots the bulk modulus (a) and shear modulus (b). In this case, there is no significant dis- crepancy between MTA and PA estimation while the discrepancy between those and FEM increases in proportion to the volume fraction of inclusions. Figure 4. A 2D 3-component unit cell with ellipse and circle inclusions The elastic properties of components are: (i) the matrix ( MK , MP ) = (1, 0.4) kN/mm2, (ii) void ellipse inclusions ( I1K , I1P ) = (0, 0) kN/mm2, (iii) circular inclusions ( I2K , I2P ) = (20, 12) kN /mm2. The ellipse inclusions (in dark) have the volume fraction of 10.1% while that of circular inclusions (in white) varies. The ratio between radius of ellipse a/b = 1/5. In FEM implementation, we used 180000 triangular elements, and the material properties of the voids phase are nearly zeros as they are in 2D 2-component examples. Fig. 5 plots the bulk modulus (a) and shear modulus (b). In this case, there is no significant discrepancy between MTA and PA estimation while the discrepancy between those and FEM increases in proportion to the volume fraction of inclusions. (a) Buck modulus (b) Shear modulus Figure 5. Comparison of the results estimated by PA, MTA, FEM of a 2D 3- component unit cell. The square sample consists of three phases: the matrix ( MK , MP ) = (40, 20) kN/mm2, ellipse inclusions ( I1K , I1P ) = (0, 0) kN/mm2, circular inclusions ( I2K , I2P ) = (10, 0.4) kN/mm2. We consider another 2D 3-component sample in which the properties of the matrix are large than those of the inclusions. The properties of the matrix ( MK , MP ), (a) Bulk modulus Figure 4. A 2D 3-component unit cell with ellipse and circle inclusions The elastic properties of components are: (i) the matrix ( MK , MP ) = (1, 0.4) kN/mm2, (ii) void ellipse inclusions ( I1K , I1P ) = (0, 0) kN/mm2, (iii) circular inclusions ( I2K , I2P ) = (20, 12) kN /mm2. The ellipse inclusions (in dark) have the volume fraction of 10.1% while that of circular inclusions (in white) varies. The ratio between radius of ellipse a/b = 1/5. In FEM implementation, we used 180000 triangular elements, and the material properties of the voids phase are nearly zeros as they are in 2D 2-component examples. . 5 plots the bulk modulus ( ) s r lus (b). In this case, there is no signif cant discrepancy between ti ation while the discrepancy between those and FEM increases i r olu e fraction of inclusions. (a) Buck modulus (b) Shear modulus Figure 5. Co parison of the results estimated by PA, MTA, FEM of a 2D 3- component unit cell. The square sample consists of three phases: the matrix ( MK , MP ) = (40, 20) kN/mm2, ellipse inclusions ( I1K , I1P ) = (0, 0) kN/mm2, circular inclusions ( I2K , I2P ) = (10, 0.4) kN/mm2. We consider another 2D 3-component sample in which the properties of the matrix are large than those of the inclusions. The properties of the matrix ( MK , MP ), (b) Shear modulus Figure 5. Comparis n of the results estimated by PA, MTA, FEM of a 2D 3-component unit cell. The square sample consists of three phases: the matrix (KM , µM) = (40, 20) kN/mm2, ellipse inclusions (KI1, µI1) = (0, 0) kN/mm2, circular inclusions (KI2, µI2) = (10, 0.4) kN/mm2 We consider another 2D 3-component sample in which the properties of the matrix are large than those of the inclusions. The properties of the matrix (KM, µM), t e ellipse inclusions (KI1, µI1) a d circular inclusions (KI2, µI2) are (1, 0.4), (0, 0), (20, 12) kN/mm2, respectively. With this set of data, the MTA lightly violates the upper HS bound while those the PA is closely under the upper bounds as shown in Fig. 6. 20 Nguyen, N., et al. / Journal of Science and Technology in Civil Engineering the ellipse inclusions ( I1K , I1P ) and circular inclusions ( I2K , I2P ) are (1, 0.4), (0, 0), (20, 12) kN/mm2, respectively. With this set of data, the MTA lightly violates the upper HS bound while those the PA is closely under the upper bounds as shown in Fig. 6. (a) Buck modulus (b) Bulk modulus (closer look of (a)) Figure 6. Comparison of elastic modulus estimated by PA and MTA of a 2D 3- component unit cell: the matrix ( MK , MP ) = (40, 20) kN/mm2, ellipse inclusions ( I1K , I1P ) =(1, 0.4) kN/mm2, circular inclusions ( I2K , I2P )= (0, 0) kN/mm2. 3.3. 3D 3-component examples In this part, we apply MTA and PA for some 3D porous media with several types of inclusion, including platelet, needle, and sphere. In the following examples, the properties of the matrix ( MK , MP ) are constant at (40, 20) kN/mm2. The bulk modulus is taken into consideration in different cases of volume fraction from low to high. Figs. 7(a-c) plot the estimation of MTA and PA for the case of ellipsoid inclusions ( I1K , I1P ) = (10, 0.4) kN/mm2 and sphere voids ( I2K , I2P ) = (0, 0) kN/mm2. We can observe that PA, MTA, and HSU nearly coincide when the volume fraction of needles is small I1X = 5% and 15%. Whereas, as can be seen in Fig.7c when the volume fraction of needles is 75%, the MTA estimation start to exceed the HSU and PA estimation still respect the upper of H-S bounds. Fig. 8(a-c) plot the estimation of MTA and PA for the case when inclusions are platelets ( I1K , I1P ) = (10, 0.4) kN/mm2 and sphere ( I2K , I2P ) = (0, 0) kN/mm2. In this case, the violation of MTA is first observed in Fig. 8(b) when the platelet phase has a volume fraction of I1X = 15%. This is more obvious in Fig. 8(c) when the sample contains a high proportion of platelet I1X = 75%. Again, the violation to H-S bounds of PA is acknowledged. Similarly, we consider the case when inclusions are platelets and ellipsoids (voids). Figs. 9 (a-c) plot the estimation of MTA and PA in the three cases of platelet inclusion volume fraction 5%, 15%, 75% respectively. The trend is not different from (a) Bulk modulus the ellipse inclusions ( I1K , I1P ) and circular inclusions ( I2K , I2P ) are (1, 0.4), (0, 0), (20, 12) kN/mm2, respectively. With this set of data, the MTA lightly violates the upper HS bound while those t e PA is closely under the upp r bounds as shown in Fig. 6. (a) Buck modulus (b) Bulk modulus (closer look of (a)) Figure 6. Comparison of elastic modulus estimated by PA and MTA of a 2D 3- component unit cell: the matrix ( MK , MP ) = (40, 20) kN/mm2, ellipse inclusions ( I1K , I1P ) =(1, 0.4) kN/mm2, circular inclusions ( I2K , I2P )= (0, 0) kN/mm2. 3.3. 3D 3-component examples In this part, we apply MTA and PA for some 3D porous media with several types of inclusion, including platelet, needle, and sphere. In the following examples, the properties of the matrix ( MK , MP ) are constant at (40, 20) kN/mm2. The bulk modulus is taken into consideration in different cases of volume fraction from low to high. Figs. 7(a-c) plot the estimation of MTA and PA for the case of ellipsoid inclusions ( I1K , I1P ) = (10, 0.4) kN/mm2 and sphere voids ( I2K , I2P ) = (0, 0) kN/mm2. We can observe that PA, MTA, and HSU nearly coincide when the volume fraction of needles is small I1X = 5% and 15%. Whereas, as can be seen in Fig.7c when the volume fraction of needles is 75%, the MTA estimation start to exceed the HSU and PA estimation still respect the upper of H-S bounds. Fig. 8(a-c) plot the estimation of MTA and PA for the case when inclusions are platelets ( I1K , I1P ) = (10, 0.4) kN/mm2 and sphere ( I2K , I2P ) = (0, 0) kN/mm2. In this case, the violation of MTA is first observed in Fig. 8(b) when the platelet phase has a volume fraction of I1X = 15%. This is more obvious in Fig. 8(c) when the sample contains a high proportion of platelet I1X = 75%. Again, the violation to H-S bounds of PA is acknowledged. Similarly, we consider the case when inclusions are platelets and ellipsoids (voids). Figs. 9 (a-c) plot the estimation of MTA and PA in the three cases of platelet inclusion volume fraction 5%, 15%, 75% respectively. The trend is not different from (b) Bulk modulus (closer look of (a)) Figure 6. Comparison of elastic modulus estimated by PA and MTA of a 2D 3-component unit cell: the matrix (KM , µM) = (40, 20) kN/mm2, ellipse inclusions (KI1, µI1) = (1, 0.4) kN/mm2, circular inclusions (KI2, µI2) = (0, 0) kN/mm2 3.3. 3D 3-component examples the case of platelet and sphere inclusion. The MTA may invade but PA always respects the H-S bounds. Note that, in these examples, I2X varies and I I1 I2X X X  . (a) I1 5%X (b) I1 15%X (c) I1 75%X Figure 7. Comparison of Bulk modulus estimated by PA and MTA of a 3D material with the matrix ( MK , MP ) = (40, 20) kN/mm2, ellipsoid inclusions ( I1K , I1P ) = (10, 0.4) kN/mm2 and sphere inclusions ( I2K , I2P ) = (0, 0) kN/mm2, I I1 I2X X X  . (a) I1 5%X (b) I1 15%X (a) υI1 = 5% the case of platelet and sphere inclusion. The MTA may invade but PA lways respects the H-S bounds. Note that, in these examples, I2X varies nd I I1 I2X X X  . ( I (b) I1 15%X (c) I1 75%X Figure 7. Comparis of Bulk modulus estimated by PA and MTA of a 3D material with the matrix ( MK , MP ) = (40, 20) kN/mm2, ellipsoid inclusions ( I1K , I1P ) = (10, 0.4) kN/mm2 and sphere inclusions ( I2K , I2P ) = (0, 0) kN/mm2, I I1 I2X X X  . (a) I1 5%X (b) I1 15%X (b) υI1 = 15% the cas of pl telet and sphere inclusion. The invade but PA always respects the H-S bounds. Note that, in these examples, I2X varies and I I1 I2X X X  . (a) I1 5%X (b) I1 15%X (c) I1 75%X Figure 7. Comparison of Bulk modulus estimated by PA and MTA of a 3D material with the matrix ( MK , MP ) = (40, 20) kN/mm2, ellipsoid inclusions ( I1K , I1P ) = (10, 0.4) kN/mm2 and sphere inclusions ( I2K , I2P ) = (0, 0) kN/mm2, I I1 I2X X X  . (a) I1 5%X (b) I1 15%X (c) υI1 75% Figure 7. Comparison of Bulk modulus estimated by PA and MTA of a 3D material with the matrix (KM , µM) = (40, 20) kN/mm2, ellipsoid inclusions (KI1, µI1) = (10, 0.4) kN/mm2 and sphere inclusions (KI2, µI2) = (0, 0) kN/mm2, υI = υI1 + υI2 1 Nguyen, N., et al. / Journal of Science and Technology in Civil Engineering