Numerical homogenization of elastic properties and plastic yield stress of rock-like materials with voids and inclusions at same scale

Journal Pre-proof Numerical homogenization of elastic properties and plastic yield stress of rock-like materials with voids and inclusions at same scale Y.J. Cao, W.Q. Shen, J.F. Shao, W. Wang PII:

S0997-7538(19)30663-1

DOI:

https://doi.org/10.1016/j.euromechsol.2020.103958

Reference:

EJMSOL 103958

To appear in:

European Journal of Mechanics / A Solids

Received Date: 15 August 2019 Revised Date:

24 December 2019

Accepted Date: 23 January 2020

Please cite this article as: Cao, Y.J., Shen, W.Q., Shao, J.F., Wang, W., Numerical homogenization of elastic properties and plastic yield stress of rock-like materials with voids and inclusions at same scale, European Journal of Mechanics / A Solids (2020), doi: https://doi.org/10.1016/ j.euromechsol.2020.103958. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Masson SAS.

Numerical homogenization of elastic properties and plastic yield stress of rock-like materials with voids and inclusions at same scale Y.J. Cao1 , W.Q. Shen2∗ , J.F. Shao1,2∗ , W. Wang1 1

Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering, HOHAI University, Nanjing 210098, China 2 Univ. Lille, CNRS, Centrale Lille, FRE 2016 - LaMcube - Laboratoire de m´ecanique multiphysique et multi´echelle, F-59000, Lille, France *Corresponding authors: [email protected], [email protected]

Abstract Mechanical properties of rock-like materials are affected by voids and mineral inclusions. Different kinds of analytical solutions have been established for the determination of effective elastic properties and plastic yield stress or failure strength of those materials. However, voids and inclusions are generally considered at different scales. By adopting the Fast Fourier Transform technique, we propose in this paper a numerical homogenization method for modeling the effective elastic properties and plastic yield stress of rock-like materials having voids and mineral inclusions at the same scale. Both are embedded in an elastic-perfectly plastic matrix. The overall elastic-plastic properties are estimated by taking into account the morphology and spatial distribution of pores and inclusions. It is found that the macroscopic elastic modulus is sensitive to both void and inclusion geometry and orientation. The overall yield stress is mainly controlled by the porosity, the inclusion volumetric fraction, the void shape, distribution and orientation. Due to the presence of voids, the influences of inclusion geometry and orientation on the overall yield stress are significantly reduced. This is an important difference with inclusion-reinforced materials without voids. Keywords: Numerical homogenization, FFT based method, Voids, Inclusions, Anisotropy, Yield stress, Rock-like materials 1. Introduction Heterogeneous geomaterials are widely encountered in many engineering fields. These materials contain different kinds of heterogeneities such as voids and mineral inclusions. Their mechanical and physical properties are generally affected by the presence of those heterogeneities. On the other hand, the porosity and mineral compositions may significantly Preprint submitted to Elsevier

January 31, 2020

vary in space, for instance with the depth. They can also be modified by different physicalchemical processes. Therefore, it is desired to establish multiscale constitutive modelings which are competent to consider the effect of heterogeneities and their evolution. Significant advance has been obtained about the estimation of macroscopic elastic property of heterogeneous materials with different kinds of micro-structures. Different kinds of estimates and bounds have been obtained either using Eshelby solution based homogenization techniques or variational principles [50, 49]. Some multi-step homogenization techniques have recently been proposed for porous materials [44]. [14] has studied the macroscopic elastic property of rock-like geomaterials. Concerning plastic behavior, a great deal of studies have also been achieved during the last decades. We mention here the most famous model proposed by Gurson [18] related to determine an explicit expression of macroscopic yield function for porous materials having an incompressible matrix which is described by a von-Mises function and spherical void. After that, the Gurson’s criterion has been widely extended, for example, by considering the effect of void shape (spheroidal voids) [15, 16, 33], the plastic anisotropy of solid phase [2, 32, 22], etc. Effective plastic yield functions for porous geomaterials having a pressuresensitive solid phase and spherical voids have also been developed in [17, 28]. The effect of void shape has further been investigated by [42]. Some investigations have been carried out on the effective strength criteria of rock-like materials reinforced by rigid inclusions [1]. In the works of [45, 48, 40], the effective mechanical properties of double porous materials have been investigated. The plastic behavior of porous materials having voids and inclusions at two different scales have also been investigated[12, 41]. In parallel with these analytical and semi-analytical approaches, significant efforts have been made on the direct full-field simulations of heterogenous materials, for instance by using the Finite Element Method [19, 46, 9, 10], Fast Fourier Transform based techniques [47, 4, 5] and other numerical solutions [20, 11, 29, 30]. These numerical models are able to take the influence of size and distribution of voids or mineral grains in considerations [3, 27, 23, 8]. Numerical simulations have also been used to verify some analytical models [13, 7]. However, in most all previous studies, voids and inclusions are generally taken into account at two separate scales. In many geomaterials like concrete and clayey rocks, mineral inclusions and voids are distributed at the same scale. Further investigations are then needed necessary for this type microstructure. In present work, a new numerical homogenization method will be proposed for studying macroscopic elasto-plastic responses of geomaterials containing voids and inclusions at the same scale. The Fast Fourier Transform technique will be adopted in this micro-macro model. By means of a battery of numerical simulations, we 2

shall investigate respectively the effect of porosity, inclusion volumetric fraction, geometry and orientation of oblate and prolate voids and inclusions on macroscopic elastic modulus and plastic yield stress. 2. Local and macroscopic behavior of studied microstructure During the last decades, significant advances have been obtained on different imaging techniques of micro-structures of heterogeneous materials. 2D SEM and 3D X-ray image technique are among the most widely used techniques. They have also been used for rock-like materials [25, 21, 26, 39]. Some specific studies have been performed on the micro-structure of claystone [38] and shale rocks [37]. According to these studies, the typical structure at the mesoscopic scale of such materials can be characterized as: mineral grains and voids are immersed in a solid matrix. As a first approximation, voids and inclusions can be assumed to be spherical or spheroidal. Based on these observations, we have chosen a simplified micro-structure representation for this class of materials, as shown in Fig.1-a. A periodic 3D micro-structure is further adopted. The simplified 3D unit cell is represented by a cube with one centered void and 1/8 inclusion in each corner of the cube, as described in Fig.1-b. The corresponding porosity and inclusion volumetric fraction are respectively denoted as f and ρ. Inclusion Pore Matrix

(a) 2D structure of studied materials

(b) Selected 3D representative volume element

Figure 1: The representative microstructure of unit cell with voids and inclusion at the same scale The mechanical property of the solid matrix at the mesoscale is assumed to be isotropic and elasto-plastic while the one of inclusion is elastic. The pressure sensitivity of the matrix is taken into account by employing the following Drucker-Prager model: F(σ) = σd + α(σm − σ0 ) ≤ 0

(1)

in which σ is the local stress whose mean part is σm and deviatoric tensor is σ0 . σd is √ calculated by σd = σ0 : σ0 . α presents the local frictional coefficient while σ0 is the 3

hydrostatic tensile strength of the solid matrix. With the incremental form, the local elasticplastic constitutive behaviors of the solid matrix is given by: ∆σ(x p ) = Ctan s (x p ) : ∆ε(x p )

