Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2021.017589

ARTICLE

Effective Elastic Properties of 3-Phase Particle Reinforced Composites with Randomly Dispersed Elastic Spherical Particles of Different Sizes

Dedicated to Professor Karl Stark Pister for his 95th birthday

Yu-Fu Ko1,* and Jiann-Wen Woody Ju2

1Civil Engineering and Construction Engineering Management, California State University, Long Beach, CA 90840-5101, USA
2Civil and Environmental Engineering, University of California, Los Angeles, CA 90095-1593, USA
*Corresponding Author: Yu-Fu Ko. Email: yu-fu.ko@csulb.edu
Received: 22 May 2021; Accepted: 05 August 2021

1  Introduction

Composite materials can considerably enhance the following material properties including strength, stiffness, thermal insulation, thermal conductivity, fatigue life, acoustical insulation, wear resistance, corrosion resistance, etc. To achieve targeted engineering performance, composite materials consisting of two or more different materials to macroscopically form new materials are needed. The “inclusions” in composites can be in the forms of whiskers, fibers, and particulates. The “matrix” in composites is the binder material. The matrix material provides the support and protection to the inclusions. The matrix material also transfers stresses and strains through inclusion/matrix interfaces under complex 3-dimensional loading. Limited mechanical properties of conventional particle reinforced composites that contain a single material type of inclusions are observed such as tensile strength, compressive strength, impact resistance.

Alternatively, hybrid particle-reinforced composites, containing several different sizes and/or materials of particles into a matrix, demonstrate superior and excellent mechanical properties and have been largely employed in engineering applications, e.g., aerospace, civil engineering, automobile industries, and military equipment. The mechanical behaviors of hybrid composites, compared with conventional composites, are improved by the weighted summation of the individual particles with different material properties and distinct sizes. Therefore, enhanced engineering performance and cost could be obtained by appropriate engineering material design [1–5].

Numerous experimental research works have exhibited that mechanical behaviors of particle reinforced composites are controlled by particle sizes, micro-structural morphology, and interfacial properties between matrix material and particles [6–17]. Several theoretical methods have been developed in literatures to derive the effective elastic properties of multiphase composites such as variational methods, effective medium methods, direct micromechanical methods, and finite element methods.

Variational methods utilize linear comparison composites or variational principles to acquire mathematical upper and lower bounds for effective elastic properties of composites [18–27]. Hashin's bounds are referred to as the 2-point bounds. Furthermore, the “improved” higher-order mathematical bounds considering the statistical micro-structural information of composites were proposed [28–33]. For instance, Silnutzer [28] proposed improved bounds, referred to as 3-point bounds, on the effective in-plane shear modulus and bulk modulus. The 3-point bounds are narrower when compared with the 2-point bounds.

Effective medium methods employ effective medium to predict the effective elastic properties of composites [34–39]. Mori-Tanaka method, differential scheme, self-consistent method, and generalized self-consistent method are among the effective medium methods. These methods considered only the volume fractions and geometries of inclusions. Conversely, probabilistic distributions or spatial locations of inclusions are not considered. Therefore, the effective medium methods are best suited for low volume fractions of inclusions or some limited particular configurations.

Direct micromechanical methods consider specific geometric configurations of inclusions dispersed in the matrix and utilize approximations to determine effective elastic properties of composites with randomly located and interacting inclusions. Eshelby [40] proposed the renowned “Eshelby's equivalence principle” based on an ellipsoidal inclusion embedded in infinite matrix. Mura [41] primarily considered rigorous “local” micromechanics. Honein [42] proposed frameworks based on Kolosov-Muskhelishvili complex potentials to study circular inclusions in plane elastostatics. Moreover, the direct micromechanical methods are further utilized and studied [43]. Nonetheless, based on the above methods, “local” (not “overall”) field solutions could be achieved. To predict the effective properties of elastic multi-phase composites, Ju et al. [44,45], based on Eshelby [40], proposed a micromechanical ensemble volume average formulation, which is higher-order in volume fractions of inclusions ϕ. It is noted that randomly dispersed ellipsoidal, and spherical inclusions were considered. In addition, both “local” and “overall” ensemble volume averaged micromechanical field equations were derived by rigorously considering the inclusion interactions in formulations and utilizing homogenization technique. Based on the proposed micromechanics formulations, effective elastic properties of composites containing randomly dispersed spherical particles or randomly located circular fibers featuring same or distinct elastic properties of inclusions as well as same inclusion sizes were studied [46–49].

