Effective elastic isotropic moduli of highly filled particulate composites with arbitrarily shaped inhomogeneities

Mechanics of Materials 135 (2019) 35–45

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

Effective elastic isotropic moduli of highly filled particulate composites with arbitrarily shaped inhomogeneities

T

⁎

J.-T. Liua, , F.-Y. Xieb, Q.-C. Hea, S.-L. Tangb, C.-W. Xiaob a b

School of Mechanical Engineering, Southwest Jiaotong University, Chengdu 610031, PR China Institute of Chemical Materials, China Academy of Engineering Physics, Mianyang 621900, PR China

A R T I C LE I N FO

A B S T R A C T

MSC: 74-00 65Z05

The purpose of present work is three-fold: (i) to model thin binding interphases within highly filled particulate composites (HFPCs) as equivalent spring-type interfaces of zero-thickness; (ii) to develop a computational approach for calculating the effective elastic moduli of HFPCs containing arbitrarily shaped inhomogeneities; (iii) to establish a numerical scheme which generates simplified models for substitution of real representative volume elements of HFPCs while ensuring linear isotropic elasticity at the macroscopic scale. At first, a spring-type interface model with explicit physical background is derived to replace the soft interphases in HFPCs, and the macroscopic quantities and effective moduli accounting for the interfacial discontinuities are reformulated. Then, the strong and weak formulations of the boundary value problems of HFPCs are constructed, and a computational approach is developed by using the extended finite element and level set methods so as to determine the effective moduli of HFPCs. A numerical procedure is further elaborated to generate macroscopically isotropic RVEs by dividing a cubic model into differently shaped subdomains, and the minimal number of cutting planes is obtained based on the necessary and sufficient conditions in He (2004). Finally, discussions are made on the inclusion size and shape effects on the overall isotropic moduli of HFPCs.

Keywords: Spring-type interface Soft binding interphase Multi-scale modeling Effective isotropic moduli Highly filled particulate composite

1. Introduction Highly filled composites (HFCs), such as polymer bonded explosives (PBXs) (Kauly et al., 1996; 1997; Banerjee and Adams, 2004; Xiao et al., 2012; He et al., 2018), novel ceramic-polymer composites (Ramakrishna et al., 2001; Bengtsson et al., 2007; Wolff et al., 2014; Babu and de With, 2014; Larson and Zok, 2018), granular media (Moshev and Garishin, 1993) and resin-based restoration composites (Stober et al., 2001; Politi et al., 2018), are materials of both practical and scientific interests. Compared with ordinary matrix-based composites, they consist of phases with high elastic stiffness contrast and very large inclusion fraction (Banerjee and Adams, 2004; He et al., 2018; Ramakrishna et al., 2001; Wolff et al., 2014; Stober et al., 2001; Kaully et al., 2008). The present study is concerned with highly filled particulate composites (HFPCs) exhibiting linear isotropic elasticity at the macroscopic scale and possessing particle volume fraction fp no less than 90%. Within these heterogeneous media, a quite limited amount of binding material exists, but its distributions and mechanical property may essentially affect the failure mechanism, the overall elastic moduli and other macroscopic properties. In this work, we focus on evaluating the effective elastic moduli of these composites, which are the

⁎

preliminary parameters for microstructure design and optimization. For HFPCs, Moshev and Garishin (1993) presented a physical discretization method to develop new structural models and verified its accuracy in evaluating the effective elastic moduli when particle concentration ranges from 0.3 to 0.5. Lemarchand et al. (2002) suggested a nonlinear elastic micromechanics approach to quantitatively estimating the frictional enhancement within HFCs under compressive loads. Banerjee and his coworkers (Banerjee and Adams, 2004; Banerjee et al., 2003) proposed a recursive cell method to compute the effective moduli of highly filled PBXs with fp ranging from 0.10 to 0.92. Heukamp et al. (2005) established a nonlinear homogenization procedure with the aid of a three phase model to account for the effects induced by interface cohesion and transition zones within HFCs. Pichler and Lackner (2009) and Pichler et al. (2012) extended the classical homogenization scheme for upscaling the viscoelastic moduli of HFCs with fp varing from 0.29 to 0.55, and applied the extended scheme to predict the viscous properties of asphalt. More recent works can be found in Refs. Fang and Tian (2018) and GonzalezGutierrez et al. (2018) and the references therein. Note that most of these works consider the interface between two adjacent constituents as being perfect and assume that the particles are of simplified similar shapes and

Corresponding author. E-mail address: [email protected] (J.-T. Liu).

https://doi.org/10.1016/j.mechmat.2019.04.022 Received 30 January 2019; Received in revised form 9 April 2019; Accepted 28 April 2019 Available online 09 May 2019 0167-6636/ © 2019 Elsevier Ltd. All rights reserved.

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

embedded in a regular way. Hence, these models and methods can only to a certain extent properly bridge the connections between microscopic and macroscopic characteristics. As shown by Refs. Kaully et al. (2008), Benveniste (1985), Hashin (1992), Azadi et al. (2009), Gu et al. (2014) and Javili et al. (2017), the macroscopic behavior of particulate composites is closely related to both the size and shape of reinforcement particles, particularly when microcavities, interfacial debondings, particle breakages and other defects are present. Gao et al. (2017) presented a micromechanical model within the framework developed in Ju and Chen (1994) to account for the inter-particle and spring-type interface effects in asphalt mixtures, and derived the closed-form estimates for the dynamic moduli by simplifying all inclusions as spherical particles. Liang et al. (2017) used the phenomenological spring-type interface model as in Qu (1993) and investigated the influence of interfaces on the effective elastic properties of cement pastes. Wang and Luo (2018) applied the spring-type interface model in Qu (1993) to analytically examining the interface effects on the effective moduli of PBXs. However, generation and discretization of a real microstructure or representative volume element of HFPCs remain a tough task due to the extremely high proportion of reinforcement particles and to the random distribution of binding material. To alleviate the difficulty in 3D geometrical modeling and discretization and to accurately capture the stress concentration induced by the high moduli contrast, we shall reformulate the thin soft binding interphases within HFPCs as spring-type interfaces on an explicit physical background and elaborate an efficient computational approach so as to rapidly construct 3D particles of arbitrary shapes and correctly compute the effective elastic moduli of these heterogeneous materials. This paper is organized as follows. The subsequent section is dedicated to recapitulating the modelling of thin binding interphases as spring-type imperfect interfaces and to defining the macroscopic quantities and effective elastic isotropic moduli while accounting the presence of imperfect interfaces. In Section 3, the strong and weak formulations of HFPC boundary value problems are provided, and a computational approach is elaborated to evaluate the effective elastic moduli of HFPCs. In Section 4, numerical examples are given and extensive discussions are made. In particular, a general method is presented to generate representative volume elements (RVEs) for HFPCs with particles of arbitrary shapes so as to ensure linear isotropic elasticity at the macroscopical scale. In Section 5, some concluding remarks are given.

Fig. 2. Replacement of a thin binding interphase with a spring-type interface of null thickness: (a) practical three-phase configuration; (b) two-phase equivalent configuration.