(2)

Ctan s denotes the tangent elastic-plastic stiffness operator of the solid matrix. According to the plastic consistency condition, its expression can be given as:   s  (F < 0 or F = 0, F˙ < 0)  C  tan ∂F ∂F (3) Cs =  ) ⊗ (C s : ∂σ ) (C s : ∂σ    ˙ C s − (F = 0, F = 0) H with:

∂F ∂F : Cs : (4) ∂σ ∂σ C s assigns as the elastic stiffness of the matrix. While the elastic stiffness of inclusions is determined by the tensor Ci . Owing to the existence of inclusions and voids, the stress and strain fields in the matrix are also heterogenous. To minimize the effect of unit cell size on the macroscopic properties, a periodic kinematic boundary condition with a macroscopic strain tensor E is applied on the studied unit cell. Therefore, the kinematic effect of voids and inclusions can be described by a periodic displacement fluctuation u∗ (x p ). Accordingly, the non-uniform strain field is expressed as the sum of the average part E and the fluctuation one ε(u∗ (x p )). Thus the determination of effective response of heterogenous materials comes to the resolution of the following local boundary values problem:    ∆σ(x p ) = Ctan ∀x p ∈ Ω s  s (x p ) : ∆ε(x p)      i   ∀x p ∈ Ωi ∆σ(x p ) = C (x p ) : ∆ε(x p ) (5)    ∗  divσ(x ) = 0 ∀x ∈ Ω, u #, σ · n − #  p p      ε(x p ) = ε(u∗ (x p )) + E ∀x p ∈ Ω H=

where # and −# are respectively the periodic and anti-periodic boundary condition. Ω denotes the whole region of the studied unit cell; Ω s and Ωi the regions for the solid matrix and inclusions. σ and ε respectively denote the local stress and strain in Ω. To investigate the macroscopic mechanical behaviors of highly heterogenous geomaterials and in particular to study the effects of voids and inclusions geometry and distribution, the numerical approach including the Fast Fourier Transform(FFT) technique, initially developed by [35], is here adopted. The key technique is the reformulation of periodic LippmannSchwinger equation which is solved by using a regular discrete grid into Fourier space. This 4

method avoids complex meshing process and is able to consider various geometrical forms of voids and inclusions [34]. Its accuracy has been demonstrated through the comparison with the Finite element method in various situations [31], [24]. It has been successfully applied to rock-like materials [4]. The description of this method for nonlinear materials is available in a number of previous papers, for instance [34] and [4]. Here the FFT-based numerical implement has presented in algorithm 1 in detail to compute the macroscopic mechanical properties for the studied heterogeneous material considering pores and inclusions. In this paper, the emphasis is put on the numerical study of the influence of voids and inclusions on macroscopic elastic and plastic behaviors.

5

Algorithm 1: Macroscopic homogenization algorithm Input: ε(tn , x p ), ∆E(tn+1 ), ∆tn+1 Output: E(tn+1 ), Σ(tn+1 ) Initialization:tn+1 = tn + ∆tn+1 ; E(tn+1 ) = E(tn ) + ∆E(tn+1 ); ε0 (tn+1 , x p ) = ε(tn , x p ) + ∆E(tn+1 ) ∀x p ∈ Ω; if x p ∈ Ω s then Using implicit radial return algorithm to compute local stress σ0 (tn+1 , x p) of the solid matrix; else if x p ∈ Ωi then σ0 (tn+1 , x p ) = Ci (x p ) : ε(x p); else σ0 (tn+1 , x p ) = 0; end end for i = 0 : Niter do The previous σ(tn ) and ε(tn ) at each point x p are known ; σ ˆ i (tn+1 , ξ p ) = F F T (σi (tn+1 , x p)); i 2 i)1/2 Convergence test Eerror = (hkξ·σkˆσˆ(ξ)k ; i (0)k −4 if error < 10 then Return; else ˆ i (tn+1 , ξ p ) εˆ i+1 (tn+1 , ξ p ) = εˆ i (tn+1 , ξ p ) − Γˆ 0 (ξ p ) : σ ∀ξ p , 0, εˆ i+1 (0) = E(tn+1 ); −1 i+1 i+1 ε (tn+1 , x p) = F F T (εˆ (tn+1 , ξ p )); if x p ∈ Ω s then Using implicit radial return algorithm to compute local stress σi (tn+1 , x p) of the solid matrix; else if x p ∈ Ωi then σ0 (tn+1 , x p ) = Ci (x p ) : ε(x p); else σ0 (tn+1 , x p ) = 0; end end i = i + 1; end end R 1 Calculate the macroscopic stress Σ(tn+1 ) = |Ω| σ(tn+1 , x p)dV Ω In this paper, the emphasis is put on the numerical study of the influence of voids and inclusions on macroscopic elastic and plastic behaviors. Analogous to laboratory test, the macroscopic properties including elastic modulus and yield stress in p − q plane can be respectively obtained by conducting a series of uniaxial and triaxial tests. In those tests, both displacement and stress conditions are generally prescribed on the sample boundaries. But the algorithm 1 only takes the prescribed macroscopic strain as the input parameter. 6

Therefore, there is a need to improve this algorithm in order to consider mixed boundary conditions in numerical calculations. As an example, the case of triaxial compression test is here considered. The macroscopic boundary conditions are given by: 12 13 Σ11 0 = Σ0 , E0 = 0, E0 = 0 22 23 Σ21 0 = 0, Σ0 = Σ0 , E0 = 0

(6)

32 33 Σ31 0 = 0, E0 = 0, E0 = E 33

The following loading steps are usually used in this kind of test. The macroscopic hydrostatic stress is first applied to a given value Σ0 in the e1 , e2 and e3 directions. This value corresponds to the confining stress for triaxial test. Then in the direction e3 , a macroscopic strain increment is prescribed, while a constant confining stress is kept in the e1 and e2 directions. To deal with this kind of mixed boundary conditions, in the framework of the Newton-Raphson’s iterative method, a specific algorithm is here used. An additional macroscopic stress balance condition is introduced in order to keep the confining stress at the selected value in the e1 , e2 directions. Further, the penalty term technique is used in order to reproduce the prescribed strain increment in the direction e3 . The modified algorithm is illustrated in 2. Algorithm 2: Newton iteration algorithm Input: ε(tn , x p ), ∆E0 (tn+1 ), Σ0 (tn+1 ); Output: δE(tn+1 ), E(tn+1 ), Σ(tn+1 ) Initialization: E(tn+1 ) = E(tn ) + ∆E(tn+1 ); for k = 0 : Niter do Call algorithm 1 to compute Σ¯ k (tn+1 ); δEk = C0−1 (Σ0 − Σ¯ k (tn+1 )); ∆Ek = δEk + ∆Ek−1 ; if ||δEk ||/||∆Ek || < 10−4 then Return; else k =k+1 end end 3. Estimation of macroscopic elastic modulus The macroscopic elastic modulus is first studied by considering the effect of inclusion and void shape and orientation. Two groups of microteructure having different distributions of heterogeneities (void and inclusion) are studied. In the first group, the unit cell contains a centered spheroidal void (oblate or prolate) and eight 1/8 spherical inclusions at the unit cell 7

corners, as shown in Figure 2(a) - 2(b). Figures 2(c) - 2(d) are for the second group, the unit cell contains a centered spheroidal inclusion (oblate or prolate) and eight 1/8 spherical voids at the corners. For the oblate voids and inclusions, θ denotes the angle from the loading direction (x3) to the minor (symmetric) axis while for the prolate ones, θ is the one from the direction x3 to the major (symmetric) axis. In each case, 150 × 150 × 150 voxels is adopted to discrete the studied representative volume element. The elastic constants of the solid matrix are: E s = 5GPa and v s = 0.15, and those for the inclusions are Ei = 200GPa and vi = 0.15. 3

q

(a) oblate pore

3

3 q

q

(b) prolate pore

(c) oblate inclusion

3 q

(d) prolate inclusion

Figure 2: Studied unit cells with different void and inclusion geometry The first group of unit cells with a centered spheroidal void is first investigated. We consider different orientations θ and aspect ratios c/a of the void (a denotes the semi-major axis while c the semi-minor one of the spheroidal void). For instance, the aspect ratio c/a ranges from 0.2 to 1.0 in the case of oblate void and 1.0 ∼ 2.25 for the prolate one. The porosity and volume fraction of inclusion are respectively selected as f = 0.1 and ρ = 0.1. The value of axial macroscopic modulus, expressed as the relative value with respect to the elastic modulus of solid matrix Ehom /E s , is estimated. For the microtructure with a centered spherical void (c/a = 1), the axial modulus is independent on the angle θ. But for a spheroidal void, the axial modulus varies with the loading angle θ, as described in Figure 3(a) for oblate void and in Figure 3(b) for the prolate one. The results indicate that the macroscopic modulus is strongly dependent on the value of c/a for the oblate voids. The elastic modulus meets its minimum value at θ = 0◦ , and then continuously increases with the orientation angle and obtains its maximum value for θ = 90◦ . This trend agrees well with some experimental data obtained in clayey rocks [36]. That means that for this class of rock-like materials, the strong anisotropy of elastic properties can be attributed to the spacial orientation of oblate voids at a microscopic scale. Compared with the case of spherical voids, the elastic modulus in the parallel direction with the symmetric axis of voids (θ = 0◦ ) is significantly weakened when the value of c/a decreases owing to the 8

increase of compressibility of voids. In the perpendicular direction to symmetric axis, the elastic modulus is only very slightly increased when the aspect ratio decreases. Furthermore, when the value of c/a is small enough, the oblate voids can be assimilated to penny-shaped open cracks. Therefore, the orientation of such cracks strongly affects the macroscopic elastic modulus. However, the macroscopic elastic modulus is much less sensitive to the value of c/a for prolate void than that for the oblate one. The elastic modulus is slightly anisotropic even for a high aspect ratio and the maximum value is obtained when θ = 0◦ . It seems that for case of prolate void, the porosity determines the macroscopic elastic modulus which is little affected by the aspect ratio of voids. 0 1 .4

E

h o m

/E

3 3 0

c /a c /a c /a c /a

3 0

s

1 .2 1 .0

3 0 0

= 1 = 1 /2 = 1 /3 = 1 /5

6 0

0 .8 0 .6 2 7 0

9 0

0 .6 0 .8 2 4 0

1 .0

1 2 0

1 .2 2 1 0

1 .4

1 5 0 1 8 0

(a) oblate pore 0 1 .4

E

h o m

/E

3 3 0

c /a c /a c /a c /a

3 0

s

1 .2 1 .0

3 0 0

= 1 = 3 /2 = 2 = 9 /4

6 0

0 .8 0 .6 2 7 0

9 0

0 .6 0 .8 1 .0

2 4 0

1 2 0

1 .2 1 .4

2 1 0

1 5 0 1 8 0

(b) prolate pore

Figure 3: Evolution of normalized effective elastic modulus related to void orientation: f = 0.1, ρ = 0.1 9

In Figure 4, we show the evolution of macroscopic elastic modulus with the aspect ratio of spheroidal inclusions. As a large difference with voids, the elastic modulus is strongly sensitive to the value of c/a for both two cases, oblate and prolate. The anisotropy of elastic modulus is significantly enhanced when c/a decreases for the oblate inclusion and increases for the prolate one. Compared with the spherical inclusion, the elastic modulus is largely enhanced and reaches its maximum value along the perpendicular direction to the symmetric axis for the oblate inclusion (θ = 90◦ ) but in the parallel direction to the symmetric axis of the prolate one. At the same time, the minimum elastic modulus is smaller than that for the spherical inclusion. But the order of reduction is less important than that of amplification for the maximum value. From these results, it is found that the macroscopic elastic modulus of composites can be strongly dependent on the aspect ratio and orientation of both oblate and prolate inclusions.

10

0 1 .8

E

h o m

3 3 0

/E

3 0

s

1 .6 1 .4

3 0 0

c /a c /a c /a c /a

= 1

c /a c /a c /a c /a

= 1

= 1 /2 = 1 /3 = 1 /5

6 0

1 .2 1 .0 2 7 0

9 0

1 .0 1 .2 2 4 0

1 .4

1 2 0

1 .6 2 1 0

1 .8

1 5 0 1 8 0

(a) oblate inclusion 0 1 .8

E

h o m

/E

3 3 0

3 0

s

1 .6 1 .4

3 0 0

= 3 /2 = 2 = 9 /4

6 0

1 .2 1 .0 2 7 0

9 0

1 .0 1 .2 1 .4

2 4 0

1 2 0

1 .6 1 .8

2 1 0

1 5 0 1 8 0

(b) prolate inclusion

Figure 4: Evolution of normalized effective elastic modulus related to inclusion orientation and aspect ratio: f = 0.1, ρ = 0.1 4. Macroscopic plastic yield stress and comparison In this section, we shall study the influence of inclusion and void on the overall plastic yield stress of heterogenous rock-like materials. For the case of non-spherical voids and inclusions, we will only provide new numerical results due to the absence of analytical macroscopic criteria. The special focus is put on the influence of inclusion and void aspect ratio and orientation on the macroscopic plastic yield stress. But for the particular case of spherical void and inclusion at the same scale, an analytical solution is selected and compared with the numerical results obtained by using the proposed method in this work. 11

In this section, only the results related to the unit cells defined in Figure 2 are presented and discussed. The results for materials having randomly distributed voids and inclusions are presented in Appendix A. For performing numerical simulations, the following mechanical properties of the solid matrix are used: E s = 5GPa, v s = 0.15, α = 0.1 ∼ 0.3, σ0 = 10MPa. The elastic parameters of the inclusion are given as Ei = 200GPa, vi = 0.15. 4.1. Influence of spheroidal inclusion geometries The macroscopic yield stress for a microstructure having a spheroidal inclusion in center and eight 1/8 spherical void in corners is firstly studied, as defined in Figure 2(c) - 2(d). A series of calculations are performed with two different aspect ratios and orientations. The obtained yield surfaces are given in Figure 5. It is interesting to observe that they are very closed and the difference is very small between the spherical and spheroidal inclusions in two different orientations. It is observed that the effects of inclusion geometry and orientation on the overall yield stress are negligible. This result may be due to the fact that the fluctuation of local stress field in the matrix caused by changes of inclusion geometry and orientation is small. As a consequence, the inclusion content ρ is the main factor influencing the overall strength except the hydrostatic case. S p h e re O b la te O b la te P ro la te P ro la te

c /a = c /a = c /a = c /a = c /a =

1 1 /2 1 /5 3 /2

6

Σ33−Σ11(M P a ) 4

2

2 Σm (M P a ) -1 5

-1 2

-9

-6

0

-3

0

3

6

-2 -4 -6

Figure 5: Macroscopic yield stress for unit cells with one centered spheroidal inclusion with different aspect ratios: α = 0.3, f = 0.1, ρ = 0.1, θ = 0◦ .

12

S p h e re O b la te O b la te P ro la te P ro la te

θ=0 θ=9 θ=0 θ=9

6

°

Σ33−Σ11(M P a )

0 ° ° 0 °

4 2 Σm (M P a )

-1 5

-1 2

-9

-6

-3

0 0

3

6

-2 -4 -6

Figure 6: Macroscopic yield stress for unit cells with one centered spheroidal inclusion for different orientation angles with α = 0.3, f = 0.1, ρ = 0.1; oblate case: c/a = 0.5; prolate case: c/a = 2.0. 4.2. Influence of spheroidal void geometries This subsection focuses on the unit cells with a centered oblate and prolate void and eight 1/8 spherical inclusion in corners , as shown in Figures 2(a)-2(b). The spheroidal void morphology influence is investigated through three parameters: porosity, aspect ratio and orientation. In Figure 7(a), we show the macroscopic yield stress of a microtructure with a centered oblate void (c/a = 1/2, θ = 0◦ ) for three different porosities. As for the spherical void, the macroscopic yield stress largely increases when the porosity decreases. But as a basic difference to the spherical void, the plastic yield surface loses the symmetry property with respect to Σm -axis owing to the void geomatry and its orientation. The effects of void aspect ratio on the overall yield surface are further studied by considering four different values. Figure 7(b) shows that the overall material strength is mightily sensitive to the aspect ratio. Especially for the both cases of hydrostatic tension and compression, the yield stress is largely reduced with the decrease of c/a. However, the reduction of yield stress due to the aspect ratio decrease is not uniform for all stress domain. For instance, there exists a stress zone in which the yield stresses are insensitive to the void aspect ratio. For the case of a high stress triaxiality, the yield stress is even decreased by the void aspect ratio decrease. Therefore, the yield surface anisotropy is with the decrease of void aspect ratio c/a. These significant void shape effects are in contrast to the case of inclusion as previously discussed. Finally, the effect of oblate void orientation is investigated. In Figure 7(c), we show the yield stress surface for three different orientations (θ = 0◦ , 45◦ and 90◦ ) and with a porosity of 13

f = 0.1. We notice that the effect of oblate void orientation is significant for the compressive stress domain but relatively small for the tensile one. What’s more, the yield stress under hydrostatic tension seems to be unconstrained by the oblate void orientation θ. The microstructure with one centered prolate void is also studied. In Figure 8(a), the effective yield stress are presented with three porosities. As for the oblate void, the yield strength is mightily affected by the porosity, particularly for the compressive stress domain. The effect of prolate void aspect ratio is studied and illustrated in Figure 8(b). The macroscopic yield stress for three different orientations of the prolate void, θ = 0◦ , θ = 45◦ and θ = 90◦ is given in Figure 8(c). By comparing Figure 8(c) with Figure 7(c), the orientation effect is more important for the prolate void than for the oblate one. Further, it seems that the yield surface has a symmetry property when θ = 0◦ and θ = 90◦ with respect to the axis Σm . This trend also holds in the case of f = 0.05 as presented in Appendix B. We shall finally examine the influence of inclusion content ρ on the overall yield stress. In Figure 9, it is reported the macroscopic yield stress respectively for the unit cells with a centered oblate and prolate void for θ = 0◦ . Three values of inclusion fraction are considered and the case with ρ = 0 corresponds to the unit cell without inclusion. By comparison with the case of unit cell with a centered inclusion, one can see that the presence of inclusion enhances the yield stress except for the hydrostatic tension and compression. Together with the result given in Figure 11 in the following subsection, it is observed that the presence of inclusion does not influence the hydrostatic strength of the material with spherical or spheroidal void. Such a strength is mainly dominated by the strength of matrix, porosity and void properties.

14

8

f= 0 .0 5 f= 0 .0 7 5 f= 0 .1

Σ33−Σ11(M P a ) 6 4 2

-2 0

-1 6

-1 2

-8

0

-4

Σm (M P a ) 0

4

8

-2 -4 -6 -8

(a) Effect of porosity with c/a = 1/2, θ = 0◦ c /a c /a c /a c /a

6

= 1

Σ33−Σ11(M P a )

= 1 /2 = 1 /3 = 1 /5 4 2

-1 2

-1 0

-8

-6

-4

Σm (M P a ) 0

-2

0

2

4

6

-2 -4 -6

(b) Effect of aspect ratio with f = 0.1, θ = 0◦ 6

Σ33−Σ11(M P a )

0 ° 4 5 ° 9 0 °

4 2

-1 2

-9

-6

0

-3

0

Σm (M P a ) 3

6

-2 -4 -6

(c) Effect of void orientation with f = 0.1, c/a = 1/2

Figure 7: Influence of f , c/a and θ of oblate void on the macroscpic yield stress, α = 0.3, ρ = 0.1 15

8

f= 0 .0 5 f= 0 .0 7 5 f= 0 .1

Σ33−Σ11(M P a ) 6 4 2

-2 0

-1 6

-1 2

-8

0

-4

Σm (M P a ) 0

4

8

-2 -4 -6 -8

(a) Effect of porosity with c/a = 2, θ = 0◦ c /a c /a c /a c /a

6

Σ33−Σ11(M P a )

= 1 = 1 .2 5 = 1 .5 = 2

4 2 Σm (M P a )

-1 2

-1 0

-8

-6

-4

0

-2

0

2

4

6

-2 -4 -6

(b) Effect of aspect ratio with f = 0.1, θ = 0◦ 0 ° 4 5 ° 9 0 °

6

Σ33−Σ11(M P a ) 4 2 Σm (M P a )

-1 2

-1 0

-8

-6

-4

0

-2

0

2

4

6

-2 -4 -6

(c) Effect of void orientation with f = 0.1,c/a = 2

Figure 8: Influence of f , c/a and θ of prolate void on the macroscpic yield stress, α = 0.3, ρ = 0.1 16

ρ= 0 . 0 ρ= 0 . 1 ρ= 0 . 2 6

Σ33−Σ11(M P a ) 4 2

-1 2

-1 0

-8

-6

-4

Σm (M P a ) 0

-2

0

2

4

6

-2 -4 -6 (a) oblate c/a = 1/2

ρ= 0 . 0 ρ= 0 . 1 ρ= 0 . 2 6

Σ33−Σ11(M P a ) 4 2

-1 2

-1 0

-8

-6

-4

-2

0 0

Σm (M P a ) 2

4

6

-2 -4 -6 (b) prolate c/a = 2

Figure 9: Macroscopic yield stress for the unit cell having a centered spheroidal void: α = 0.3, f = 0.1 and different inclusion volume fractions ρ 4.3. Comparison with analytical solution for the case of spherical void and inclusion In the following, the macroscopic yield stress for the microstructure having a spherical inclusion(or void) in the center and eight 1/8 spherical voids (or inclusions) on the corners is studied. In addition, the close-form criterion proposed by [43] is selected for the purpose 17

of comparison with numerical results. This criterion was established by using the modified secant method as proposed in [6], and by considering a solid matrix exhibiting anqarbitrary

1 0 second-order yield criterion such as f (σ) = σ ¯ 2d + aσm + bσ2m − c, with σ ¯d = σ : σ0 . 2 For the materials studied in this work, the behavior of matrix is governed by a DruckerPrager criterion. This is a special case of the second-order yield surface. Consequently, the macroscopic yield criterion considering the void and inclusion at the same scale can be given in the following form:

F(Σ, f, ρ) = Θ

Σ2eq

Σ2m 3 2 3 9 2 Σm f − α (1 − ρ)] − (1 − f − ρ)[ (1 − f ) − 4ρα2 ]α2 (7) + [ + 3(1 − f − ρ)α 2 2 4 2 σ0 2 σ0 σ0 2

(1−ρ+ 6−6α f )(1− f − 2 α2 ρ)

3 9−4α2 in which Θ = ; Σeq is the macroscopic equivalent stress, with Σeq = 2 1− f + 9−4α2 ρ 6−6α q 3 Σ : Σ, and Σm denotes the macroscopic mean stress. With this criterion in hand, the 2 influence of porosity and inclusion fraction on the macroscopic yield stress can be addressed. Moreover, the accuracy of this criterion can be also assessed by the numerical homogenization method proposed in this work. In Figure 10, one presents the comparisons of the macroscopic yield stress respectively predicted by the proposed FFT-based method and the closed-form criterion for different values of porosity f . Here the volume fraction of inclusion is fixed as ρ = 0.1. It is clear that the porosity f weakens significantly the overall yield stress of the heterogenous materials in the whole stress domain. However, due to the tension-compression asymmetry of the solid matrix, the influence of porosity is much less sensitive in the tensile domain than in the compressive one. This feature is well captured by both the FFT-based solution and the analytical criterion Eq.(7). The predicted macroscopic yield stresses are quite close to each other at the tensile domain. However, the macroscopic yield stress obtained by the criterion Eq.(7) is overestimated in the compressive domain. It is a common shortcoming of the analytical criteria issued from the modified secant method. The main reason is the fact that a uniform average stress field is assumed for the whole solid matrix instead of considering a highly heterogeneous stress field. This would fail to capture the accurate local strain field. Contrary, the FFT-based solution can well describe such a heterogenous local stress and strain field, so that the obtained results are very closed to those of FEM simulations presented in [31].

18

f = 0 f = 0 f = 0 f = 0

.0 5 .1 0 .1 5 .2 0

-F F -F F -F F -F F

T

f = 0 f = 0 f = 0 f = 0 T T T

.0 5 .1 0 .1 5 .2 0

-C r -C r -C r -C r

ite ite ite ite

rio rio rio rio

n E q n E q n E q n E q

.( 7 .( 7 .( 7 .( 7

-1 2 )

Σ33−Σ11(M P a )

) ) )

-1 0 -8 -6 -4 -2 Σm (M P a )

-2 7

-2 4

-2 1

-1 8

-1 5

-1 2

-9

-6

-3

0 0

3

6

Figure 10: Comparison of macroscopic yield stress between criterion Eq.(7) and FFT-based results in case of spherical void and inclusion for different porosities with α = 0.3, ρ = 0.1 In addition, the influence of inclusion content ρ on the macroscopic yield stress is presented in Figure 11. According to the FFT-based results, one can see that the yield stress under hydrostatic tension and compression is nearly not affected by the inclusion fraction ranging from ρ = 0 to ρ = 0.2. It seems that the plastic yield stress for these two particular stress states is mainly controlled by the porosity and very little influenced by the rigid inclusion. The macroscopic shear yield stress with high triaxiality is significantly enhanced by the inclusion content ρ. Significant differences are obtained with the predictions given by the criterion Eq. (7), in particular for the stress states around the hydrostatic compression path. It seems that the interaction between pores and inclusions at the same scale are not properly taken into account in the selected analytical criterion.

19

ρ= ρ= ρ= ρ=

0 .0 0 .0 0 .1 0 .2

0 -F 5 -F 0 -F 0 -F

ρ= ρ= ρ= ρ=

F T F T F T F T

0 .0 0 .0 0 .1 0 .2

0 -C 5 -C 0 -C 0 -C

rite rite rite rite

rio rio rio rio

n E q n E q n E q n E q

.( 7 .( 7 .( 7 .( 7

-8 )

Σ33−Σ11(M P a )

) ) )

-6

-4

-2 Σm (M P a ) -3 2

-2 8

-2 4

-2 0

-1 6

-1 2

-8

0

-4

0

4

8

Figure 11: Comparison of macroscopic yield stress between criterion Eq.(7) and FFT-based results in case of spherical void and inclusion for different values of ρ with α = 0.3, f = 0.1

α= 0 . 1 - F F T α= 0 . 2 - F F T α= 0 . 3 - F F T

α= 0 . 1 - C r i t e r i o n E q . ( 7 ) α= 0 . 2 - C r i t e r i o n E q . ( 7 ) α= 0 . 3 - C r i t e r i o n E q . ( 7 )

-7

Σ33−Σ11(M P a )

-6 -5 -4 -3 -2 -1 Σm (M P a )

-2 8

-2 4

-2 0

-1 6

-1 2

-8

-4

0 0

4

8

Figure 12: Comparison of macroscopic yield stress between criterion Eq.(7) and FFT-based results in case of spherical void and inclusion for different values of α with f = 0.1, ρ = 0.1 Finally, the influence of friction coefficient α of solid matrix on the macroscopic yield stress is investigated. In Figure 12, one shows the evolution of macroscopic yield stress for different values of α ranging from 0.1 to 0.3. The FFT-based results are compared with the criterion Eq.(7). The macroscopic yield stress is weakened when α is reduced. There is a good agreement of the macroscopic yield stress in the tensile domain between the numerical and analytical solutions. But the macroscopic yield stress in the compression domain is largely overestimated for the criterion Eq.(7), in particular for a high value of α. 20

In view of the comparisons presented above, the macroscopic yield stress of materials with voids and pores at the same scale is not correctly predicted by the analytical criterion. The main reason is the fact that the non-uniform local stress field and the interaction between voids and pores at the same scale are probably not correctly taken into account in the formulation of analytical criteria. The FFT-based numerical homogenization method provides some reference solutions for the evaluation and improvement of analytical criteria. 5. Conclusions A new FFT-based numerical homogenization modeling is established in the present study to investigate the elastoplastic behavior of rock-like geomaterials having voids and inclusions at the same scale. This model is able to investigate not only the effect of porosity or inclusion volume fraction, but also the effect of aspect ratio and orientation of voids and inclusions. This model provides a significant advantage with respect to most analytical and semi-analytical homogenization models. Based on a series of numerical simulations, some interesting remarks can be formulated. The macroscopic elastic modulus becomes dependent on loading orientation due to the presence of spheroidal inclusions and voids. The anisotropy of elastic modulus is strongly enhanced by the decrease of aspect ratio of an oblate void, but only slightly influenced by that of a prolate void. On the other hand, this elastic modulus anisotropy is very sensitive to the aspect ratio of spheroidal inclusion. For geomaterials with the presence of voids and inclusion at the same scale, voids play a dominant role on the macroscopic yield stress or failure strength with respect to inclusions. The macroscopic yield stresses are mightily reduced by the increase of porosity f . They are more sensitive to the value of c/a of oblate void than to the one of prolate void. The macroscopic yield stress becomes strongly dependent on loading direction because of the existence of spheroidal void. The reinforcement of inclusion content ρ for the studied heterogenous geomaterial is only for the zone with high triaxiality. For hydrostatic loading, the overall strength is nearly independent on the inclusion fraction. The yield stress is a slightly affected by the inclusion shape and orientations. In future studies, the real microstructure of rock-like geomaterials will be investigated by experimental method for the numerical study.

21

Acknowledgement This work was partially supported by the National Key RD Program of China (Grant 2017YFC1501101) and the Programme B18019 of Discipline Expertise to Universities MOE & MST. Appendix A. Effect of randomly distributed voids and inclusions on the macroscopic yield stress - Unit cell with spherical voids and inclusions We have also carried out a series of numerical calculations on the unit cell containing randomly distributed voids and inclusions at the same scale (void number N1 = 50, inclusion number N2 = 50). The obtained results are then compared with those for the microstructure with a single centered void or inclusion (as in Figure 2). We have further investigated the effect of frictional coefficient α of the solid matrix, inclusion content ρ and porosity f on the macroscopic yield stress for those two kinds of unit cells. The obtained results are illustrated in Figures A.13 to A.14. One can found that for the same value of porosity and inclusion fraction, the macroscopic yield stress of the microstructure having randomly distributed voids and inclusions is systematically lower than that with a single centered void and 1/8 inclusion at the corners as shown in Figure A.13. According to Figure A.14, this conclusion is also suitable for the case of a single centered inclusion and 1/8 void at the corners. Moreover, this weakening effect of the random distribution of voids or inclusions is more important for the compressive stress zone than for the tensile one. The difference of macroscopic yield stress can probably be attributed to the interaction effect between voids or inclusions for the cases with a random distribution. For the effect of frictional parameter, its increase makes an important reinforcement for the overall yield surface as shown in Figure A.15. Further, it is interesting to observe that the difference of yield stress between the two unit cells is also enhanced by the increase of frictional coefficient.

22

f = 0 f = 0 f = 0 f = 0

.0 5 .1 0 .1 5 .2 0

-o n -o n -o n -o n

e

f = 0 f = 0 f = 0 f = 0 e e e

.0 5 .1 0 .1 5 .2 0

-ra -ra -ra -ra

-7

n d o n d o n d o n d o

m

Σ33−Σ11(M P a )

m

-6 m m

-5 -4

α= 0 . 3 , ρ= 0 . 1

-3 -2 -1

-2 4

-2 0

-1 6

-1 2

-8

0

Σm (M P a )

-4

0

4

8

Figure A.13: Macroscopic yield stress for three values of porosity: comparison between microstructures having randomly distributed voids and a single centered void

-6 ρ= 0 . 0 - o n e ρ= 0 . 1 - o n e ρ= 0 . 2 - o n e

ρ= 0 . 0 - r a n d o m ρ= 0 . 1 - r a n d o m ρ= 0 . 2 - r a n d o m

Σ33−Σ11(M P a )

-5 α= 0 . 3 , f = 0 . 1 -4 -3 -2 -1 Σm (M P a )

-1 5

-1 2

-9

-6

-3

0 0

3

6

Figure A.14: Macroscopic yield stress for different values of ρ: comparison between microstructures having randomly distributed inclusions and a single centered inclusion

23

α= 0 . 3 - o n e α= 0 . 2 - o n e α= 0 . 1 - o n e

-6

α= 0 . 3 - r a n d o m α= 0 . 2 - r a n d o m α= 0 . 1 - r a n d o m

Σ33−Σ11(M P a )

-5 -4 f = 0 . 1 , ρ= 0 . 1 -3 -2 -1 Σm (M P a )

-1 5

-1 2

-9

-6

-3

0 0

3

6

Figure A.15: Influence of frictional coefficient α on macroscopic yield stress: comparison between microstructures having randomly distributed voids and a single centered void - Unit cell with spheroidal voids and spherical inclusions We consider now a microstructure having randomly distributed oblate or prolate voids, all having the same orientation. To avoid the shape effect of inclusions, only randomly distributed spherical inclusions are embedded in the studied microstructure, as presented in Figure A.16.

(a) oblate voids

(b) prolate voids

Figure A.16: Unit cells with randomly distributed spheroidal voids and spherical inclusions These results obtained from the random distribution of oblate voids are given in Figure A.17, for different values of c/a and loading orientation θ. Figure A.18 is for the case of 24

prolate voids. Similarly to the case with a single centered void, the aspect ratio of randomly distributed spheroidal voids has an important influence on the material strength. There is an anisotropy of macroscopic yield stress that is dependent on orientation angle theta. As the case of spherical voids in above subsection, the strength of the microtructure with randomly distributed spheroidal pores is lower than that with one centered spheroidal void for the same aspect ratio, in particular for the compressive stress zone. c /a = 1 c /a = 1 /2 c /a = 1 /5 6

Σ33−Σ11(M P a ) 4 2

-1 8

-1 5

-1 2

-9

-6

Σm (M P a ) 0

-3

0

3

6

-2 -4 -6 (a) θ = 0◦ °

0

4 5 9 0

6

Σ33−Σ11(M P a )

° °

4 2 Σm (M P a )

-1 8

-1 5

-1 2

-9

-6

0

-3

0

3

6

-2 -4 -6

(b) c/a = 1/2

Figure A.17: Effect of aspect ratio and orientation for the microstructure having randomly distributed oblate voids on macroscopic yield stress with α = 0.3, ρ = 0.1 and f = 0.05 25

6

c /a = 1 c /a = 2 c /a = 5 /2

Σ33−Σ11(M P a ) 4 2 Σm (M P a )

-1 8

-1 5

-1 2

-9

-6

0

-3

0

3

6

-2 -4 -6 (a) θ = 0◦ °

0

4 5 9 0

6

Σ33−Σ11(M P a )

° °

4 2 Σm (M P a )

-1 8

-1 5

-1 2

-9

-6

0

-3

0

3

6

-2 -4 -6

(b) c/a = 2

Figure A.18: Effect of aspect ratio and orientation for the microstructure having randomly distributed prolate voids on macroscopic yield stress with α = 0.3, ρ = 0.1 and f = 0.05 Appendix B. Macroscopic yield stress for f = 0.05 We complete here the results presented in Section 4.2 with a porosity of f = 0.05.

26

c /a c /a c /a c /a

= 1

8

= 1 /2 = 1 /3 = 1 /5

6

Σ33−Σ11(M P a )

4 2 -2 4

-2 0

-1 6

-1 2

-8

Σm (M P a ) 0

-4

0

4

8

-2 -4 -6 -8 (a) θ = 0◦

0 ° 4 5 ° 9 0 °

8

Σ33−Σ11(M P a ) 6 4 2

-2 0

-1 6

-1 2

-8

0

-4

0

Σm (M P a ) 4

8

-2 -4 -6 -8 (b) c/a = 1/3

Figure B.19: Effect of aspect ratio and orientation of oblate void on macroscopic yield stress with α = 0.3, ρ = 0.1 and f = 0.05

27

c /a c /a c /a c /a

= 1

8

= 3 /2 = 2 = 3

6

Σ33−Σ11(M P a )

4 2 -2 4

-2 0

-1 6

-1 2

-8

Σm (M P a ) 0

-4

0

4

8

-2 -4 -6 -8 (a) θ = 0◦

0 ° 4 5 ° 9 0 °

8

Σ33−Σ11(M P a ) 6 4 2

-2 0

-1 6

-1 2

-8

-4

0 0

Σm (M P a ) 4

8

-2 -4 -6 -8 (b) c/a = 3

Figure B.20: Effect of aspect ratio and orientation of prolate void on macroscopic yield stress with α = 0.3, ρ = 0.1 and f = 0.05 References [1] Jean-Fran¸cois Barth´el´emy and Luc Dormieux. A micromechanical approach to the strength criterion of drucker-prager materials reinforced by rigid inclusions. International Journal for Numerical and

28

[2] [3]

[4] [5]

[6] [7]

[8]

[9] [10] [11]

[12] [13]

[14] [15]

[16]

[17] [18]

Analytical Methods in Geomechanics, 28(7-8):565–582, 2004. Ahmed Amine Benzerga and Jacques Besson. Plastic potentials for anisotropic porous solids. European Journal of Mechanics-A/Solids, 20(3):397–434, 2001. Nicolas Bilger, Fran¸cois Auslender, Michel Bornert, Jean-Claude Michel, Herv´e Moulinec, Pierre Suquet, and Andr´e Zaoui. Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis. International Journal of Solids and Structures, 42(2):517–538, 2005. Y.J. Cao, W.Q. Shen, N. Burlion, and J.F. Shao. Effects of inclusions and pores on plastic and viscoplastic deformation of rock-like materials. International Journal of Plasticity, 108:107 – 124, 2018. Y.J. Cao, W.Q. Shen, J.F. Shao, and N. Burlion. Influences of micro-pores and meso-pores on elastic and plastic properties of porous materials. European Journal of Mechanics - A/Solids, 72:407 – 423, 2018. Luc Dormieux, Djim´edo Kondo, and Franz-Josef Ulm. Microporomechanics. John Wiley & Sons, 2006. Borys Drach, Igor Tsukrov, and Anton Trofimov. Comparison of full field and single pore approaches to homogenization of linearly elastic materials with pores of regular and irregular shapes. International Journal of Solids and Structures, 96:48–63, 2016. Julie Escoda, Fran¸cois Willot, Dominique Jeulin, Julien Sanahuja, and Charles Toulemonde. Influence of the multiscale distribution of particles on elastic properties of concrete. International Journal of Engineering Science, 98:60–71, 2016. Felix Fritzen, Samuel Forest, Thomas B¨ohlke, Djimedo Kondo, and Toufik Kanit. Computational homogenization of elasto-plastic porous metals. International Journal of Plasticity, 29:102 – 119, 2012. Felix Fritzen, Samuel Forest, Djimedo Kondo, and Thomas B¨ohlke. Computational homogenization of porous materials of green type. Computational Mechanics, 52(1):121–134, 2013. Felix Fritzen and Matthias Leuschner. Reduced basis hybrid computational homogenization based on a mixed incremental formulation. Computer Methods in Applied Mechanics and Engineering, 260:143– 154, 2013. M G˘ ar˘ ajeu and P Suquet. Effective properties of porous ideally plastic or viscoplastic materials containing rigid particles. Journal of the Mechanics and Physics of Solids, 45(6):873–902, 1997. Elias Ghossein and Martin L´evesque. A fully automated numerical tool for a comprehensive validation of homogenization models and its application to spherical particles reinforced composites. International Journal of Solids and Structures, 49(11-12):1387–1398, 2012. A Giraud, NB Nguyen, and D Grgic. Effective poroelastic coefficients of isotropic oolitic rocks with micro and meso porosities. International Journal of Engineering Science, 58:57–77, 2012. Mihai Gologanu, Jean-Baptiste Leblond, and Josette Devaux. Approximate models for ductile metals containing non-spherical voids case of axisymmetric prolate ellipsoidal cavities. Journal of the Mechanics and Physics of Solids, 41(11):1723–1754, 1993. Mihai Gologanu, Jean-Baptiste Leblond, and Josette Devaux. Approximate models for ductile metals containing nonspherical voids case of axisymmetric oblate ellipsoidal cavities. Journal of Engineering Materials and Technology, 116(3):290–297, 1994. TF Guo, Jonas Faleskog, and CF Shih. Continuum modeling of a porous solid with pressure-sensitive dilatant matrix. Journal of the Mechanics and Physics of Solids, 56(6):2188–2212, 2008. Arthur L Gurson. Continuum theory of ductile rupture by void nucleation and growth: Part i–yield criteria and flow rules for porous ductile media. Journal of engineering materials and technology,

29

[19]

[20] [21] [22] [23]

[24] [25] [26]

[27]

[28] [29]

[30]

[31]

[32]

[33]

[34]

99(1):2–15, 1977. J´erˆ ome Julien, Mihail G˘ ar˘ ajeu, and Jean-Claude Michel. A semi-analytical model for the behavior of saturated viscoplastic materials containing two populations of voids of different sizes. International Journal of Solids and Structures, 48(10):1485–1498, 2011. S Kanaun. An efficient homogenization method for composite materials with elasto-plastic components. International Journal of Engineering Science, 57:36–49, 2012. Shaina Kelly, Hesham El-Sobky, Carlos Torres-Verd´ın, and Matthew T Balhoff. Assessing the utility of fib-sem images for shale digital rock physics. Advances in water resources, 95:302–316, 2016. SM Keralavarma and AA Benzerga. A constitutive model for plastically anisotropic solids with nonspherical voids. Journal of the Mechanics and Physics of Solids, 58(6):874–901, 2010. Younis-Khalid Khdir, Toufik Kanit, Fahmi Za¨ıri, and Moussa Na¨ıt-Abdelaziz. Computational homogenization of plastic porous media with two populations of voids. Materials Science and Engineering: A, 597:324–330, 2014. MY Li, YJ Cao, WQ Shen, and JF Shao. A damage model of mechanical behavior of porous materials: Application to sandstone. International Journal of Damage Mechanics, page 1056789516685379, 2016. Laurent Louis, P Baud, and T-F Wong. Characterization of pore-space heterogeneity in sandstone by x-ray computed tomography. Geological Society, London, Special Publications, 284(1):127–146, 2007. Lin Ma, Kevin G Taylor, Patrick J Dowey, Loic Courtois, Ali Gholinia, and Peter D Lee. Multi-scale 3d characterisation of porosity and organic matter in shales with variable toc content and thermal maturity: Examples from the lublin and baltic basins, poland and lithuania. International Journal of Coal Geology, 180:100–112, 2017. Komlanvi Madou and Jean-Baptiste Leblond. Numerical studies of porous ductile materials containing arbitrary ellipsoidal voids–i: Yield surfaces of representative cells. European Journal of MechanicsA/Solids, 42:480–489, 2013. S Maghous, L Dormieux, and JF Barth´el´emy. Micromechanical approach to the strength properties of frictional geomaterials. European Journal of Mechanics-A/Solids, 28(1):179–188, 2009. Qing-Xiang Meng, Huan-Ling Wang, Wei-Ya Xu, and Yu-Long Chen. Numerical homogenization study on the effects of columnar jointed structure on the mechanical properties of rock mass. International Journal of Rock Mechanics and Mining Sciences, 124:104127, 2019. QX Meng, HL Wang, WY Xu, M Cai, Jianrong Xu, and Q Zhang. Multiscale strength reduction method for heterogeneous slope using hierarchical fem/dem modeling. Computers and Geotechnics, 115:103164, 2019. Jean-Claude Michel, Herv´e Moulinec, and P Suquet. Effective properties of composite materials with periodic microstructure: a computational approach. Computer methods in applied mechanics and engineering, 172(1-4):109–143, 1999. Vincent Monchiet, Oana Cazacu, Eric Charkaluk, and Djimedo Kondo. Macroscopic yield criteria for plastic anisotropic materials containing spheroidal voids. International Journal of Plasticity, 24(7):1158–1189, 2008. Vincent Monchiet, Eric Charkaluk, and Djimedo Kondo. Macroscopic yield criteria for ductile materials containing spheroidal voids: an eshelby-like velocity fields approach. Mechanics of Materials, 72:1–18, 2014. H Moulinec and P Suquet. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Computer methods in applied mechanics and engineering,

30

[35]

[36]

[37] [38]

[39]

[40]

[41] [42] [43] [44]

[45]

[46]

[47]

[48]

[49] [50]

157(1):69–94, 1998. Herve Moulinec and Pierre Suquet. A fast numerical method for computing the linear and nonlinear mechanical properties of composites. Comptes rendus de l’Acad´emie des sciences. S´erie II, M´ecanique, physique, chimie, astronomie, 318(11):1417–1423, 1994. H Niandou, JF Shao, JP Henry, and D Fourmaintraux. Laboratory investigation of the mechanical behaviour of tournemire shale. International Journal of Rock Mechanics and Mining Sciences, 34(1):3– 16, 1997. Audrey Ougier-Simonin, Fran¸cois Renard, Claudine Boehm, and Sandrine Vidal-Gilbert. Microfracturing and microporosity in shales. Earth-Science Reviews, 162:198–226, 2016. Jean-Charles Robinet. Min´eralogie, porosit´e et diffusion des solut´es dans l’argilite du Callovo-Oxfordien de Bure (Meuse, Haute-Marne, France) de l’´echelle centim´etrique ` a microm´etrique. PhD thesis, Poitiers, 2008. Tarik Saif, Qingyang Lin, Alan R Butcher, Branko Bijeljic, and Martin J Blunt. Multi-scale multidimensional microstructure imaging of oil shale pyrolysis using x-ray micro-tomography, automated ultra-high resolution sem, maps mineralogy and fib-sem. Applied energy, 202:628–647, 2017. Wanqing Shen, Zheng He, Luc Dormieux, and Djim´edo Kondo. Effective strength of saturated double porous media with a drucker–prager solid phase. International journal for numerical and analytical methods in geomechanics, 38(3):281–296, 2014. WQ Shen, D Kondo, Luc Dormieux, and JF Shao. A closed-form three scale model for ductile rocks with a plastically compressible porous matrix. Mechanics of Materials, 59:73–86, 2013. WQ Shen, J Zhang, JF Shao, and D Kondo. Approximate macroscopic yield criteria for drucker-prager type solids with spheroidal voids. International Journal of Plasticity, 99:221–247, 2017. Roland Traxl and Roman Lackner. Multi-level homogenization of strength properties of hierarchicalorganized matrix–inclusion materials. Mechanics of Materials, 89:98–118, 2015. A Trofimov, S Abaimov, I Akhatov, and I Sevostianov. On the bounds of applicability of two-step homogenization technique for porous materials. International Journal of Engineering Science, 123:117– 126, 2018. Pierre-Guy Vincent, Yann Monerie, and Pierre Suquet. Porous materials with two populations of voids under internal pressure: I. instantaneous constitutive relations. International Journal of Solids and Structures, 46(3-4):480–506, 2009. Pierre-Guy Vincent, Yann Monerie, and Pierre Suquet. Porous materials with two populations of voids under internal pressure: Ii. growth and coalescence of voids. International Journal of Solids and Structures, 46(3):507 – 526, 2009. Pierre-Guy Vincent, Pierre Suquet, Yann Monerie, and Herv?Moulinec. Effective flow surface of porous materials with two populations of voids under internal pressure: Ii. full-field simulations. International Journal of Plasticity, 56:74 – 98, 2014. Pierre-Guy Vincent, Pierre Suquet, Yann Monerie, and Herv´e Moulinec. Effective flow surface of porous materials with two populations of voids under internal pressure: I. a gtn model. International Journal of Plasticity, 56:45–73, 2014. J. R. Willis. The Overall Elastic Response of Composite Materials. Journal of Applied Mechanics, 50(4b):1202–1209, 12 1983. John R Willis. Variational and related methods for the overall properties of composites. In Advances in applied mechanics, volume 21, pages 1–78. Elsevier, 1981.

31

Highlights: -

A new FFT-based homogenization method is developed for materials with voids and inclusions at the same scale; The macroscopic elastic properties and plastic yield stress are determined; The effects of void and inclusion geometry and orientation as well as their interactions are investigated; The proposed method provides some reference results for the evaluation of approximate analytical criteria.

Laboratoire de Mécanique Multiphysique Multiéchelle, CNRS FRE2016 Avenue Paul Langevin – Cité scientifique F – 59655 Villeneuve d’Ascq cedex T. +33 (0)3 20 33 71 52 – F. +33 (0)3 20 33 71 53

Editors, European Journal of Mechanics A/Solids December 24, 2019

Declaration of interest of conflict I am pleased to submit our paper “Numerical homogenization of elastic properties and plastic yield stress of rock-like materials with voids and inclusions at same scale” in «European Journal of Mechanics A/Solids ». This is an original work, which has not been published or submitted elsewhere. We declare that there is no conflict of interest. Jianfu Shao On behalf of all co-authors E-mail: [email protected] Phone: (+33)320337182

Laboratoire de Mécanique Multiphysique Multiéchelle, CNRS FRE2016 Avenue Paul Langevin – Cité scientifique F – 59655 Villeneuve d’Ascq cedex T. +33 (0)3 20 33 71 52 – F. +33 (0)3 20 33 71 53