Spectral Analysis of Inhomogeneities Shows that the Elastic Stiffness of Random Composites Decreases with Increasing Heterogeneity
11 Spectral Analysis of Inhomogeneities Shows that the Elastic Stiffness of Random Composites Decreases with Increasing Heterogeneity
Ehsan Ban
Department of Mechanical, Aerospace, and Nuclear Engineering, Rensselaer Polytechnic Institute, Jonsson Engineering Center, Room 2049, 110 8th Street, Troy, NY 12180; Corresponding author, e-mail: [email protected] . Currently at University of Pennsylvania.
Abstract
Investigation of inhomogeneities has wide applications in different areas of mechanics including the study of composite materials. Here, we analytically study an arbitrarily-shaped isotropic inhomogeneity embedded in a finite-sized heterogeneous medium. By modal decomposition of the influence of the inhomogeneity on the deformation of the composite, an exact relation is presented that determines the variation of effective elastic stiffness caused by the presence of the inhomogeneity. This relation indicates that the effective elastic stiffness of a composite is always a concave function of the elastic modulus of the inhomogeneity, embedded inside the composite. Therefore, as the heterogeneity of elastic random composites increases, the rate of increase in effective stiffness caused by the stiffer constituents is smaller than the rate of its decrease due to the softer constitutions. So, weakly heterogeneous random composites become softer and less conductive with increasing heterogeneity at the same mean of constituent properties. We numerically evaluated the effective properties of about ten thousand composites to empirically support these results, characterize the nonaffinity of the displacement field in random composites, characterize the sample to sample variability of the effective properties and extend the results to the thermal conductivity of composites.
1. Introduction
The study of inhomogeneities has wide applications in mechanics. For instance, in determination of the effects of healthy and cancerous cells on the deformation of soft tissues (Schwarz and Safran, 2013; Mills et al., 2014; Shin et al., 2016) as well as predicting the properties of rubber in the presence of filler particles (Smallwood, 1944; Lopez-Pamies et al., 2013). Such composite materials are ubiquitously found in nature, whereas they can also be engineered to produce materials with desired properties. Stiff, lightweight composites are used as structural materials in cars and airplanes (Matthews and Rawlings, 1999); composites with enhanced thermal conductivity are used to cool electronic devices (Chung, 2010). Composites with given volume fractions of constituents exhibit a range of effective properties depending on the spatial arrangement of the constituent phases. In elasticity, if the strain field is uniform over the composite, the effective elastic stiffness equals the weighted arithmetic mean of the constituent moduli. This case represents the upper bound of the effective properties corresponding to the respective volume fraction of constituents. The lower bound is represented by the harmonic weighted average of the constituent properties (Hill, 1963; Milton, 2002). For any composition, specific geometric arrangements of phases can be found that match the properties corresponding to these bounds, indicating that these bounds are optimal (Christensen, 2005). The effective properties of dilute mixtures can be approximated using the solution for the case of a single inhomogeneity embedded in a homogeneous medium. The inhomogeneity problems have been studied in various transport phenomena including elasticity (Smallwood, 1944; Eshelby, 1957; Zhou et al., 2013) and electromagnetism (Maxwell, 1881). Inhomogeneity problems have also been investigated in coupled transport phenomena (Ru, 2000) as well as in materials with nonlinear constitutive behaviors (Suquet, 1997). At larger volume fractions of inhomogeneities, however, the interactions between inhomogeneities causes the dilute approximations to be inaccurate. This difficulty has been addressed using various approximation methods (Christensen, 2005; Nemat-Nasser and Hori, 2013; Willis, 1981). Exact solutions exist for special cases such as for composites with specific constituent geometries and composites with specific distributions of constituents (Milton, 2002; Torquato, 1991). The estimates for random composites can be improved by accounting for the statistical moments of the spatial distribution of constituents, represented by multipoint correlation functions (Torquato, 2005, 1991). Similar approaches have been taken in studies of weakly heterogeneous composites where the differences between the constituent properties are small. Such methods have been utilized to study the bulk modulus of weakly heterogeneous random composites (Molyneux, J and Beran, M, 1965), evaluate the effective thermal conductivity of polycrystals using weak contrast expansions (Avellaneda and Bruno, 1990), and provide solutions for the effective stiffness and toughness of biomimetic random composites (Dimas et al., 2015a, 2015b). Similar approaches have been proposed to evaluate the effective nonlinear properties of composites (CastaΓ±eda, 1991). In this paper, we first present the analytical spectral decomposition of the influence of an arbitrarily-shaped inhomogeneity on the effective elastic stiffness of a finite-sized elastic medium. The modal relation is illustrated and validated using the example of a spherical inhomogeneity embedded in a volumetrically expanding composite. Next, we used the analytical inhomogeneity results together with a weak contrast series expansion to study the effective properties of random elastic and conductive composites. Finally, numerical calculations are presented that empirically support the theoretical results, characterize the displacement of the random composites and its variability, and extend the presented results to the study of the thermal conductivity of composites.
2. Theory
This section describes an exact relation for the variation of the effective elastic stiffness, πΆ eff , of a finite-sized heterogeneous material caused by the addition of an arbitrarily-shaped isotropic inhomogeneity using the modal decomposition of the influence of the inhomogeneity. The inhomogeneity relation is then validated and illustrated using the example case of a spherical inhomogeneity embedded in a volumetrically expanding composite. The modal relation is used together with a weak contrast expansion to study the change in effective properties of random composites at various degrees of the microstructural heterogeneity of the composite. Uniaxial tests are considered to quantify the elastic stiffness of heterogeneous samples. A cubic heterogeneous elastic sample of dimensions πΏ and volume π = πΏ is considered. The effective elastic stiffness is measured by applying a small constant displacement π in the direction π₯ to one of the cube face, π R , while holding the opposite face fixed in π₯ . Few other degrees of freedom are constrained to eliminate rigid body motions. The application of the boundary conditions produces a distribution of boundary traction π = {π€ , π€ , π€ } . The π₯ component of the traction π. π = π€ can be integrated over π R to evaluate the Throughout the article bold italic symbols refer to vector fields. total boundary reaction . Here, π = {1,0,0} denotes the unit vector in the direction π₯ . The effective small-strain elastic stiffness resisting uniaxial deformation in the 1 direction, can be evaluated as πΆ eff = β« π. π ππ π R /(π΄ π ) =β« π π. π ππ π R /(ππ ) , where π = π /πΏ is the applied normal strain in the π₯ direction, and π΄ is the area of π R at rest. We denote this state as the reference state, state 0, with the boundary traction π (0) and effective elastic stiffness πΆ eff(0) . In this state, the deformation of the sample is described by the displacement field, π (0) (Fig. 1a). Throughout the text, parenthesized superscripts refer to the state at which the deformation π is applied. State 0 (Fig. 1a) refers to the state where the inhomogeneity is absent. At state 1 (Fig. 1b-d) an inhomogeneity with elastic modulus πΈ inh(1) is present, and at state 2 the elastic modulus of the inhomogeneity is πΈ inh(2) . For the two vectors π = {π , π , π } and π = {π , π , π } , the componentwise dot product π. π = π π +π π + π π . Fig. 1. (Color online) Schematic of the operation used to evaluate the variation in effective elastic stiffness caused by the presence of an arbitrarily-shaped inhomogeneity (central gray shape). (a) The reference state, state 0, where the inhomogeneity is absent. Effective elastic stiffness πΆ eff is measured in a uniaxial tension test by displacing the right sample boundary in the π₯ direction by π (b) State 1, where the uniaxial extension test is performed while an inhomogeneity of elastic modulus πΈ inh(1) and Poissonβs ratio π is added. (c) Free body diagram of the inhomogeneity with elastic modulus πΈ inh(1) in state 1. The inhomogeneity is deformed by the traction force π (1) . (d) The inhomogeneity is replaced by the reactions to the boundary traction that the inhomogeneity is experiencing, βπ (1) . The inhomogeneity is a Lipschitz domain Ξ© in β with smooth boundary Ξ . In this section, the influence of an inhomogeneity on effective elastic stiffness is evaluated using a modal decomposition and superposition operations. The effect of the inhomogeneity is replaced by the tractions that it exchanges with the matrix while under the influence of the displacements to the compositeβs boundary. The boundary tractions are decomposed into orthonormal components. Then, the change in boundary traction caused by the application of the force component corresponding to each eigenfunction is evaluated using the reciprocal theorem (Love, 2011). It is found that invariants exist for individual modes, which determine the change in sample deformation as a function of the elastic modulus of the inhomogeneity. The invariants are then used to evaluate the variation of effective elastic stiffness caused by the addition of the inhomogeneity. The result is an exact analytical relation between the variation of effective properties of a finite-sized initially heterogeneous medium and the elastic modulus of the inhomogeneity, showing that effective elastic stiffness is a concave function of the properties of an inhomogeneity embedded inside the composite. This result is used later to demonstrate that in random composites the mean effective elastic stiffness decreases with increasing the composite βs heterogeneity. We begin by considering the case where an inhomogeneity with the elastic modulus πΈ inh(1) is added to the heterogeneous cubic sample (Fig. 1b). If the boundary conditions described in state 0 (Fig. 1a) are applied to this sample, the boundary traction changes to π (1) . The effective elastic stiffness and the displacement field in this case are denoted by πΆ eff(1) and π (1) , respectively. We refer to this state of the system as state 1. If we consider the free body diagram of the inhomogeneity in state 1, its boundary is deformed by the displacement πΜ (1) = {π’Μ , π’Μ , π’Μ } and the traction force π (1) = {π‘ , π‘ , π‘ } . As commonly practiced in small-strain elasticity, we replace the inhomogeneity by the forces that it applies to the surrounding medium, βπ (1) . Using the reciprocal theorem (Love, 2011) between the states 0 and 1 yields β« π π (1) . π ππ π R = β« π π (0) . π ππ π R + β« π (1) . πΜ (0) πΞ Ξ . (1) Therefore, the change in the effective elastic stiffness can be expressed as πΆ eff(1) = πΆ eff(0) + 1ππ β« π (1) . πΜ (0) πΞ Ξ . (2) In the next section, we use the modal decomposition of the force and displacement in the inhomogeneity to rewrite Eq. (2) in terms of πΈ inh . Spectral decompositions have been thoroughly studied in problems of elasticity (Ciarlet, 2013; Lebedev and Vorovich, 2003) and are commonly used in studying the normal modes of elastic vibrations (Thomson, 1996). While, eigenvectors form an orthonormal basis in the spectral form of stiffness matrices for frames and discrete structures (Bathe, 1982), eigenfunctions are present in the spectral analysis of continua. In this section, we use eigenfunctions to analyze the spectral decomposition of the displacement and traction at the boundary of an isotropic inhomogeneity in the absence of body forces (Fig. 1c). We begin by considering the tractions and boundary displacement of an isotropic inhomogeneity with an elastic modulus πΈ inh = 1 and Poissonβs ratio π . We denote the traction of the inhomogeneity with unit elastic modulus by πβ² . Following the work of Eshelby (Eshelby, 1959), the
Greenβs function πΊ ππ can be used to find displacement at the boundary of the inhomogeneity due to the boundary traction as π’Μ π = β« πΊ ππ π‘ πβ² πΞ Ξ . Here, the functions πβ² and πΜ are elements of a Hilbert space. The integral equation relating πβ² and πΜ can be represented by πΜ = π πβ² where π is a linear operator associated with the inner product (πβ² (1) , πβ² (2) ) = β« πβ² (1) . πβ² (2) πΞ Ξ (Reed and Simon, 1981). π is positive because (πβ², π πβ²) = β« πΜ . πβ² πΞ Ξ β₯ 0 , since β« πΜ . πβ² πΞ Ξ is twice the external work expended to deform the inhomogeneity by πΜ , which is stored as elastic strain energy and is non-negative. The reciprocal theorem of elasticity shows that (ππβ² (1) , πβ² (2) ) =(πβ² (1) , π πβ² (2) ) , and therefore π is also self-adjoint. Furthermore, over the compact domain Ξ© the integral representation of π indicates that π is compact (Reed and Simon, 1981). Therefore, the spectral theorem of compact self-adjoint operators can be used to decompose the boundary traction and displacements as πβ² =β πΌ π π πππ=1 and πΜ = β π π πΌ π π πππ=1 (Reed and Simon, 1981). Here, π π and π π are eigenvalues and eigenfunctions, πΌ π = (πβ², π π ) , and π is the number of eigenfunctions, which is not necessarily finite in this case. The basis set, π π functions, can be chosen such that the eigenfunctions are orthonormal, i.e. (π π , π π ) = 0 for π β π and (π π , π π ) = 1 for π = π . π π π = π π π π , and therefore, the positive property of π indicates that (π π , π π π ) = π π β₯ 0 . Next, we consider the change in the orthonormal decomposition when the inhomogeneity βs elastic modulus is changed to πΈ inh at the same Poissonβs ratio . We recall that boundary traction is related to stress as π = π π , where π is the stress tensor, and π denotes the unit outward normal at the boundary. For an isotropic material, Hookβs law can be used to express stress as πΈ inh [ π(1+π)(1β2π) π β. πΜ + (βπΜ + (βπΜ ) T )] where π is the identity tensor. Therefore, π is related to πΜ as π = πΈ inh [ π(1+π)(1β2π) π β. πΜ + (βπΜ + (βπΜ ) T )] π . This relation also shows that if πΈ inh = 1 then πβ² = [ π(1+π)(1β2π) π β. πΜ + (βπΜ + (βπΜ ) T )] π . These equations indicate that any traction π and boundary displacement πΜ for an inhomogeneity with modulus πΈ inh can be expressed as π =πΈ inh πβ² and πΜ , where πβ² and πΜ correspond to an inhomogeneity with πΈ inh = 1 and the same π . Therefore, using the spectral decomposition of πβ² and πΜ , for an inhomogeneity with modulus πΈ inh , π and πΜ can be expressed as π = πΈ inh πβ² =πΈ inh β πΌ π π πππ=1 and πΜ = β π π πΌ π π πππ=1 . In this Section, we use Eq. (2) to find πΆ eff(1) at state 1. If we consider that only the reaction to the traction related to the π βth eigenfunction, βπ π(1) = βπΈ inh(1) πΌ π(1) π π , is present in the stretch tests, then we can rewrite Eq. (1) as β« π π π(1) . π ππ π R = β« π π π(0) . π ππ π R + β« π π(1) . πΜ (0) πΞ Ξ , (3) where π π(1) is the traction at the boundary π R due to the displacement π and traction βπ π(1) . The orthonormal property of the eigenmodes, π π , yields that β« π π(1) . πΜ (0) πΞ Ξ = πΈ inh(1) πΌ π(1) πΌ π(0) π π . Therefore, β« π π π(1) . π ππ π R = β« π π π(0) . π ππ π R + πΈ inh(1) πΌ π(1) πΌ π(0) π π . (4) Now, we consider another state where the inhomogeneity has a modulus πΈ inh(2) at the same Poissonβs ratio, and the reaction to the traction from the π βth mode is applied βπ π(2) = πΈ inh(2) πΌ π(2) π π . We denote this state as state 2. If the reciprocal theorem is applied between states 1 and 2 when we apply the traction from the π βth mode , π π , we obtain β« π π π(1) . π ππ π R + β« π π(2) . πΜ (1) πΞ Ξ = β« π π π(2) . π ππ π R + β« π π(1) . πΜ (2) πΞ Ξ . (5) By using the orthogonality of the modes, we can rewrite Eq. (5) as β« π π π(1) . π ππ π R + πΈ inh(2) πΌ π(2) πΌ π(1) π π = β« π π π(2) . π ππ π R + πΈ inh(1) πΌ π(1) πΌ π(2) π π . (6) If the inhomogeneity is not deformed in the π βth mode in state 0, increasing its stiffness does not alter the effective properties. However, if the element deforms, combining Eqs. (4 and 6) and the counterpart of Eq. (4) corresponding to states 2 and 0, gives (πΌ π(0) β πΌ π(1) )πΈ inh(1) πΌ π(1) = (πΌ π(0) β πΌ π(2) )πΈ inh(2) πΌ π(2) = π π . (7) Equation (7) indicates that the quantity (πΌ π(0) β πΌ π(π) )/(πΈ inh(π) πΌ π(π) ) is invariant of πΈ inh(π) . We denote this invariant by π π . We note that π π = 0 when the deformation of the inhomogeneity is prescribed by the boundary conditions and will not change by the application of π π(1) . It can be noted that π π is invariant of πΈ inh(π) , whereas it depends on the prescribed boundary conditions. For a nonzero π π , by replacing πΌ π(1) in Eq. (4) and using the definition of π π from Eq. (7) we obtain β« π π π(1) . π ππ π R = β« π π π(0) . π ππ π R + π π (πΌ π(0) ) π π + 1/πΈ inh(1) . (8) Equations (3 and 8) indicate that β« π π(1) . πΜ (0) πΞ Ξ = π π (πΌ π(0) ) /(π π + 1/πΈ inh(1) ) . The integral β« π (1) . πΜ (0) πΞ Ξ can be decomposed as β β« π π(1) . πΜ (0) πΞ Ξππ=1 . Therefore, Eq. (2) can be rewritten as πΆ eff(1) = πΆ eff(0) + 1ππ β π π (πΌ π(0) ) π π + 1/πΈ inh(1)ππ=1 , (9) where πΆ eff(0) and πΆ eff(1) are the effective stiffness in states 0 and 1. Equation (9) provides a relation between the change of the effective elastic stiffness of the composite and the variation of elastic modulus of an inhomogeneity. We note that π π is necessarily non-negative since otherwise, it is possible to find a variation of the modulus πΈ inh(1) = β1/π π , that leads to the divergence of πΆ eff(1) , which is physically impossible. Equation (9) indicates that effective elastic modulus is, in general, a concave function of the properties of the inhomogeneity: π πΆ eff ππΈ inh2 β€ 0. (10) The exceptional case where equality holds is described in the Discussion Section. The inhomogeneity relation of Eq. (9) is different from most previous results by the fact that we did not assume a specific shape for the inhomogeneity. Many previous exact theories hold for inhomogeneities with specific geometric shapes, for example, spheres, ellipsoids, or polygons. Furthermore, we did not assume a homogeneous medium for the sample containing the inhomogeneity. Our theoretical derivations hold even if the inhomogeneity is part of a finite-sized composite consisting of phases with highly different properties. This assumption distinguishes our derivations from the results obtained using the T-matrix approaches (Gubernatis and Krumhansl, 1975) which are based on the assumption of small fluctuations of the strain field. Additionally, we assumed a finite-sized domain for the composite. Finally, the condition of the continuity of the medium is not required for our analysis. The reciprocal theorem for frame structures (Maxwell, 1864) can be used to derive Eq. (9) in the case of discrete structures such as cellular solids, networks of fibers (Ban et al., 2016), and metamaterials. The spectral theorems linear algebra and the eigendecomposition of compliance matrices can be used in that case. Special cases of Eq. (9) for inhomogeneities with one and two modes of deformation such as elastic rods and Euler-Bernoulli beams were explored in networks of beams (Ban et al., 2016). To validate the modal approach and apply it to a specific case, we evaluated the change in the effective bulk modulus of a spherical composite that includes a spherical inhomogeneity. The inhomogeneity and the surrounding matrix have elastic moduli πΈ inh and πΈ m and Poisson βs ratio π . The two spheres have radii π i and π o and are centered at the same point π (Fig. 2). The effective bulk modulus was measured by displacing the outer surface of the composite, π o , by a small displacement, π , in the radial direction, π π . We then measured the resulting traction π over π o and evaluated the effective bulk modulus as πΎ eff =β« π. π π ππ π o /(π πΏ ) = β« π π. π π ππ π O /(ππΏ ) . Here, the volumetric strain πΏ =3π /π o , and π is the initial surface area of the spherical composite. Fig. 2. (Color online) The illustrative example and the operations used to evaluate the variation of the effective bulk modulus of a spherical elastic composite caused by the presence of a spherical inhomogeneity (gray). (a) The reference state, state 0, where the inhomogeneity is absent. Effective bulk modulus πΎ eff is measured in volumetric expansion tests by the radial displacement of the outer surface of the composite by π . The configurations before and after the application of deformation are plotted using dashed and solid lines, respectively. (b) State 1, where the volumetric expansion test is performed while a spherical inhomogeneity of elastic modulus πΈ inh(1) and Poissonβs ratio π is added to the composite. (c) Free body diagram of the inhomogeneity of elastic modulus πΈ inh(1) in state 1. The inhomogeneity is deformed by the traction force π (1) . We denote the domain and boundary of the inhomogeneity by Ξ© and Ξ . (d) The inhomogeneity is replaced by the reactions to the boundary traction that the inhomogeneity is experiencing, βπ (1) . A planar representation of the traction force is shown for simplicity. Similar to the case of the arbitrarily-shaped inhomogeneity, we consider the two states 0 and 1 where the inhomogeneity is absent (Fig 2a) and where an inhomogeneity of elastic modulus πΈ inh(1) is present, respectively (Fig. 2b). The free-body diagram of the inhomogeneity shows that at state 1 it will have uniform surface traction and displacement of π (1) and πΜ (1) , respectively (Fig. 2c). Without a loss of generality, we assume that the composite is fixed at the point π . We take out few additional degrees of freedom to ensure that the composite does not rotate about π . The deformation of the composite is, therefore, spherically symmetric, and the inhomogeneity only deforms in a single mode, volumetric expansion. The associated surface traction and displacement are π (1) = π‘ π π and πΜ (1) = π’Μ π π . In this case, the modal decomposition, presented in Section 2.2.1, can be employed to express π (1) in terms of an eigenfunction, π (1) =πΈ inh(1) πΌ π , where π is the only eigenfunction. The normal property of the modes requires that the inner product (π , π ) = 1 . Therefore, π = (1/2π i βπ)π π . This eigenfunction represents the volumetric expansion of the inhomogeneity, and the factor (1/2π i βπ)π π ensures that the inner product equals 1. The eigendecomposition also indicates that πΜ (1) = π πΌ π , where π is the only present eigenvalue. Next, the bulk modulus of the inhomogeneity, πΎ inh(1) can be used to relate the magnitude of the traction π‘ and volumetric strain πΏ inh(1) as π‘ = πΎ inh(1) πΏ inh(1) . For an isotropic inhomogeneity πΎ inh(1) = πΈ inh(1) /(3(1 β 2π)) . Therefore, if a displacement π’Μ is present at the boundary of the inhomogeneity, π‘ = πΈπ’Μ /((1 β 2π)π i ) . In this case, comparison with the modal form of π indicates that for a displacement π’Μ , πΌ = 2π’Μ βπ/(1 β 2π) . Noting that in this case, πΜ (1) = π’Μ π π = π πΌ π the eigenvalue can be evaluated as π = π i (1 β 2π). Equation (8) can be in this case rewritten as β« π π . π π ππ π o = β« π π . π π ππ π o + π (πΌ ) π + 1/πΈ inh(1) , (11) where πΌ corresponds to the deformation of the inhomogeneity in state 0. πΌ and the invariant π can be evaluated for this problem using the solution for a pressurized hollow sphere (Bower, 2009). In general, the solution for the displacement in the radial direction takes the form π’(π) = π΄π + π΅/π where π΄ and π΅ depend on the geometry, material properties, and boundary conditions. We begin by πΌ in the absence of internal pressure, or traction in the inner surface of the hollow sphere. By prescribing the displacement π at the outer boundary of the composite and considering the constitutive equation, we evaluate the constants π΄ and π΅ . We then find the displacement at the inner surface of the sphere π’Μ , and evaluate πΌ = 2π’Μ βπ/(1 β 2π) as πΌ = 6π i π o2 βππ (1 β π)(1 β 2π) (2π o3 (1 β 2π) + π i3 (1 + π)) . (12) Similarly, we prescribed a uniform radial traction force of magnitude π‘ at the interior surface of the hollow sphere and evaluated the displacement. We then used the expression π = (π’ β π’ )/(π π‘ ) and evaluated the invariant π as π = (π o3 β π i3 )(1 + π)πΈ m (2π o3 (1 β 2π) + π i3 (1 + π)) . (13) Equations (11 β
13) and the relation between the boundary tractions and πΎ eff can be used to evaluate the change in πΎ eff by the addition of the inhomogeneity as πΎ eff(1) = πΎ eff(0) + 1ππΏ π (πΌ ) π + 1/πΈ inh(1) . (14) Equation (14) demonstrates that πΎ eff is a concave function of πΈ inh(1) . Equation (14) was evaluated for inhomogeneities of different stiffnesses and radii (Fig. 3). The numerical values used to plot Eq. (14) are described in Section 3.2. Fig. 3. (Color online) Effective bulk modulus of a spherical elastic composite containing a spherical inhomogeneity. In all cases, the effective bulk modulus was a concave function of the elastic modulus of the inhomogeneity. The elastic moduli of the inhomogeneity and the surrounding matrix are πΈ inh and πΈ m , respectively. The effective bulk modulus of the composite is plotted for various values of πΈ inh . The effective bulk moduli obtained using numerical calculations, dilute approximation, using Eshelbyβs approach (Christensen, 2005), and the presented spectral method (Eq. 14) are plotted for inhomogeneity over composite radii of (a) 0.2, (b) 0.5, and (c) 0.9. The values of the effective bulk moduli are normalized by the bulk modulus of the corresponding homogeneous solid, where πΈ inh = πΈ m . We used Eq. (10) along with a series expansion to show that random composites become softer and less thermally conductive with increasing heterogeneity. In our analysis composites were made of a large number of inhomogeneities.
We studied composites made by the coalescence of π S homogenous constituents. We denote the properties (elastic modulus or thermal conductivity) of the π 'th constituent by π΄ π . The π΄ π values can be considered as components of a vector, π , which describes the microstructural properties of the composite. We take the case of a homogenous material where all π΄ π βs are the same, as a reference. We refer to the effective properties and the properties of the constituents of the homogeneous material by π΄ eff(0) and π , respectively. In a generic composite, we express the effective properties and microstructural properties by π΄ eff and π , respectively. We can express the effective properties of a composite with respect to the homogeneous sample using a series expansion as π΄ eff (π ) = π΄ eff(0) + [(πΉπ . β π )π΄ eff + 12 (πΉπ . β π ) π΄ eff ] π=π + β―, (15) where πΉπ = π β π . Similar weak contrast series expansions have been used in studying composite systems (Milton, 2002) and random networks of resistors (Luck, 1991). In a random composite, multiple realizations of the composite and hence π , produce a distribution of π΄ eff values. The variance of π΄ eff can be found by calculating the variance of the series truncating the terms beyond the first order in πΉπ , as π π΄ eff β π π΄ π [β π π΄ eff ] π=π . Here, π π΄ π is the variance of the microstructural properties, which grows with increasing composite heterogeneity. Similarly, averaging the series truncated beyond the second order terms in πΉπ , gives π΄ eff(0) β β©π΄ eff βͺ β β(1/2)π π΄ π [ β π2 π΄ eff ] π=π , where β©π΄ eff βͺ denotes the averaged effective properties over an ensemble of composites with the same geometry, but with different realizations of microstructural properties, π . Equation (10)(10) indicates that β π2 π΄ eff is negative in elastic composites. Therefore, β©π΄ eff βͺ < π΄ eff(0) , and the effective elastic stiffness of the composites decrease by increasing heterogeneity. Intuitively, the concave characteristic of the effective stiffness indicates that if the composite is comprised of two types of materials, with increasing the variance of the microstructural properties, the rate of gain in effective elastic stiffness from the stiffer constituents is smaller than the rate of its loss from the softer constituents.
3. Comparison with numerical results
This section presents a comparison of the theoretical results with numerical calculations. First, the numerical methods are presented. Next, the illustrative example of the effective bulk modulus of a composite containing an individual spherical inhomogeneity is described. Finally, two random composites are investigated: composites made by the coalescence of cubic blocks, and composites consisting of spherical inclusions in a homogenous matrix. The composite problems can be expressed as partial differential equations with spatially varying constants. The partial differential equations of elasticity and Fourier heat conduction were solved using the Abaqus (Hibbsett et al., 1998) and COMSOL Multiphysics software packages (COMSOL, 2014). Geometric models and unstructured meshes were generated using Abaqus and COMSOL. The problem domains were discretized into tetrahedral and hexahedral elements. Small enough elements were used to eliminate the mesh sensitivity of the reported results. Small-strain elasticity and steady state Fourier heat conduction were modeled by prescribing small displacement and temperature differences at the boundaries of the samples. The resulting forces and heat fluxes were integrated at the boundaries of the composites and were used for the evaluation of the effective properties. Implicit solution methods were used, and geometric nonlinearities were not considered.
Numerical calculations were performed to verify the modal relation for the bulk modulus of a composite containing a spherical inhomogeneity. A geometric model was made comprising of two concentric spheres: the inner sphere represented the inhomogeneity, and the outer hollow sphere represented the surrounding matrix. The inhomogeneity and the surrounding matrix had elastic moduli πΈ inh and πΈ m , radii π i and π o , and Poissonβs ratio πΈ inh /πΈ m ratios 10 -12 , 5 and few values between them, and π i /π o values of 0.2, 0.5, and 0.9 were tested. In every test, a small outward radial displacement, π , was prescribed to the outer surface of the composite, and the resulting traction was evaluated in the outer surface of the sample. A radial strain of 0.5% was prescribed in all cases. The effective bulk modulus was evaluated as the ratio of the mean surface traction to the volumetric strain. Small-strain deformation was considered, and the volumetric strain was evaluated as /π o . The effective bulk modulus increased with increasing the stiffness of the inhomogeneity and exhibited a concave curve (Fig. 3). The same trend was observed when testing inhomogeneities of various radii. As expected, the change in effective bulk modulus was larger in composites with larger inhomogeneities. These results are compared with the theoretical results in the Discussion Section. To support the theoretical results on elastic random composites (Section 2.4) and empirically extend them to thermally conductive composites, numerical calculations were performed. Cubic representative volume elements were made by the coalescence of cubic homogeneous constituents of equal size. To demonstrate the effect of heterogeneity on the mean composite properties, the material properties of the constituents were sampled from statistical distributions of mean π΄ eff(0) and variance, π π΄2 , where π π΄2 was varied in calculations. For each value of π π΄2 , an ensemble of five hundred replicas of the composite were tested. Effective properties were measured in each of the five hundred samples. The mean and variance of the resulting distributions of effective properties, β©π΄ eff βͺ and π π΄ eff , were computed in terms of the variance of local properties, π A2 . The effective thermal conductivity was measured by the application of a small temperature difference, βπ (0) , to two opposing faces of the composite. Insulating boundary conditions were prescribed at the other sample boundaries. The effective thermal conductivity, π eff , was evaluated as βππΏ/βπ (0) where π is the heat flux per unit area over one of the surfaces of the composite with prescribed temperature. The numerical calculations indicated that for an arbitrary π΄ eff(0) , the mean of the distribution of composite properties decreases from the reference π΄ eff(0) in proportion to the variance of the microstructural properties: π΄ eff(0) β β©π΄ eff βͺ~π π΄2 (Fig. 4a). In addition, the variance of effective properties increases linearly with the variance of microstructural properties: π π΄ eff ~π π΄2 (Fig. 4b). Figure 4 includes the results for both elastic stiffness and thermal conductivity and for bimodal and lognormal distributions of microstructural properties. We observed that increasing the number of homogeneous constituents in the model leaves the results on composites unchanged. The data sets are described in the figure caption. Fig. 4. (Color online) Application of the inhomogeneity relation to random composites made by the coalescence of cubic blocks shows that effective elastic stiffness and thermal conductivity decrease with increasing heterogeneity in microstructural properties. Numerical results show (a) the reduction of the mean effective properties β©π΄ eff βͺ and (b) increase of the variance of effective properties π π΄ eff with the variance of the microstructural properties, π π΄2 /π π΄2 . The symbols represent data corresponding to the effective elastic stiffness of composites with lognormal (outlined squares) and bimodal (orange triangles) distributions of microstructural properties, and the effective thermal conductivity of composites with bimodal distributions of local thermal conductivity (red circles). Effective properties are plotted for composites that are twice (plus signs) and three times (open diamonds) larger than the original composite. In composites with bimodal distributions of properties, two types of constituents are present that have equal probabilities of presence. The data points and dashed lines correspond to the numerical results and analytical power laws, respectively. We also tested the bulk modulus of the random composites made by the coalescence of cubic blocks (Fig. 5). The elastic moduli of the blocks were sampled from identical bimodal distributions of varying standard deviations. An ensemble of 50 samples with different realizations of microstructural properties were tested at each value of standard deviation. Mean effective bulk moduli were evaluated for each group of composites. The mean bulk modulus of the composites decreased with increasing microstructural heterogeneity, and the amount of decrease was proportional to the variance of the microstructural properties (Fig. 5). The resulting mean effective bulk moduli resided between the upper and lower Voigt-Reuss bounds (Hill, 1963; Milton, 2002) and also the upper and lower Hashin-Shtrikman bounds (Hashin and Shtrikman, 1963) (Fig. 5). Fig. 5.
Comparison of the mean effective bulk modulus of the random composites made by the coalescence of cubic blocks with the Hashin-Shtrikman and Voigt-Reuss bounds. The elastic moduli of the cubic blocks were sampled from identical bimodal distributions. The inset displays a snapshot of a tested representative volume element, where black and white colors represent constituents with two different elastic moduli.
Moreover, we tested random composites made from different numbers of constituents. The numerical results indicated that increasing the number constituents in the composite, leaves the decrease of effective properties unchanged (Fig. 4a). However, the variance of the effective properties is inversely related to the number of microstructural constituents, π S (Fig. 6a). This result is in agreement with the previous observations in elastic composites (Jeulin et al., 2004) and is a consequence of the central limit theorem (Dekking et al., 2007). Since the variance of effective properties becomes smaller as composites become larger, we conclude that in composites made of a large number of constituents, the mean of the effective properties is representative of the effective properties of all of the random composites in the ensemble. In the absence of heterogeneity, strain is spatially uniform over the composite. In contrast, in random composites, the distribution of strain is nonuniform, and the variance of the strain field grows with increasing the heterogeneity of the composites (Fig. 6b). We calculated the variance of the strain field using the metric π = β©π’ xβ²2 βͺ V /π where π’ xβ² = π’ β π’ uniform is the non-uniform component of displacement in the direction of the uniaxially applied macroscopic strain. π’ and π’ uniform are the total and uniform components of the displacement. Here, β©π’ xβ²2 βͺ V denotes the mean of π’ xβ²2 over the compositeβs volume. In conclusion, this result indicates that the growth of heterogeneity in the displacement field with microstructural heterogeneity agrees with the reported softening effect marking the departure of effective properties from the Voigt bound as heterogeneity increases.
Fig. 6. (Color online) In the random composites made by the coalescence of cubic blocks, the variance of effective properties is inversely proportional to the number of homogeneous microstructural constituents, and the heterogeneity in the displacement field grows in proportion to the variance of the microstructural properties. The increase in the variance of displacement marks a departure from a uniform state of strain with increasing heterogeneity. All quantities were normalized by the respective values corresponding to the results of Fig. 4. The data points and dashed lines correspond to the numerical results and analytical power laws, respectively. We also tested the effective elastic stiffness and thermal conductivity of random composites consisting of spherical inclusions randomly placed in homogenous matrices. Spheres of uniform volume and random properties were placed in a cubic representative volume element by the random choice of the Cartesian coordinates of their centers. In different sets of tests, inclusions were modeled that each occupied 0.002 or 0.003 of the compositeβs volume.
The elastic moduli or thermal conductivities of the inclusions were independently sampled from identical bimodal distributions. Both the elastic inclusions and the surrounding matrix had Poissonβs ratio of 0.3 . The mean elastic modulus or thermal conductivity of the inclusions was 20 times that of the surrounding matrix. Tests were performed by prescribing a small displacement or temperature difference at the sample boundaries. The same boundary conditions used in the case of cubic composites were used here, and the effective properties were evaluated using the resulting boundary traction forces or heat fluxes. Ensembles of 50 composites were tested at each variance of inclusion properties, where different realizations of inclusion properties were present in each test. The matrix properties were the same in all tests. The thermally conductive composites became less conductive on average with increasing the variance of the thermal conductivity of the inclusions. At small to moderate variances of the properties of the inclusions, the decrease in the effective thermal conductivity was proportional to the variance of the properties of the inclusions (Fig. 7). Similarly, we tested elastic composites with inclusions that in total occupied 1% and 10% of the sample volume. In the case with 10% volume fraction, π π /π , we tested composites with different numbers and sizes of individual inclusions. The first case had a smaller number of large inclusions, and the second case had more inclusions of a smaller size. In all cases, effective stiffness decreased with increasing the variance of the elastic moduli of the inclusions. Furthermore, at the small and moderate variances of the elastic modulus of the inclusions, the decrease in effective stiffness was proportional to the variance of the moduli of the inclusions (Fig. 7). The effective stiffness of the composites with a higher volume fraction of inclusions decreased more than the composites with a smaller volume fraction of inclusions. At the same total volume fraction, however, the relative decrease of effective stiffness was independent of the number and size of the individual inclusions. Taken together these tests demonstrate that the decrease in effective properties is also present in random composites with spherical inclusions of random properties. Fig. 7. (Color online) The effective elastic stiffness and heat conductivity of random composites with spherical inclusions with random properties decreases in proportion to the variance of the properties of the inclusions. The relative decrease of the effective elastic stiffness depends on the total volume fraction of inclusions, π π /π . But at the same total volume fraction of inclusions, the relative decrease is independent of the number of inclusions and their size. The data points and dashed lines correspond to the numerical results and analytical power laws, respectively.
4. Discussion
To illustrate the modal approach and validate it, we evaluated the change in the effective bulk modulus of a composite containing an individual spherical inhomogeneity. The results showed a close agreement between the numerical calculations and the modal approach (Fig. 3). We also compared the result with a dilute approximation of effective bulk modulus based on the solution for spherical inhomogeneities (Fig. 3, red dashed curves) , derived using Eshelbyβs approach (Christensen, 2005). As expected the dilute approximation agrees with the effective bulk moduli of the composites with small volume fractions of inhomogeneities. In the systems where the inhomogeneity radius is 90% of the composite radius, however, the dilute approximation predicts effective bulk moduli that are smaller than the bulk moduli obtained using numerical calculations and the modal relation (Fig. 3). The difference between the results increases with increasing the difference between the properties of the inhomogeneity and the surrounding matrix. In contrast, the numerical calculations and the modal relation closely agree even in the tests performed at a large volume fraction of the inhomogeneity (Fig. 3). This example demonstrates the validity of the modal derivation in finite-sized composites. The studied composites have constituents whose properties are independently sampled from identical statistical distributions. Our results are, therefore, expected to apply to composites without long-range correlations of microstructural properties. It is expected that if long-range correlations are introduced to bias the sampling of the microstructural properties, the effective elastic stiffness still decreases with increasing heterogeneity. The averaging calculations, however, will no longer be valid, and we expect that the amount of decrease in stiffness will no longer be proportional to the variance of the microstructural properties. Moreover, in cases where a gradient of material properties is present, we expect a decrease of effective stiffness with increasing heterogeneity, at the same mean of properties, whereas we do not expect the amount of decrease to be proportional to the variance of the microstructural properties. Future works may extend the present results to heterogeneous media with spatially correlated microstructural properties. The modal relation shows that the effective stiffness is a concave function of the elastic modulus of an inhomogeneity inside a composite (Eq. (10)). The derivation also points to an exceptional case where the second derivative of effective stiffness with respect to the elastic modulus of the inhomogeneity vanishes. A system of elastic rods (or Hookean springs) connected in parallel, is an example of this exceptional case, where the deformation of the inhomogeneity is prescribed by the boundary conditions, and all π π βs vanish. Therefore, π πΆ eff /ππΈ inh2 = 0 , and the effective stiffness of the system does not change with increasing the variance of the elastic modulus of the rods, at the same mean of constituent properties. In this exceptional system, after external displacement is prescribed to evaluate the effective properties, the individual constituents cannot be further deformed by the application of external forces that correspond to the modes that comprise the deformation of the constituent. The developed modal approach presents insights into the mechanics of inhomogeneities and random composites. Its extension to different classes of composites, however, requires the determination of the invariants π π for specific composite geometries. The π π parameters are invariant of the inhomogeneity βs modulus, whereas they depend on the shape of the inhomogeneity, the response function of the surrounding medium, and the prescribed boundary conditions. Further developments of the modal approach may lead to novel solutions for inhomogeneity problems and composites.
5. Conclusions This paper presented a modal analysis of the influence of an inhomogeneity on the effective elastic stiffness of a composite. The analysis showed that effective stiffness is always a concave function of the elastic modulus of an inhomogeneity inside the composite. The concave property indicates that weakly heterogenous random composites become softer with increasing microstructural heterogeneity; the loss of the stiffness by the softer constituents is larger than the gain from the stiffer constituents. The modal relation was illustrated and validated in the example case of a spherical inhomogeneity in a volumetrically expanding composite. The results on random composites were validated using numerical calculations and were compared with the Voigt-Reuss and Hashin-Shtrikman bounds. Moreover, the numerical results showed that random composites including spherical inclusions of random properties become softer with increasing the variance of the elastic moduli of the inhomogeneities. Finally, the results were numerically extended to thermally conductive composites, indicating a decrease in effective thermal conductivity with increasing microstructural heterogeneity.
Acknowledgments
We thank Javad Heydari (formerly at RPI) and Poorya Mirkhosravi (UCSD) for fruitful discussions. We thank Professor Catalin R. Picu (RPI) for guidance.
Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.
References
Avellaneda, M., Bruno, O., 1990. Effective conductivity and average polarizability of random polycrystals. J. Math. Phys. 31, 2047. https://doi.org/10.1063/1.528656 Ban, E., Barocas, V.H., Shephard, M.S., Picu, R.C., 2016. Softening in random networks of non-identical beams. J. Mech. Phys. Solids 87, 38 β
50. https://doi.org/10.1016/j.jmps.2015.11.001 Bathe, K.J., 1982. Finite Element Procedures. Prentice Hall, New York, NY. Bower, A.F., 2009. Applied Mechanics of Solids. CRC Press. CastaΓ±eda, P.P., 1991. The effective mechanical properties of nonlinear isotropic composites. J. Mech. Phys. Solids 39, 45 β
71. https://doi.org/10.1016/0022-5096(91)90030-R Christensen, R., 2005. Mechanics of Composite Materials. Dover Publications, Mineola, N.Y. Chung, D., 2010. Composite Materials: Science and Applications, 2nd ed. 2010 edition. ed. Springer, London; New York. Ciarlet, P.G., 2013. Linear and Nonlinear Functional Analysis with Applications. SIAM, Philadelphia, PA. COMSOL Multiphysics Reference Manual, version 5, 2014. Dekking, F.M., Kraaikamp, C., LopuhaΓ€, H.P., Meester, L.E., 2007. A Modern Introduction to Probability and Statistics: Understanding Why and How. Springer, London. Dimas, L.S., Veneziano, D., Giesa, T., Buehler, M.J., 2015a. Probability distribution of fracture elongation, strength and toughness of notched rectangular blocks with lognormal Youngβs modulus. J. M ech. Phys. Solids 84, 116 β Properties of Heterogeneous Rectangular Blocks With Lognormal Youngβs
Modulus: Effective Moduli. J. Appl. Mech. 82, 011003 β β β β β Hibbett, Karlsson, Sorensen, 1998. ABAQUS/standard: Userβs Manual. Hibbitt,
Karlsson & Sorensen. Hill, R., 1963. Elastic properties of reinforced solids: Some theoretical principles. J. Mech. Phys. Solids 11, 357 β β
27. https://doi.org/10.1007/978-1-4020-2316-3_5 Lebedev, L., Vorovich, I., n.d. Functional Analysis in Mechanics. Springer, New York. Lopez-Pamies, O., Goudarzi, T., Nakamura, T., 2013. The nonlinear elastic response of suspensions of rigid inclusions in rubber: I β An exact result for dilute suspensions. J. Mech. Phys. Solids 61, 1 β
18. https://doi.org/10.1016/j.jmps.2012.08.010 Love, A.E.H., 2011. A Treatise on the Mathematical Theory of Elasticity, 4th Revised ed. edition. ed. Dover Publications, New York. Luck, J.M., 1991. Conductivity of random resistor networks: An investigation of the accuracy of the effective-medium approximation. Phys. Rev. B 43, 3933 β β β lbyβs problem for twoβ dimensional piezoelectric inclusions of arbitrary shape. Proc. R. Soc. Lond. Math. Phys. Eng. Sci. 456, 1051 β β Shin, B., Gopaul, D., Fienberg, S., Kwon, H.J., 2016. Application of Eshelbyβs
Solution to Elastography for Diagnosis of Breast Cancer. Ultrason. Imaging 38, 115 β β β β
76. https://doi.org/10.1115/1.3119494 Willis, J.R., 1981. Variational and Related Methods for the Overall Properties of Composites. Adv. Appl. Mech. 21, 1 β
78. https://doi.org/10.1016/S0065-2156(08)70330-2 Zhou, K., Hoh, H.J., Wang, X., Keer, L.M., Pang, J.H.L., Song, B., Wang, Q.J., 2013. A review of recent works on inclusions. Mech. Mater. 60, 144 ββ