Wolff et al., 2014; Kaully et al., 2008; Banerjee et al., 2003). We are here concerned with the particular but very important case where the inclusion volume fraction fp of HFPCs is equal to or larger than 90%. Thus, the thickness of the binding interphase is much smaller than the characteristic dimensions of the reinforcement particles. The following hypotheses are made about the distribution of the binding interphases: (i) each binding interphase between two adjacent particles is uniform in thickness (Fig. 1(b)); (ii) binding interphases between different particles may have different thickness values. In what follows, we shall first model these binding interphases in HFPCs as spring-type interfaces of null thickness. After that, the effective quantities of HFPCs are defined via a representative volume element (RVE). 2.1. Interfacial relations characterizing the thin binding interphases Consider a generic binding interphase Ω(0) of small thickness h (Fig. 2). The middle surface of Ω(0) and its unit normal vector are ^ (1) symbolized by Γ and n, respectively. Two reinforcement particles Ω 0

^ (2) are supposed to be perfectly bonded to Ω(0) by the interfaces S and Ω 1 ^ (1) and Ω ^ (2) are assumed to be and S . The materials forming Ω(0), Ω 2

individually homogeneous and linearly elastic isotropic (Banerjee and Adams, 2004; Ramakrishna et al., 2001; Banerjee et al., 2003). Thus, their constitutive behavior is governed by the isotropic Hooke’s law:

σ (r ) = (r ) ɛ (r ), (r = 0, 1, 2).

(1)

Above, σ and ε are the second-order stress and strain tensors defined over Ω(r); the strain tensor ε(r) is related to the displacement vector u(r) over Ω(r) by (r)

1 [∇u(r ) + (∇u(r ) )T ]; 2

ɛ (r ) =

2. Modelling of the thin binding interphases and definition of the effective elastic isotropic moduli of HFPCs

(r)

(2)

(r )

the symbol denotes the fourth-order isotropic elastic stiffness tensor of the material forming Ω(r) and is specified by