Finite element methods, i.e., numerical solutions, are obtained based on the “unit cell model” as well as assume particular periodic arrays of inclusions [50–53].

Finally, Ju and co-researchers, following the micromechanical frameworks [44,45], continue to study the effective elastoplastic with/without damage behaviors of advanced composite materials [54–81]. Additionally, Ko et al. [82,83] proposed new higher-order bounds on effective transverse elastic properties of 3-phase hybrid fiber-reinforced composites. The hybrid fiber-reinforced composites contain randomly located and unidirectionally aligned circular fibers featuring distinct fiber sizes and different elastic material properties.

Majority of research works on prediction of the effective elastic properties of particle reinforced composites focus on conventional composites containing one type of particle and did not consider the effects of different particle sizes. In addition, the “unit cell model,” commonly adopted within the 3-D RVE finite element models, assumes periodic arrays of particles, not randomly distributed.

The predictions of effective elastic properties of particle reinforced composites with different particle sizes are still limited in research works. The objective of this paper is to propose analytical frameworks based on the direct micromechanical methods to derive the effective elastic properties of 3-phase spherical particle reinforced composites featuring different particle sizes and same elastic material property of particle. All particles are considered randomly dispersed, non-intersection, and with perfect interfaces between matrix and particles. In addition, both “local” and “overall” ensemble volume averaged micromechanical field equations will be derived by rigorously considering the spherical particle interactions and the probabilistic spatial distribution of particles in formulations as well as utilizing homogenization technique with ensemble volume average approach.

Section 2 presents analytical local solutions of spherical particle interactions. The spherical particles are assumed to be elastic and randomly dispersed in the matrix material. 3-phase composites contain 2 inclusion phases featuring different particle sizes and same elastic properties. Consequently, Section 3 presents the derivations of the ensemble volume averaged eigenstrains based on the probabilistic particle interactions mechanism considering uniform radial distribution function (URDF). Additionally, 2 formulations with different orders, Formulation II and Formulation I, are proposed. In Section 4, the effective elastic properties of 3-phase composites consisting of spherical particles, which are randomly dispersed as well as featuring different particle sizes and same elastic material properties, are analytically obtained. In Section 5, our theoretical predictions are compared with available experimental data together with numerical examples to demonstrate the competence of the proposed frameworks. In Section 6, conclusions are presented.

2  Analytical Local Solutions of Spherical Particle Interactions

Consider a 3-phase particle reinforced composite consisting of randomly dispersed and spherical shaped elastic particles with different particle sizes embedded in an isotropic elastic matrix as shown in Fig. 1. The matrix material, spherical particles with radius a1, spherical particles with radius a2 are denoted as the 0th phase, the 1st phase, and the 2nd phase, respectively. μ0 and k0 are the shear modulus and bulk modulus of the matrix material. μ1 and k1 are the shear modulus and bulk modulus of the 1st phase particles. μ2 and k2 are the shear modulus and bulk modulus of the 2nd phase particles. Since elastic material properties of the 1st phase and the 2nd phase are the same, k=k1=k2 and μ=μ1=μ2 are assumed. Additionally, the linear elastic isotropic stiffness tensors for each phase can be expressed as

(Cζ)ijkl=λΩδijδkl+μζ(δikδjl+δilδjk),ζ=0,1,2(1)

where λζ and μζ denote the local Lamé coefficients for the ζth phase particles.

Figure 1: 2-particle interaction in a 3-phase composite

According to the eigenstrain concept of Eshelby [40], with the replacement of the spherical particles by the matrix material, the perturbed strain field ε′(x) generated by the spherical particles can be expressed in terms of the specified eigenstrain ε∗(x). Therefore, the linear elastic stiffness tensor Cζ for the ζth phase particles can be expressed as

Cζ:[ε0+ε′(x)]=C0:[ε0+ε′(x)−ε∗(x)],ζ=1,2(2)

where ε0 denotes the far-field loading induced uniform strain field in a homogeneous matrix material (without inhomogeneities). In this paper, the symbol “∙” is the dot product between two 4th-order tensors; the symbol “:” is the double-dot product between the 4th-order tensor and the 2nd-order tensor.

Due to the presence of the distributed eigenstrain ε∗(x), the perturbed strain generated in V per Eshelby [40] is

ε′(x)=∫VG(x−x′):ε∗(x′)dx′(3)

where V is volume of a representative volume element (RVE), x, x′∈V.

In addition, the 4th-order tensor G per Mura [41] is expressed as

Gijkl=18π(1−ν0)r′3Fijkl(−15,3ν0,3,3−6ν0,−1+2ν0,1−2ν0),i,j,k,l=1,2,3with r′≡x−x′;r′=‖r′‖(4)

where ν0 is the Poisson's ratio of the homogeneous matrix material.

We also define the components of the 4th-order tensor F in the Cartesian coordinates as (m = 1–6)

Fijkl(Bm)≡B1ni′nj′nk′nl′+B2(δiknj′nl′+δilnj′nk′+δjkni′nl′+δjlni′nk′)+B3δijnk′nl′     +B4δklni′nj′+B5δijδkl+B6(δikδjl+δilδjk)(5)

where n′≡r′/r′ and δij represents the Kronecker delta. The summation convention applies.

Making use of Eqs. (2) and (3),

−Ai:ε∗(x)=ε0+∫VG(x−x′):ε∗(x′)dx′,x∈V(6)

Ai=[Ci−C0]−1⋅C0(7)

By considering spherical particle interactions, Eq. (6) are rearranged as

−Ai:ε(i)∗(x)=ε0+∫ΩiG(x−x′):ε(i)∗(x′)dx′+∫ΩjG(x−x′):ε(j)∗(x′)dx′x∈Ωi,i≠j,i,j=1,2,(8)

where ε(i)∗(x′) represents the eigenstrain located at x′ in the ith phase spherical particles inside Ωi domain.

To obtain the 1st-order solution for the eigenstrain ε(i)∗0 for the ith phase spherical particles, spherical particle interactions are dropped in Eq. (8). The 1st-order formulation yields [44]

−Ai:ε(i)∗0=ε0+S:ε(i)∗0(9)

In addition, the 4th-order tensor S, known as interior-point Eshelby tensor of a spherical particle, is [41,58,59]

S=∫ΩiG(x−x′)dx′,x, x′∈Ωi(10)

The components of S for a spherical particle are [41]

Sijkl=115(1−v0)[(5v0−1)δijδkl+(4−5v0)(δikδjl+δilδjk)](11)

where v0 is the Poisson's ratio of the matrix material.

To obtain the effects of spherical particle interactions, Eq. (8) is subtracted by Eq. (9). Then, solving the integral equation

−Ai:d(i)∗(x)=∫ΩjG(x−x′)dx′:ε(j)∗0+∫ΩiG(x−x′):d(i)∗(x′)dx′       +∫ΩjG(x−x′):d(j)∗(x′)dx′, for x∈Ωi,  i≠j and i,j=1,2(12)

In addition, we define

d(i)∗(x)=ε(i)∗(x)−ε(i)∗0(13)

In order to find the correction of ε(i)∗(x) to account for higher-order spherical particle interactions, the 4th-order tensor G(x−x′) are expanded with respect to its center point xj inside the Ωj domain;

G(x−x′)=G(x−xj)−(x′−xj):[∇x⊗G(x−xj)]     +12[(x′−xj)⊗(x′−xj)]:[∇x⊗∇x⊗G(x−xj)]+⋯(14)

where the following equalities are employed.

∇x′⊗G(x−x′)=−∇x⊗G(x−x′)(15)

Substitution of Eq. (14) to Eq. (12) yields

−Ai:d(i)∗(x)=∫ΩjG(x−x′)dx′:ε(j)∗0+∫ΩiG(x−x′):d(i)∗(x′)dx′       +ΩjG(x−xj):d¯(j)∗(xj)−Ωjaj{∇x⊗G(x−xj)}:P¯(j)∗       +12Ωjaj2{∇x⊗∇x⊗G(x−xj)}:Q¯(j)∗+⋯with x∈Ωi, i,j=1,2(16)