HFPCs can be considered as homogeneous media at the macroscopic scale (Fig. 1(a)), even though, at the microscopic level, they consist of crystals of arbitrary shapes embedded randomly inside the binder material according to the SEM observations (Banerjee and Adams, 2004;

(r ) = 3κ (r )  + 2μ(r )  with κ

(r)

and μ

(r)

(3)

the bulk and shear moduli. The phases 1 and 2 are the

Fig. 1. Schematic diagram of the HFPC constitution from the macro-scale to the micro-scale: (a) Macroscopic homogeneous media; (b) Microstructure with thin binding interphases; (c) RVE with equivalent spring-type interfaces. 36

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

same, so that the following relations hold:

(•)′ = (•) −

κ (1) = κ (2) = κp, μ(1) = μ(2) = μp .

(4)

1 I ⊗ I, 3

=−

(11)

According to the micromechanics, the foregoing moduli κ* and μ* can be properly determined with the aid of an RVE for a HFPC in question, as schematically shown by the spherical domain in Fig. 1(b). This RVE domain is further replaced by the one in Fig. 1(c) where all interphases are modelled as the spring-type interfaces presented in the last subsection. The macroscopic strain and stress tensors involved in (10) are defined via the RVE Ω in Fig. 1(c) as follows:

In (3),  and  are two complementary orthogonal projectors defined by

=

1 tr(•) I . 3

(5)

where I and  separately represent the second- and fourth-order identity tensors. The binding material forming the interphases is much softer than the one constituting the reinforcement particles, and the thickness h of a generic interphase is very small compared with the characteristic dimension lc of the reinforcement particles. Thus, it is appropriate to model the interphases as spring-type interfaces of null thickness so as to alleviate the difficulty in numerically simulating the complex microstructure in question. The practical three-phase model is shown in Fig. 2(a). Eliminating the binder interphase Ω(0), substituting back a zero-thickness interface geometrically located at the interphase middle ^ (1) and Ω ^ (2) down and surface Γ and extending the materials forming Ω

Σ=

1 2V

∫∂Ω [(σ·nR ) ⊗

E=

1 2V

∫∂Ω [u ⊗

x + x ⊗ (σ ·nR )] dS,

nR + nR ⊗ u] dS,

0

(14)

u = ɛ 0x.

[[u]] = (ωn N + ωT T )· t + 0 (h),

(6)

The elastic energy W stored in Ω is given by

[[t ]] = 0 (h),

(7)

W=

With (14), it follows directly from (13) that Further, due to the fact the microscopic stress field σ is continuous over Ω, the definition (12) of Σ amounts to

Σ=

W=

Above, [[ · ]] and ⟨ · ⟩ are the interfacial jump and average operators defined by [[·]] = ·(2) − ·(1) and · = (·(2) + ·(1) )/2 ; t(r) is the traction vector given by t (r ) = σ (r ) n ; Eb and vb are Young’s modulus and Poisson’s ratio of the binder material and satisfy the relations

Eb =

3κ (0) + μ(0)

,

vb =

3κ (0)

(9)

The interfacial relations (6) and (7) are based on a physically meaningful configuration and derived in a mathematically rigorous way (see Gu et al., 2014 for details). In addition, the parameters ωn and ωT involved have the explicit physical meanings given by (8) and (9). In the literatures (Benveniste, 1985; 2006; Zhu et al., 2011), it is known that the relations (6) and (7) characterize the widely used spring-layer interface model. Up to now, the thin binding interphases in HFPCs have been described by the spring-type interface model defined by (6) and (7). This model will be applied in evaluating the effective moduli of HFPCs with high inclusion volume fractions.

(16)

tr(Σ)tr(E ) 1 Σ′: E′ Σ: E = + 2 6 2

(17)

1 κ *tr(E )tr(E ) + μ*E′: E′. 2

(18)

(a) ɛ10 = ε0 I ,

(19)

⎡0 1 0⎤ (b) ɛ20 = ε0 ⎢ 1 0 0 ⎥ ⎣0 0 0⎦

(20)

where ε0 is a given scalar. Then, the effective moduli κ* and μ* can be calculated by

κ* =

tr(Σ) , 3tr(E )

μ* =

1 E′: Σ′ . 2 E′: E′

(21)

In the next section, the extended finite element method and level set method will be applied to computing κ* and μ*. 3. Numerical computation of the effective isotropic moduli of HFPCs containing arbitrarily-shaped particles HFPCs in question are made of numerous reinforcement crystals of arbitrary shapes distributed randomly inside the binder material. Therefore, numerical approaches become an indispensable tool to predict their elastic response and to evaluate the effective moduli to this end. Among various strategies (Dolbow and Harari, 2009; Harari and Dolbow, 2010; Fries, 2008; Belytschko et al., 2009), the extended finite element method (XFEM) (Belytschko and Black, 1999; Moës et al., 1999) combined with the level set method (LSM) (Osher, 1993) has already been successfully employed to deal with the matrix-based composite problems involving the interfacial traction discontinuity

2.2. Macroscopic quantities and effective isotropic moduli of HFPCs In this study, the macroscopical response of 3D HFPCs under investigation is also assumed to be linearly elastic isotropic, so that it is characterized by (Gu et al., 2014; Duan et al., 2005; 2007a)

Σ = (3κ * + 2μ*) E =κ *tr(E ) I + 2μ*E′.

∫Ω σdV .

To determine κ* and μ*, the following matrices of ε0 in (14) are prescribed successively:

2μ(0)

− . 6κ (0) + 2μ(0)

1 V

where Σ′ and E′ denote the respective deviator parts of Σ and E defined by (11). Introducing (10) into (17), we obtain

N = I − T = n ⊗ n. (8)

9κ (0) μ(0)

(15)

E = ɛ 0.

with

(1 + vb)(1 − 2vb) 1 + vb , ωT = 2h , (1 − vb) Eb Eb

(13)

where V is the volume of Ω and nR denotes the outward unit vector normal to the external boundary ∂Ω of Ω. Recall that the definitions (12) and (13) hold even when the microscopic strain and stress fields over Ω are discontinuous. Next, we prescribe the uniform strain boundary condition on ∂Ω:

up to Γ0, we obtain the two-phase configuration in Fig. 2(b). The newly generated domains are named as Ω(1) and Ω(2), respectively. The equivalent interfacial relations across Γ0 in the two-phase configuration can be derived by employing the same procedure as in Refs. Gu et al. (2014), Bövik (1994) and Benveniste (2006) and taking into account the fact that (0) ≪ (1) = (2) (Banerjee and Adams, 2004; Ramakrishna et al., 2001; Banerjee et al., 2003). Precisely, they are given by

ωn = h

(12)

(10)

Above, κ* and μ* separately represent the effective bulk and shear moduli; Σ and E are the macroscopic stress and strain tensors; (•)′ refers to the deviator part of (•) defined by 37

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

into account the constitutive Eq. (1) and the geometrical relation (2), we finally obtain m

θr

1 1 (k ) ⎞ [[u]](k ) ·⎛ N (k ) + ·[[δ u]](k ) dS T ωT ⎝ ωn ⎠

∫Ω ∇u: : ∇δudV + ∑ ∑ ∫

⎜

r = 1 k = 1 Γr _ k

=

⎟

∫ t ·δudS.

(26)

∂Ωt

Up to now, the weak formulation of the BVP in question has been established. It will be applied in generating a computational approach. Fig. 3. Boundary value problem of a HFPC RVE.

3.2. Extended finite element and level set methods (Benvenuti, 2008; Yvonnet et al., 2008), the interfacial displacement discontinuity (Zhu et al., 2011; Belytschko et al., 2009) or both of them (Farsad et al., 2010; Liu et al., 2015). In what follows, we shall use this strategy to solve the boundary value problems (BVPs) of HFPC RVEs with the uniform strain boundary conditions (14) and to compute the effective elastic isotropic moduli.

We now turn to the elaboration of a computational approach to predict the effective elastic moduli of HFPCs. The RVE shown in Fig. 3 is used to this end. Prior to running into the details, each particle domain Ω(r) has to be described in a computationally efficient mathematical manner. In this work, the level set method (LSM) is adopted. Recall that Ω(r) is a polyhedron and its surface Γr formed by θr flat surface elements is characterized as follows:

3.1. Weak formulation

Γr = {x ∈ 3|fr _ i (x ) = 0, 1 ≤ i ≤ θr }.

Consider an RVE of HFPC, Ω, containing m reinforcement crystals Ω(r) (r = 1, 2, …, m ) as shown in Fig. 3. The subdomain Ω(r) occupied by crystal r is a polyhedron bounded by its surface Γr formed of a finite number of flat surfaces or interfaces Γr _ i such that Γr = ∪i Γr _ i . These spring-type interfaces of zero thickness replace the binder interphases mentioned in the last section and are characterized by the interfacial relations (6) and (7). Further, nr _ i denotes the outward unit vector normal to the flat interface Γr _ i . On the external surface ∂Ω of the RVE Ω, only the Dirichlet condition (14) is prescribed. In this work, body forces are assumed to be negligible. The strong formulation of the foregoing BVP is given by the field Eqs. (1)–(3) and

In the above equation, fr _ i (x ) is the signed distance function defined by Belytschko and Black (1999) and Moës et al. (1999)

fr _ i (x ) = sign[nr _ i ·(x − x *)] min

xsf ∈ Γr _ i

xsf ∈ Γr _ i

Ω(r )

Next, the particle domain Ω

=

∫

n0

u (x ) =

m r = 1 Γr

where ∂Ωt is the part of Ω where the traction vector t is imposed and the effects of the fact that Ω(r) may intersect ∂Ω are assumed to be negligible. Note that, when the uniform strain (14) is prescribed on the whole boundary ∂Ω of Ω, we have ∂Ωt = ⌀ . Invoking the interfacial relation (6) for Γr, which can equivalently be expressed as

1 1 (i) ⎞ = ⎛ N (i) + T ·[[u]](i) , ωT ⎝ ωn ⎠ ⎜

∼ Nℓ1 (x ) Hℓ1 (x ) u͠ ℓ1.

(31)

ℓ1= 1

Hℓ1 (x ) = sign[fr (x )] − sign[fr (x ℓ1)].

(23)

∂Ωt

(30)

Above, uℓ and u͠ ℓ1 represent the standard and enriched nodal displacements, respectively; symbols n0 and n1 separately denote the corre∼ sponding total numbers of nodes; Nℓ(x) and Nℓ1 (x ) are the standard shape functions; the enrichment function Hℓ1 (x ) aims to capture the strong discontinuity in the displacement field and takes the form

(22)

Γr

n1

∑ Nℓ (x ) uℓ + ∑ ℓ=1

t (r )·δu(r ) dS.

∫ σ: δɛdV + ∑ ∫ t·[[δu]] dS = ∫ t ·δudS

(i)

is described by

As implied by the interfacial relations (6) and (7), the displacement field presents a jump across each interface element Γr _ i , while the stress field stays continuous within the whole domain Ω. Thus, the function below is chosen to approximate the displacement field over Ω:

Adding (22) with r = 1, 2, …, m , we get

t

(29) (r)

Ω(r ) = {x ∈ 3|fr _ i (x ) ≤ 0, 1 ≤ i ≤ θr }.

δɛ (r ) dV

Ω

(28)

x * = arg min x − xsf .

which hold over each Ω(r) (r = 1, 2, …, m ), by the interfacial relations (6)-(7) for every Γr (r = 1, 2, …, m ) and by the boundary conditions (14) on ∂Ω. It is now ready to construct the weak formulation of the above BVP. A piecewise continuous virtual displacement vector δu satisfying δu = 0 on ∂Ω is first introduced over Ω. Then, applying the virtual work theorem to each subdomain Ω(r) leads to

∫

x − xsf ,

where sign(•) refers to the signed function (see, e.g., Belytschko and Black, 1999) and x* is given by

divσ (r ) = 0

σ (r ):

(27)

(32)

Substituting (31) into the weak formulation (26) gives the discrete equilibrium equation of the BVP in question:

(KV + KS ) U = F .

(33)

Here, the generalized displacement vector U consists of all nodal displacements uℓ and u͠ ℓ1 to be determined; KV and KS are the stiffness matrices related to the bulk elasticity and to the spring-type interface elasticity. These matrices are calculated by

⎟

(24)

nem

KV =

(23) can be rewritten as

∑ ∫Ωe BT BdV , j=1

m

θr

∫ σ: δɛdV + ∑ ∑ ∫ r = 1 i = 1 Γr _ i

Ω

=

∫ t ·δudS, ∂Ωt

m

1 1 (i) ⎞ [[u]](i) ·⎛ N (i) + T ·[[δ u]](i) dS ω ω n T ⎝ ⎠ ⎜

KS =

⎟

θr

ns

1

∑ ∑ ∑ ∫Γr _i_s B T ⎛ ω ⎜

r=1 i=1 j=1

⎝

n

N (i) +

1 (i) ⎞ T B dS ωT ⎠ ⎟

(34)

with B and B being determined by ∇u = and [[δu]] = B ue . Quantity ue is the nodal displacement vector of an element, nem is the total number of the bulk elements constituting Ω, and ns is the number

Bue

(25)

where θr stands for the number of the flat surfaces forming Γr. Taking 38

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

of the interface elements forming Γr _ i . The external force vector F is computed by

F=

∫ N T ·tdS ∂Ωt

x10 = [0, 0, R 0]T ,

(35)

x ζ0

where N is the shape function matrix. After numerically determining the displacement field by (33), we can calculate the strain and stress fields by the strain-displacement relation (2) and the constitutive Eq. (1) together with (3).

⎡ cos ⎛ 2(ζ1 − 1) π ⎞ sin ⎛ (ζ2 − 0.5) π ⎞⎤ ⎢ ⎥ λ1 λ2 ⎠ ⎝ ⎠⎥ ⎢ ⎝ ⎢ ⎛ 2(ζ1 − 1) π ⎞ ⎛ (ζ2 − 0.5) π ⎞ ⎥ = R 0 ⎢ sin sin ⎥, λ1 λ2 ⎠ ⎝ ⎠⎥ ⎢ ⎝ ⎢ ⎥ (ζ − 0.5) π ⎞ ⎢ ⎥ cos ⎛ 2 λ2 ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ ⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

xλ0 = [0, 0, −R 0]T ,

(36)

λ2 ⩾ 2; ζ1 = 1, 2, …, λ1 with λ1 ≥ 3 and and ζ2 = 1, 2, …, λ2 ; ζ = λ1 ζ2 + ζ1 + 1 and λ = λ1 λ2 + 2 . Accordingly, the cutting planes Γ♭ are described in terms of the foregoing points as follows:

3.3. Computation of the effective elastic isotropic moduli of HFPCs To calculate the two effective elastic isotropic moduli of HFPCs, the remote uniform strain boundary conditions (14) with ε0 specified by (19) and (20) are considered and prescribed on the whole external surface ∂Ω of the RVE Ω. The displacement field is then determined by solving the discrete governing Eq. (33). Hereafter, we can calculate the macroscopic stress and strain tensors through (15) and (16). Finally, the effective bulk and shear moduli κ* and μ* are obtained by using (21).

Γ1 ∈ £1{x10 , x 20 , x30}, Γ2 ∈ £2 {x10 , x30 , x40}, …, Γλ1 ∈ £ λ1 {x10 , x λ01+ 1, x 20} , Γλ1+ 1 ∈ £ λ1+ 1 {x 20 , x λ01+ 2 , x λ01+ 3, x30} , Γλ1+ 2 ∈ £ λ1+ 2 {x30 , x λ01+ 3, x λ01+ 4 , x40}, …, Γ2λ1 ∈ £2λ1 {x λ01+ 1, x 20λ1+ 1, x λ01+ 2 , x 20}, …, Γλ1 λ2 + 1 ∈ £ λ1 λ2 + 1 {x λ01 (λ2 − 1) + 2 , xλ0 , x λ01 (λ2 − 1) + 3} ,

4. Numerical examples and discussions

Γλ1 λ2 + 2 ∈ £ λ1 λ2 + 2 {x λ01 (λ2 − 1) + 3, xλ0 , x λ01 (λ2 − 1) + 4}, …,

In this section, we first elaborate an efficient and robust numerical procedure to generate RVEs for HFPCs such that the linearly elastic isotropy is ensured at the macroscopic scale. Then, the size- and shapedependent effects of HFPCs are examined.

Γm ∈ £m {x λ01 λ2 + 1, xλ0 , x λ01 (λ2 − 1) + 2} , where, for example, £ λ1+ 1{x 20 , sisting of the points x 20 , x λ01+ 2 ,

x λ01+ 2 , x λ01+ 3, x30} x λ01+ 3 and x30 .

(37) represents a plane con-

Further, each Γ♭ mathematically verifies the relation

n♭·(x − x♭0) = 0 for any x ∈ £♭.

4.1. A numerical procedure of generating RVEs for HFPCs

(38)

Above, x♭0 refers to the center point of Γ♭, and n♭ denotes its unit normal vector satisfying the following equation:

Let us consider a cubic domain Ω of side length w bounded by its outer surface ∂Ω (Fig. 4). The outward unit vector normal to ∂Ω is symbolized by nR. The origin O of the global Cartesian coordinate system is chosen to coincide with the center of Ω. The whole domain Ω is constituted by crystal materials of a certain kind.

n♭ =

(£♭ {2} − £♭{1}) × (£♭ {3} − £♭{1}) [(£♭{2} − £♭{1}) × (£♭{3} − £♭{1})]·[(£♭{2} − £♭{1}) × (£♭{3} − £♭{1})] (39)

with × and · separately standing for the cross and dot products. All interfaces Γ♭ are assumed to be characterized by the interfacial relations (6) and (7).

4.1.1. Selection of the cutting planes We first use m finite-sized planes Γ♭ (♭ = 1, 2, …, m ) to divide the crystal domain Ω shown in Fig. 4 so as to produce a simplified but essentially similar substitution for a real RVE of a HFPC. It is known that a cubic model consisting of a spherical isotropic inhomogeneity embedded centrally inside a homogeneous isotropic medium exhibits linear isotropic elasticity. Therefore, we define a fictive spherical surface Γ f of radius R0 smaller than 0.5w, whose center is located at the one of Ω. Then, λ points are extracted from Γ f in the following manner such that the overall isotropy is approximately guaranteed:

4.1.2. Evaluation of the thickness of each binding interface Since the exact thickness of each binding interphase is hardly to be determined analytically or experimentally, we now establish a simplified way to identify its value based on the aforementioned uniformthickness hypothesis and the inclusion volume fraction fp. Let us revisit the crystal domain Ω mentioned above, and define a reference thickness href for all embedded interfaces Γ♭. The thickness of interface elements in the rth particle Ωr is expressed as

hr _ i = cr _ i href

(40)

with cr _ i denoting the relative ratio of the thickness hr _ i of Γr _ i to href. The reference thickness href is assessed by invoking the definition of fp and performing some simple manipulations:

href =

(1 − fp ) V m

θ

∑r = 1 ∑i =r 1 cr _ i Ar _ i

. (41)

where Ar _ i represents the area of Γr _ i and V is the volume of Ω. Finally, the thickness of each interface element is retrieved by using (40). In what follows, we apply the just presented strategy of constructing an RVE with cutting planes for a certain HFPC and prescribe the boundary conditions on ∂Ω for estimating the overall elastic moduli of composites (Hornby et al., 1994; Nemat-Nasser and Hori, 1999).

Fig. 4. A cubic model composed of crystal materials of a certain kind. 39

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

4.1.3. Isotropy examination with the quantitative error indicators Now, we propose a general method to check the isotropy of the effective elastic stiffness tensor * of Ω on the basis of the necessary and sufficient conditions in He (2004). First, we impose the uniform strain boundary condition (14) on ∂Ω with

ɛ0

Table 1 Mechanical moduli of the HMX crystals (Banerjee and Adams, 2004).

3κ * (n ) =

) sin ( )⎤⎥ ⎥ ) sin ( )⎥⎥, (ι = (η ⎥ cos ( ) ⎥ ⎦

2(η1 − 1) π ϱ1

η2 π ϱ2

17.7

0.21

2(η1 − 1) π

η2 π

ϱ1

ϱ2

2

− 1)ϱ1 + η1 + 1).