where volume of a spherical particle in the ith phase is Ωi=4πai3/3 with the spherical particle radius ai; volume of a spherical particle in the jth phase is Ωj=4πaj3/3 with the spherical particle radius aj.

Moreover, the average fields shown in Eq. (16) are defined as

d¯(j)∗=1Ωj∫Ωjd(j)∗(x′)dx′(17)

P¯(j)∗=1Ωjaj∫Ωj(x′−xj)⊗d(j)∗(x′)dx′(18)

Q¯(j)∗=1Ωjaj2∫Ωj(x′−xj)⊗(x′−xj)⊗d(j)∗(x′)dx′(19)

The 3rd-order tensor P¯(j)∗ and the 4th-order tensor Q¯(j)∗ are the dipole and quadrapole of d¯(j)∗ in the Ωj domain, respectively.

It is noted that the leading order of P¯(j)∗ is the order of O(ρj4) instead of O(ρj3). It can be proved by substituting Eq. (18) into Eq. (16) and making use of the symmetry of spherical particles. In addition, ρj=aj/r with r defined as the distance between the centers of spherical particles in the ith phase and spherical particles in the jth phase.

By conducting volume averaging of Eq. (16) over the Ωi domain and dropping the terms with higher order moments, the  d¯(j)∗ accounting for local spherical particle interactions is approximately obtained.

By setting i equal to 1 and j equal to 2, we obtained Eqs. (20)–(23).

−A1:d¯(1)∗=G¯21:ε(2)∗0+S:d¯(1)∗+G¯22:d¯(2)∗+O(ρ8)(20)

in which

d¯(1)∗=1Ω1∫Ω1d(1)∗(x)dx;S=∫Ω1G(x−x′)dx′;for x∈Ω1, x′∈Ω1(21)

G¯21=1Ω1∫Ω1∫Ω2G(x−x′)dx′dx=130(1−ν0)[ρ23H1+(ρ25+ρ23ρ12)H2],for x∈Ω1,x′∈Ω2(22)

G¯22=∫Ω1G(x−x2)dx=130(1−ν0)(ρ23H1+ρ25H2),for x∈Ω1,x2∈Ω2(23)

Similarly, by setting i equal to 2 and j equal to 1, we obtained Eqs. (24)–(27).

−A2:d¯(2)∗=G¯12:ε(1)∗0+S:d¯(2)∗+G¯11:d¯(1)∗+O(ρ8)(24)

where

d¯(2)∗=1Ω2∫Ω2d(2)∗(x′)dx′;S=∫Ω2G(x−x′)dx′;for x∈Ω2,x′∈Ω2(25)

G¯12=1Ω2∫Ω2∫Ω1G(x−x′)dx′dx=130(1−ν0)[ρ13H1+(ρ15+ρ13ρ22)H2],for x∈Ω2,x′∈Ω1(26)

G¯11=∫Ω2G(x−x1)dx=130(1−ν0)(ρ13H1+ρ15H2),with x∈Ω2,x1∈Ω1(27)

and the components of H1 and H2 are rendered by

Hijkl1(x1−x2)≡5Fijkl(−15,3ν0,3,3−6ν0,−1+2ν0,1−2ν0)Hijkl2(x1−x2)≡3Fijkl(35,−5,−5,−5,1,1)(28)

Moreover, we define ρ1=a1/r, ρ2=a2/r, Ω1=4πa13/3 and Ω2=4πa23/3. It is noted that G¯11 in Eq. (27) and  G¯22 in Eq. (23) are the exterior-point Eshelby tensors of spherical particles. The leading-order is of the order O(ρ8) in Eqs. (20) and (24). It can be proved after truncating the terms with higher order moments and making use the fact that both Ωjaj{∇x⊗G(x−xj)} and P¯(j)∗ are of the order O(ρ4).

Furthermore, Eqs. (20) and (24) are further rearranged:

(A1+S1):d¯(1)∗+G¯22:d¯(2)∗=−G¯21:ε(2)∗0G¯11:d¯(1)∗+(A2+S2):d¯(2)∗=−G¯12:ε(1)∗0(29)

Eq. (29) can be solved to obtain