η2 π ϱ2

(43)

1 1 = , * ni nj nk nl N : *: N Mijkl

(44)

1 1 = , (i, j, k, l = 1, 2, 3), * nk nl I : *: N Miikl

(45)

max ι (E * (n(ι) )) − minι (E * (n(ι) )) , ∼ E*

(48)

The RVE Ω is meshed with the tetrahedrons of size 0.07w. The remote boundary conditions (14) with different ε0 (Hornby et al., 1994; NematNasser and Hori, 1999) are prescribed on the external surfaces ∂Ω of Ω to evaluate the effective elastic tensor *. Afterwards, the isotropy examination of Ω is carried out by the procedure given in the last subsection. In particular, use is made of ϱ1 = ϱ2 = 40 . Fig. 5 plots the resulting relative error indicators defined by (46) and (47) versus the different numbers of cutting planes. It is observed that, as the number of cutting planes increases, both of the relative error indicators eE and eκ diminish dramatically at first, and then show much slower reduction after reaching a critical point, say the point {λ1 = 9, λ2 = 8} in Fig. 5 (d). Further, eE and eκ for all instances under consideration are smaller than 1% when λ1 is equal to or larger than 12. Thus, we can consider the models in question as being suitable for the HFPCs exhibiting the isotropic elasticity at the macroscopic scale. To qualitatively examine if the predicted effective moduli of the foregoing RVEs are correct, we have used the analytical estimates for the effective moduli of a cubic HMX model of side length w containing a spherical particle of radius R0 by the generalized self-consistent method (Christensen and Lo, 1979; Duan et al., 2007b). Precisely, the spherical particle is bonded to the outer domain through the spring-type interface characterized by (6) and (7). Both the numerical and analytical estimates for the effective moduli of all previous four instances are tabu∼ lated in Table 2. It is observed that E * and κ͠ * of each instance are pretty close to the relevant analytical ones where the spherical particle intervenes. Further, they vary in a similar way as Ea* and κa* when different binder materials are considered. The correctness of the numerical ∼ predictions is thus verified. Moreover, little difference between E * and Ea* (or κ͠ * and κa* ) in the different cases examined is observed. Fig. 6 shows the variations of the Young’s and bulk moduli, E* and κ*, with the orientation defined by the unit vector n(ι), when the binder material of Case 1 is used. From Fig. 5(a), it is known that the maximum relative error of E* (or κ*) among all considered directions takes place when the parameters {λ1 = 12, λ2 = 11} are selected and its value is less than 7% (or 4%). Hence, the orientation distribution functions of E* and κ* in Fig. 6 are graphically close to being spherical. Note however that the orientation distributions of E* and κ* with {λ1 = 15, λ2 = 14 } are closer to a sphere than the ones with {λ1 = 3, λ2 = 2 } since its maximum relative error is smaller (c.f. Fig. 5(a)), even though this is hard to be seen graphically. In regard to the macroscopic elastic isotropy, the highly filled RVE model generated with {λ1 = 15, λ2 = 14 } is a better one for HFPCs in question.