d¯(1)∗=[(G¯22)−1⋅(A1+S1)−G¯11⋅(A2+S2)−1]−1⋅[(A2+S2)−1⋅G¯12:ε(1)∗0−(G¯22)−1⋅G¯21:ε(2)∗0](30)

d¯(2)∗=[G¯22⋅(A1+S1)−1−(G¯11)−1⋅(A2+S2)]−1⋅[(G¯11)−1⋅G¯12:ε(1)∗0−G¯21⋅(A1+S1)−1:ε(2)∗0](31)

It can be proved that the leading order is of O(ρ3) for G¯11(A2+S)−1 and the leading order is of O(ρ−3) for (G¯22)−1(A1+S) in Eq. (30). Therefore, G¯11(A2+S)−1 is truncated.

In addition, ρ1≤1/2, ρ2≤1/2, and ρ1+ρ2≤1. Therefore, we obtained

d¯(1)∗=[(A1+S1)−1⋅G¯22⋅(A2+S2)−1⋅G¯12:ε(1)∗0−(A1+S1)−1⋅G¯21:ε(2)∗0](32)

Likewise, Eq. (31) can be expressed as

d¯(2)∗=[(A2+S2)−1⋅G¯11⋅(A1+S1)−1⋅G¯21:ε(2)∗0−(A2+S2)−1⋅G¯12:ε(1)∗0](33)

3  Averaged Eigenstrains of Spherical Particles by Ensemble Volume Average Approach

To find the solutions of averaged eigenstrains of spherical particles considering spherical particle interactions, ensemble volume average approach is adopted. The ensemble volume average process is stated as

⟨d¯(i)∗⟩(xi)=∫V−Ωid¯(i)∗(xi−xj)P(xj∣xi)dxj,i≠j(34)

i=1,j=2:⟨d¯(1)∗⟩(x1)=∫V−Ω1d¯(1)∗(x1−x2)P(x2∣x1)dx2(35)

i=2,j=1:⟨d¯(2)∗⟩(x2)=∫V−Ω2d¯(2)∗(x2−x1)P(x1∣x2)dx1(36)

where P(xj∣xi) denotes the conditional probability density function to discover the spherical particles centered at xj in the jth phase given the spherical particles centered at xi in the ith phase. In addition, P(xj∣xi), a 2-point probability density function, which is homogenous and 3-dimensional statistically isotropic, is considered in this paper. Thus, integration domain V in Eq. (34) can be determined as a spherical shape. Furthermore, Ωi is defined as the probabilistic “exclusion zone” for the spherical particles centered at xj in the jth phase. ⟨⋅⟩ represents the ensemble volume average.

In this paper, we adopt a statistically isotropic and uniform 2-point conditional probability density function known as uniform radial distribution function (URDF).