with the fourth-order effective elastic compliance tensor * = (*)−1 and N = n ⊗ n . As proved in He (2004), a necessary and sufficient condition for * or * to be isotropic is that both E*(n) and κ*(n) are independent of n. However, it is quite difficult to generate a perfectly isotropic model through introducing cutting planes as presented above. Therefore, the following two error indicators, namely eE and eκ, are employed to quantitatively characterize the relative differences of the effective Young’s and bulk moduli when varying the direction n(ι):

eE =

HMX

Case 2: Eb = 5 × 10−4Ep, vb = 0.45; Eb = 0.5Ep, vb = 0.1; Case 3: Eb = 0.5Ep, vb = 0.45. Case 4:

Above, ϱ1 and ϱ2, respectively, specify the numbers of points along the longitudinal and latitudinal directions; symbols η1 and η2 are non-negative integers varying separately in the ranges of [1, ϱ1] and [0, ϱ2]. As shown in He (2004), in the general case where * is anisotropic, the orientation distribution functions of the effective Young’s modulus E* and bulk modulus κ* can be defined as follows:

E * (n ) =

Poisson’s ratio vp

Case 1: Eb = 5 × 10−4Ep, vb = 0.1;

where ε0 is a non-zero scalar and n is a unit vector defining the direction along which the uniaxial strain is applied. In the present paper, the subsequent set of directions is chosen to test the isotropy of *:

( (

Young’s modulus Ep (GPa)

(42)

= ε0 (n ⊗ n )

⎡ cos ⎢ ⎢ n(ι) = ⎢ sin ⎢ ⎢ ⎢ ⎣

Crystal

(46)

max ι (κ * (n(ι) )) − minι (κ * (n(ι) )) , (47) κ͠ * ∼* η η η= E = ∑υ = 1 E * (n(υ) )/ η , κ͠ * = ∑υ = 1 κ * (n(υ) )/ η and with ϱ1ϱ2 − 2(ϱ1 − 1). By now, all preliminaries to generating and to justifying an RVE model showing approximately isotropy at the macro-scale have been presented. eκ =

4.2. Discussions on the isotropy of * The RVE Ω in the last subsection is now used for a certain HFPCto discuss the isotropy of *. The model parameters are selected as w = 10−5 m, R 0 = 0.25w and fp = 95%. The number of cutting planes in Eq. (37) is gradually increased so that the isotropy of * is verified extensively. Here, λ1 varies from 3 to 15, and λ2 = λ1 − 1. All embedded interfaces are characterized by the interfacial relations (6) and (7) where the thickness parameter h is taken to be the same and determined via the method presented in the last subsection. The RVE Ω is supposed to be made of the HMX material whose elastic moduli are listed in Table 1. To cover a large portion of the binder materials of practical interest (Banerjee and Adams, 2004; Ramakrishna et al., 2001; Wolff et al., 2014; Stober et al., 2001; Kaully et al., 2008), the elastic moduli of Γ♭s are selected as

4.3. Particle size effects on the effective isotropic moduli of HFPCs As concluded in the last subsection, a cubic domain intercepted with λ1 ≥ 12 can be considered approximately as an isotropic RVE of the HFPC. The HMX domain Ω of side length w with the cutting parameters {λ1 = 12, λ2 = 11} (Fig. 7) is employed to study the size effect of HFPCs in what follows. The side length w this time increases from 10−8 m to 10−3 m, and the radius of the fictive sphere is specified as R 0 = 0.25w . All cutting planes Γ♭ are characterized by the interfacial relations (7) and (8) with the same thickness parameter h and are determined on the basis of (36)–(39). To test the size-dependent behavior of HFPCs, two 40

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

Fig. 5. Relative errors of the indicators defined in (46) and (47) with respect to the different numbers of cutting planes and various binder materials.

negligible. As fp enlarges, say fp = 95%, κ* and μ* become accentuated simultaneously in view of the stiffness contributions from the reinforcement crystals. In summary, the crystal particle size has negligible effect on the overall isotropic moduli of HFPCs when the inclusion volume fraction fp is unchanged, and this particular property can be adopted to improve the intrinsic behaviors of practical HFPCs.

Table 2 Numerically predicted effective moduli of the designated RVE together with the analytical ones of an identical RVE containing a spherical particle.

Material

Interception

Case Case Case Case

λ1 = λ1 = λ1 = λ1 =

1 2 3 4

12, 12, 12, 12,

λ2 λ2 λ2 λ2

= = = =

11 11 11 11

Numerical moduli

Analytical moduli

∼ E * (GPa)

κ͠ * (GPa)

Ea* (GPa)

κa* (GPa)

16.3786 16.3790 16.5644 16.7473

9.1517 9.1522 9.2953 9.5109

15.5224 15.5331 17.1024 17.2461

8.9019 8.9095 9.7779 10.0349

4.4. Particle shape effects on the effective isotropic moduli of HFPCs The crystal particles of practical importance are arbitrarily shaped and numerous particles exist simultaneously in HFPCs. We now evaluate the effective isotropic moduli of these heterogeneous materials by the elaborated computational approach. The cubic domain Ω in Fig. 4 is used, and m planes Γ♭ (♭ = 1, 2, …, m ) defined by Eq. (38) are applied to divide Ω. More precisely, all cutting planes are considered to be infinite this time (Fig. 9), and characterized by Eqs. (6) and (7). The relevant geometrical parameters of each Γ♭ satisfy (36) and (39). The thickness h involved in characterizing Γ♭ is given by h♭ = ♭href / m , and determined by the previous numerical procedure. The model parameters are specified as w = 10−5 m, R 0 = 0.25w, fp = 90% and 95%, {λ1 = 12, λ2=11} and {λ1 = 15, λ2=14} . All bulk subdomains are assumed to be constituted by the HMX material specified previously, and all Γ♭s are supposed as being made of the binder material of the same kind. Similar to the former subsection, the elastic moduli of the binder material are also chosen as Eb = (5 × 10−4 ∼ 0.5) Ep and vb = 0.1. The entire domain of Ω is meshed with the tetrahedral elements of size 0.07w. The remote boundary conditions (14) with (19) and (20) are enforced on the external surfaces ∂Ω of Ω, and the effective isotropic

different inclusion volume fractions, i.e. fp = 0.9 and 0.95, are considered for each configuration of the RVE in question. The whole bulk domain Ω is considered as being made of the HMX crystal whose elastic coefficients are specified in Table 1. The elastic moduli of the binder material forming all embedded interfaces Γ♭ are set as Eb = (5 × 10−4 ∼ 0.5) Ep and vb = 0.1. The entire RVE Ω is discretized by the tetrahedral elements of size 0.07w, and the displacement-controlled boundary conditions defined by (14) with (19) and (20) are successively enforced on ∂Ω of Ω. Then, the boundary value problem under consideration is solved by the computational approach elaborated, and the overall isotropic moduli are calculated based on the numerical procedure developed in the last section. The resulting predictions of the overall isotropic moduli κ* and μ* against the RVE side length w are presented in Fig. 8. It is seen that the particle size shows so limited influence on both κ* and μ* when the inclusion volume fraction fp remains constant that its effect is 41

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

Fig. 6. The Young’s and bulk moduli orientation distribution functions E* (n(ι) ) and κ* (n(ι) ) of the designed RVE with the binder material of Case 1 and two distinct cutting strategies: these functions are almost independent of the orientation and, thus, their geometric representations are close to being spherical.

e κ* =

|κ * {λ1 = 12, λ1 = 11} − κ * {λ1 = 15, λ1 = 14}| , κ * {λ1 = 12, λ1 = 11}

(49)

e μ* =

|μ* {λ1 = 12, λ1 = 11} − μ* {λ1 = 15, λ1 = 14}| . μ* {λ1 = 12, λ1 = 11}

(50)

The relevant relative error indicator values are tabulated in Table 4. It is seen that both κ* and μ* are strongly dependent on the crystal shape when Eb is much smaller than Ep (e.g. Eb = 5 × 10−4Ep ), and such effect gradually diminishes with Eb enlarging and eventually vanishes when Eb ≈ Ep due to the spring-like interfacial behavior governed by (6) and (7). Hence, the overall elasticity of HFPCs is strongly affected not only by the crystals but also by the binder materials due to Eb/Ep ≪ 1. Further, the shape-dependent effect of HFPCs with fp = 90%, is more obvious than those with fp = 95%. This is ascribed to the fact the shape-dependent effect of HFPCs is induced by the binder material. The lower the inclusion volume fraction fp is, the more the binder material in HFPCs exists and thus, the stronger the shape-dependent effect of HFPCs will show. In brief, attention should be paid to the influence of the particle shapes on the overall moduli of HFPCs when synthesizing new HFPCs or modify their inherent properties.

Fig. 7. Schematic diagram of the highly filled RVE to study the size effect.

moduli of Ω are calculated by the computational approach elaborated. The resulting predictions of κ* and μ* with different interfacial parameters are listed in Table 3. It is observed that both κ* and μ* are closely related to the HMX crystal shape as fp is kept constant. To describe this shape-dependent effect explicitly, we introduce the relative error indicators below to characterize the difference of the effective moduli obtained by the foregoing two cutting strategies. 42

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

Fig. 8. Particle size effects on the overall isotropic moduli of HFPCs with different fbs and various Ebs. Table 3 Effective elastic isotropic moduli of the RVEs containing differently-shaped particles. Effective moduli fp

Interception

90%

λ1 = 12, λ2 = 11

Eb/Ep

κ*/κp

μ*/μp

5 × 10−4

0.0128

0.0197

5 × 10−2 0.5

0.4390

0.5620

0.8374 0.0109

0.9026 0.0166

5 × 10−4

λ1 = 15, λ2 = 14

5 × 10−2 0.5

5 × 10−4 95%

λ1 = 12, λ2 = 11

5 × 10−2 0.5

5 × 10−4

λ1 = 15, λ2 = 14

Fig. 9. Schematic diagram of a highly filled RVE containing multiple particles of different shapes.

5 × 10−2 0.5

0.4274

0.5498

0.8287 0.0229

0.8973 0.0349

0.5848

0.7038

0.8965 0.0208

0.9410 0.0311

0.5729

0.6932

0.8894

0.9373

Table 4 Values of the relative difference indicators by (49) and (50). Effective moduli

5. Concluding remarks fp

In this study, the spring-type imperfect interface has been recalled with its physical background, and used to characterize the thin binding interphases in HFPCs with the particle volume fraction higher than 90%. The macroscopic quantities and effective moduli over an RVE accounting for the interfacial displacement discontinuity are reformulated. Hereafter, a computational approach has been developed with the aid of the extended finite element and level set methods to

90%

Eb/Ep

e κ*

e μ*

5 × 10−4

0.1484

0.1574

5 × 10−2 0.5

0.0264

0.0217

0.0104 0.0917

0.0059 0.1089

0.0203

0.0151

0.0079

0.0039

5 × 10−4 95%

43

5 × 10−2 0.5

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

evaluate the effective elastic moduli of HFPCs containing crystal particles of arbitrary shapes. A numerical scheme has in particular been proposed to generate RVEs so as to ensure their isotropy of macroscopic elasticity. With the help of the computational approach elaborated, the size and shape effects of crystal particles on the effective elastic isotropic moduli of HFPCs have been examined in detail. It is found that the effective isotropic moduli of HFPCs show negligible dependence on the particle size as the inclusion volume fraction fp is kept constant, and both of the effective elastic isotropic moduli gradually increase with the increment of fp. Concerning the HFPCs with the identical fp but composed of differently shaped particles, the effective bulk and shear moduli of HFPCs vary dramatically, particularly when the Young’s modulus Eb of the binder phase is much smaller than the one Ep of the particle phase. Hence, the overall moduli of HFPCs are strongly shapedependent. It is of particular importance to take these features into account when synthesizing new HFPCs and improving the properties of existing ones.

Duan, H.L., Yi, X., Huang, Z.P., Wang, J., 2007. A unified scheme for prediction of effective moduli of multiphase composites with interface effects. Part i: theoretical framework. Mech. Mater. 39 (1), 81–93. https://doi.org/10.1016/j.mechmat.2006. 02.009. Duan, H.L., Yi, X., Huang, Z.P., Wang, J., 2007. A unified scheme for prediction of effective moduli of multiphase composites with interface effects: part II - application and scaling laws. Mech. Mater. 39 (1), 94–103. https://doi.org/10.1016/j.mechmat. 2006.02.010. Fang, X.Q., Tian, J.Y., 2018. Elastic-adhesive interface effect on effective elastic moduli of particulate-reinforced asphalt concrete with large deformation. Int. J. Eng. Sci. 130, 1–11. https://doi.org/10.1016/j.ijengsci.2018.05.007. Farsad, M., Vernerey, F.J., Park, H.S., 2010. An extended finite element/level set method to study surface effects on the mechanical behavior and properties of nanomaterials. Int. J. Numer. MethodsEng. 84 (12), 1466–1489. https://doi.org/10.1002/nme. 2946. Fries, T.P., 2008. A corrected XFEM approximation without problems in blending elements. Int. J. Numer. MethodsEng. 75 (5), 503–532. https://doi.org/10.1002/nme. 2259. Gao, X., Fan, Z., Zhang, J., Liu, S., 2017. Micromechanical model for asphalt mixture coupling inter-particle effect and imperfect interface. Construct. Build. Mater. 148, 696–703. https://doi.org/10.1016/j.conbuildmat.2017.05.015. Gonzalez-Gutierrez, J., Cano, S., Schuschnigg, S., Kukla, C., Sapkota, J., Holzer, C., 2018. Additive manufacturing of metallic and ceramic components by the material extrusion of highly-filled polymers: a review and future perspectives. Materials 11 (5), 840. https://doi.org/10.3390/ma11050840. Gu, S.T., Liu, J.T., He, Q.C., 2014. Size-dependent effective elastic moduli of particulate composites with interfacial displacement and traction discontinuities. Int. J. Solids Struct. 51 (13), 2283–2296. https://doi.org/10.1016/j.ijsolstr.2014.02.033. Harari, I., Dolbow, J., 2010. Analysis of an efficient finite element method for embedded interface problems. Comput. Mech. 46 (1), 205–211. https://doi.org/10.1002/nme. 3345. Hashin, Z., 1992. Extremum principles for elastic heterogenous media with imperfect interfaces and their application to bounding of effective moduli. J. Mech. Phys. Solids 40 (4), 767–781. https://doi.org/10.1016/0022-5096(92)90003-K. He, G., Liu, J., Gong, F., Lin, C., Yang, Z., 2018. Bioinspired mechanical and thermal conductivity reinforcement of highly explosive-filled polymer composites. Compos. Part A 107, 1–9. https://doi.org/10.1016/j.compositesa.2017.12.012. He, Q.C., 2004. Characterization of the anisotropic materials capable of exhibiting an isotropic young or shear or area modulus. Int. J. Eng. Sci. 42 (19-20), 2107–2118. https://doi.org/10.1016/j.ijengsci.2004.04.009. Heukamp, F.H., Lemarchand, E., Ulm, F.J., 2005. The effect of interfacial properties on the cohesion of highly filled composite materials. Int. J. Solids Struct. 42 (1), 287–305. https://doi.org/10.1016/j.ijsolstr.2004.07.007. Hornby, B.E., Schwartz, L.M., Hudson, J.A., 1994. Anisotropic effective-medium modelling of the elastic properties of shales. Geophysics 59 (10), 1570–1583. https://doi. org/10.1190/1.1443546. Javili, A., Steinmann, P., Mosler, J., 2017. Micro-to-macro transition accounting for general imperfect interfaces. Comput. Methods Appl. Mech.Eng. 317, 274–317. https://doi.org/10.1016/j.cma.2016.12.025. Ju, J.W., Chen, T.M., 1994. Effective elastic moduli of two-phase composites containing randomly dispersed spherical inhomogeneities. Acta Mechanica 103 (1–4), 123–144. https://doi.org/10.1007/BF01180222. Kaully, T., Siegmann, A., Shacham, D., 2008. Mechanical behavior of highly filled natural CaCO3 composites: effect of particle size distribution and interface interactions. Polymer Compos. 29 (4), 396–408. https://doi.org/10.1002/pc.20435. Kauly, T., Keren, B., Siegmann, A., Narkis, M., 1996. Highly filled thermoplastic composites. II: effects of particle size distribution on some properties. Polymer Compos. 17 (6), 806–815. https://doi.org/10.1002/pc.10673. Kauly, T., Keren, B., Siegmann, A., Narkis, M., 1997. Highly filled particulate thermoplastic composites: part I packing density of irregularly shaped particles. J. Mater. Sci. 32 (3), 693–699. https://doi.org/10.1023/A:1018543920255. Larson, N.M., Zok, F.W., 2018. In-situ 3D visualization of composite microstructure during polymer-to-ceramic conversion. Acta Materialia 144 (1), 579–589. https:// doi.org/10.1016/j.actamat.2017.10.054. Lemarchand, E., Ulm, F.J., Dormieux, L., 2002. Effect of inclusions on friction coefficient of highly filled composite materials. J. Eng. Mech. 128 (8), 876–884. https://doi.org/ 10.1061/(ASCE)0733-9399(2002)128:8(876). Liang, S., Wei, Y., Wu, Z., 2017. Multiscale modeling elastic properties of cement-based materials considering imperfect interface effect. Construct. Build. Mater. 154, 567–579. https://doi.org/10.1016/j.conbuildmat.2017.07.196. Liu, J.T., Gu, S.T., He, Q.C., 2015. A computational approach for evaluating the effective elastic moduli of non-spherical particle reinforced composites with interfacial displacement and traction jumps. Int. J. Multiscale Comput.Eng. 13 (2), 123–143. https://doi.org/10.1615/IntJMultCompEng.2014011640. Moës, N., Dolbow, J., Belytschko, T., 1999. A finite element method for crack growth without remeshing. Int. J. Numer. MethodsEng. 46 (1), 131–150. https://doi.org/10. 1002/(SICI)1097-0207(19990910)46:1<131::AID-NME726>3.0.CO;2-J. Moshev, V.V., Garishin, O.C., 1993. Physical discretization approach to evaluation of elastic moduli of highly filled granular composites. Int. J. Solids Struct. 30 (17), 2347–2355. https://doi.org/10.1016/0020-7683(93)90122-n. Nemat-Nasser, S., Hori, M., 1999. Micromechanics: Overall Properties of Heterogeneous materials (Second Revised edition). Amsterdam. : Elsevier press Osher, S., 1993. A level set formulation for the solution of the Dirichlet problem for Hamilton-Jacobi equations. SIAM J. Math. Anal. 24 (5), 1145–1152. https://doi.org/ 10.1137/0524066. Pichler, C., Lackner, R., 2009. Upscaling of viscoelastic properties of highly-filled

Acknowledgements This work was funded by the National Natural Science Foundation of China (grant number 51505441, 51775459, 11572266 and 51305362) and the Young Scientific Innovation Team of Science and Technology of Sichuan province (Grand No. 2017TD0017). Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.mechmat.2019.04.022 References Azadi, P., Nunnari, S., Farnood, R., Kortschot, M., Yan, N., 2009. High speed compression of highly filled thin composites: effect of binder content and stiffness. Progr. Organic Coatings 64 (4), 356–360. https://doi.org/10.1016/j.porgcoat.2008.07.021. Babu, I., de With, G., 2014. Highly flexible piezoelectric 0–3 PZT-PDMS composites with high filler content. Compos. Sci. Technol. 91 (31), 91–97. https://doi.org/10.1016/j. compscitech.2013.11.027. Banerjee, B., Adams, D.O., 2004. On predicting the effective elastic properties of polymer bonded explosives using the recursive cell method. Int. J. Solids Struct. 41 (2), 481–509. https://doi.org/10.1016/j.ijsolstr.2003.09.016. Banerjee, B., Cady, C.M., Adams, D.O., 2003. Micromechanics simulations of glass-estane mock polymer bonded explosives. Modell. Simul. Mater. Sci.Eng. 11 (4), 457–475. https://doi.org/10.1088/0965-0393/11/4/304. Belytschko, T., Black, T., 1999. Elastic crack growth in finite elements with minimal remeshing. Int. J. Numer. MethodsEng. 45 (5), 601–620. https://doi.org/10.1002/ (sici)1097-0207(19990620)45:5<601::aid-nme598>3.0.co;2-s. Belytschko, T., Gracie, R., Ventura, G., 2009. A review of extended/generalized finite element methods for material modeling. Modell. Simul. Mater. Sci.Eng. 17, 043001. https://doi.org/10.1088/0965-0393/17/4/043001. Bengtsson, M., Baillif, M.L., Oksman, K., 2007. Extrusion and mechanical properties of higly filled cellulose fibre-polypropylene composites. Compos. 38 (8), 1922–1931. https://doi.org/10.1016/j.compositesa.2007.03.004. Benveniste, Y., 1985. The effective mechanical behaviour of composite materials with imperfect contact between the constituents. Mech. Mater. 4 (2), 197–208. https:// doi.org/10.1016/0167-6636(85)90016-x. Benveniste, Y., 2006. A general interface model for a three-dimensional curved thin anisotropic interphase between two anisotropic media. J. Mech. Phys. Solids 54 (4), 708–734. https://doi.org/10.1016/j.jmps.2005.10.009. Benvenuti, E., 2008. A regularized XFEM framework for embedded cohesive interfaces. Comput. Methods Appl. Mech.Eng. 197 (49–50), 4367–4378. https://doi.org/10. 1016/j.cma.2008.05.012. Bövik, P., 1994. On the modelling of thin interface layers in elastic and acoustic scattering problems. Q. J. Mech. Appl.Math. 47 (1), 17–42. https://doi.org/10.1093/qjmam/ 47.1.17. Christensen, R.M., Lo, K.H., 1979. Solutions for effective shear properties in three phase sphere and cylinder models. J. Mech. Phys. Solids 27 (4), 315–330. https://doi.org/ 10.1016/0022-5096(79)90032-2. Dolbow, J., Harari, I., 2009. An efficient finite element method for embedded interface problems. Int. J. Numer. MethodsEng. 78, 229–252. https://doi.org/10.1002/nme. 2486. Duan, H.L., Wang, J., Huang, Z.P., Karihaloo, B.L., 2005. Size-dependent effective elastic constants of solids containing nano-inhomogeneities with interface stress. J. Mech. Phys. Solids 53 (7), 1574–1596. https://doi.org/10.1016/j.jmps.2005.02.009.

44

Mechanics of Materials 135 (2019) 35–45

J.-T. Liu, et al.

Wang, J., Luo, J., 2018. Micromechanical study on the effective elastic moduli of polymer-bonded explosives with imperfect interfaces. J. Energetic Mater. 36 (3), 278–290. https://doi.org/10.1080/07370652.2017.1392648. Wolff, M.F.H., Salikov, V., Antonyuk, S., Heinrich, S., Schneider, G.A., 2014. Novel, highly-filled ceramic-polymer composites synthesized by a spouted bed spray granulation process. Compos. Sci. Technol. 90 (10), 154–159. https://doi.org/10.1016/j. compscitech.2013.11.006. Xiao, J.J., Wang, W.R., Chen, J., Ji, G.F., Zhu, W., Xiao, H.M., 2012. Study on structure, sensitivity and mechanical properties of HMX and HMX-based PBXs with molecular dynamics simulation. Comput. Theor. Chem. 999 (1), 21–27. https://doi.org/10. 1016/j.comptc.2012.08.006. Yvonnet, J., Le-Quang, H., He, Q.C., 2008. An XFEM/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites. Comput. Mech. 42 (1), 119–131. https://doi.org/10.1007/s00466008-0241-y. Zhu, Q.Z., Gu, S.T., Yvonnet, J., Shao, J.F., He, Q.C., 2011. Three-dimensional numerical modelling by XFEM of spring-layer imperfect curved interfaces with applications to linearly elastic composite materials. Int. J. Numer. MethodsEng. 88 (4), 307–328. https://doi.org/10.1002/nme.3175.

composites: investigation of matrix-inclusion-type morphologies with power-law viscoelastic material response. Compos. Sci. Technol. 69 (14), 2410–2420. https:// doi.org/10.1016/j.compscitech.2009.06.008. Pichler, C., Lackner, R., Aigner, E., 2012. Generalized self-consistent scheme for upscaling of viscoelastic properties of highly-filled matrix-inclusion composites - application in the context of multiscale modeling of bituminous mixtures. Compos. 43 (2), 457–464. https://doi.org/10.1016/j.compositesb.2011.05.034. Politi, I., McHugh, L.E.J., Al-Fodeh, R.S., Fleming, G.J.P., 2018. Modification of the restoration protocol for resin-based composite (RBC) restoratives (conventional and bulk fill) on cuspal movement and microleakage score in molar teeth. Dental Mater. 34 (9), 1271–1277. https://doi.org/10.1016/j.dental.2018.05.010. Qu, J., 1993. The effect of slightly weakened interfaces on the overall elastic properties of composite materials. Mech. Mater. 14 (4), 269–281. https://doi.org/10.1016/01676636(93)90082-3. Ramakrishna, S., Mayer, J., Wintermantel, E., Leong, K.W., 2001. Biomedical applications of polymer-composite materials: a review. Compos. Sci. Technol. 61 (9), 1189–1224. https://doi.org/10.1016/S0266-3538(00)00241-4. Stober, T., Gilde, H., Lenz, P., 2001. Color stability of highly filled composite resin materials for facings. Dental Mater. 17 (1), 87–94. https://doi.org/10.1016/S01095641(00)00065-8.

45