P(xj∣xi)={(Ni+Nj)V,if r^≥1,where r^≡r/(a1+a2),r>(a1+a2)0,otherwise.i,j=1,2,i≠j(37)

where r is the distance between the centers of the spherical particles in the ith phase and the spherical particles in the jth phase. Ni and Nj are the numbers of the ith phase spherical particles and the jth phase spherical particles, respectively. Ni/V and Nj/V are the number density of the ith phase spherical particles and the jth phase spherical particles in a composite, respectively.

After substituting Eq. (32) into Eq. (35),

〈d¯(1)∗〉(x1)=∫V−Ω1d¯(1)∗(x1−x2)P(x2∣x1)dx2    =∫V−Ω1[(A1+S1)−1⋅G¯22⋅(A2+S2)−1⋅G¯12:ε(1)∗0]P(x2∣x1)dx2      −∫V−Ω1[(A1+S1)−1⋅G¯22:ε(2)∗0]P(x2∣x1)dx2(38)

〈d¯(1)∗〉(x1)=∫(a1+a2)∞∫Ξ[P(x2∣x1)(A1+S1)−1⋅G¯22⋅(A2+S2)−1⋅G¯12dΞdr]:ε(1)∗0      −∫(a1+a2)∞∫Ξ[P(x2∣x1)(A1+S1)−1⋅G¯21dΞdr]:ε(2)∗0    =∫a1+a2∞∫Ξ[P(x2∣x1)(A1+S1)−1⋅G¯22⋅(A2+S2)−1⋅G¯12dΞdr]:ε(1)∗0      +∫2a1∞∫Ξ[P(x1∣x1)(A1+S1)−1⋅G¯22⋅(A2+S2)−1⋅G¯12dΞdr]:ε(1)∗0(39)

where Ξ is the spherical surface of radius r. We can also prove that

∫Ξ[(A1+S1)−1⋅G¯21]dΞ=0(40)

where ∫Ξ2H1(n)dΞ2=0; ∫Ξ2H2(n)dΞ2=0.

Additionally, we can prove that

∫ΞninjdΞ=4πr23δij;∫ΞninjnknldΞ=4πr215(δijδkl+δikδjl+δilδjk)(41)

in which, n=r/r is the normal vector at a point on Ξ with r=x2−x1.

Likewise, after substituting Eq. (33) into Eq. (36), ⟨d¯(2)∗⟩(x2) can be stated as

〈d¯(2)∗〉(x2)=∫V−Ω2d¯(2)∗(x2−x1)P(x1∣x2)dx1     =∫V−Ω2[(A2+S2)−1⋅G¯11⋅(A1+S1)−1⋅G¯21:ε(2)∗0]P(x1∣x2)dx1      −∫V−Ω2[(A2+S2)−1⋅G¯12:ε(1)∗0]P(x1∣x2)dx1(42)

〈d¯(2)∗〉(x2)=∫(a1+a2)∞∫Ξ[P(x1∣x2)(A2+S2)−1⋅G¯11⋅(A1+S1)−1⋅G¯21dΞdr]:ε(2)∗0      −∫(a1+a2)∞∫Ξ[P(x1∣x2)(A2+S2)−1⋅G¯12dΞdr]:ε(1)∗0    =∫(a1+a2)∞∫Ξ[P(x1∣x2)(A2+S2)−1⋅G¯11⋅(A1+S1)−1⋅G¯21dΞdr]:ε(2)∗0      +∫(2a2)∞∫Ξ[P(x2∣x2)(A2+S2)−1⋅G¯11⋅(A1+S1)−1⋅G¯21dΞdr]:ε(2)∗0(43)

Similarly, we can also prove that

∫Ξ[(A2+S2)−1⋅G¯12]dΞ=0(44)

In the following sections, Formulation II (higher order) and Formulation I (lower order) will be presented to obtain the effective elastic properties of 3-phase particle reinforced composites. In this paper, Superscript II indicates Formulation II and Superscript I indicates Formulation I. It is noted that Formulation II represents higher-order formulation than Formulation I. In addition, with same elastic material property and same shape of spherical particles in the 1st phase and the 2nd phase, A2=A1, S2=S1 are adopted.

3.1 Formulation II

With lengthy algebra together with employing identities Eqs. (18), (37) and (41), ⟨ε¯(1)∗⟩ is the ensemble volume averaged eigenstrain tensor and is approximately obtained as

⟨ε¯(1)∗⟩UII=Γ(1)UII:ε(1)∗0(45)

where the isotropic tensor Γ∼(1)UII is

Γ∼(1)UII=γ11UIIδijδkl+γ21UII(δikδjl+δilδjk)(46)

γ11UII=ϕ(2)t1UII+ϕ(1)t3UII;γ21UII=12+ϕ(2)t2UII+ϕ(1)t4UII(47)

t1UII=n1×λ3(1+λ)3+125β12×λ5+λ3(1+λ)5×m12+125β12×λ3(1+λ)5×m13+135β12×λ5+λ3(1+λ)7×m14(48)

t3UII=n3+125β12×116×m12+125β12×132×m13+135β12×164×m14(49)

t2UII=n2×λ3(1+λ)3+125β12×λ5+λ3(1+λ)5×m22+125β12×λ3(1+λ)5×m23+135β12×λ5+λ3(1+λ)7×m24(50)

t4UII=n4+125β12×116×m22+125β12×132×m23+135β12×164×m24(51)

n1=m1115β12;n2=m2115β12;n3=n18;n4=n28;λ=a2a1(52)

m11=−150{2β1[2+ν0(−2+5ν0)]+α1[10+ν0(−10+7ν0)]}3α1+2β1(53)

m21=75{2β1[11+ν0(−11+5ν0)]+3α1[10+ν0(−10+7ν0)]}3α1+2β1(54)