## Abstract

Elastic multi-phase materials with a phase having appropriately tuned non-positive-definite elastic moduli have been shown theoretically to permit extreme increases in multiple desirable material properties. Stability analyses of such composites were only recently initiated. Here, we provide a thorough stability analysis for general composites when one phase violates positive-definiteness. We first investigate the dynamic deformation modes leading to instability in the fundamental two-phase solids of a coated cylinder (two dimensions) and a coated sphere (three dimensions), from which we derive *closed-form analytical* sufficient stability conditions for the full range of coating thicknesses. Next, we apply the energy method to derive a general correlation between composite stability limit and composite bulk modulus that enables determination of *closed-form analytical* sufficient stability conditions for arbitrary multi-phase materials by employing effective modulus formulas coupled with a numerical finite-element stability analysis. We demonstrate and confirm this new approach by applying it to (i) the two basic two-phase solids already analysed dynamically; and (ii) a more geometrically complex matrix/distributed-inclusions composite. The specific new analytical stability results, and new methods presented, provide a basis for creation of novel, stable composite materials.

## 1. Introduction

Materials with counterintuitive physical properties have recently attracted increasing attention, resulting in novel solids having negative Poisson's ratio (Lakes 1987), negative thermal expansion (Mary *et al.* 1996), negative refractive index (Shelby *et al.* 2001), negative effective mass density (Liu *et al.* 2005) and negative elastic stiffness in a constrained composite phase (Jaglinski *et al.* 2007).

Composite materials have long offered a convenient approach to create materials with enhanced performance by combining several materials on various length scales. The overall properties of a composite result from the individual specifics of the constituent materials (often called *phases*), their shape, geometrical arrangement and bonding (Milton 2001). While classical composite materials theory contemplated combining materials all having ‘positive’ properties, and the bounds thereby obtained provided perhaps unsurprising limits on the achievable composite response, the possibility of combining materials having positive and negative properties offers the tantalizing prospect of dramatically improved composite behaviour that can far exceed the classical bounds by relaxing an assumption on which all these bounds rely. Indeed, classical composite materials formulas (as well as exact elasticity solutions) for the overall response of a composite with an appropriately tuned negative-property phase predict a response far exceeding the classical bounds, which were long regarded as the absolute limits of the achievable. A prime example is a linear elastic solid, whose elastic moduli are constrained by upper and lower bounds in terms of, for conventional materials, the constituents' moduli and volume fractions (Voigt 1889; Reuss 1929; Paul 1960; Hill 1963; Hashin & Shtrikman 1963). When incorporating one phase having appropriately tuned non-positive-definite elastic moduli (i.e. *negative stiffness*), however, the classical bounds can be dramatically surpassed (Lakes & Drugan 2002). (Models incorporating viscoelasticity and finite deformation confirm the bound-exceeding performance.) Multiple additional composite properties can achieve extreme values via a tuned negative-stiffness phase, e.g. damping capacity, thermal expansion, piezoelectricity and pyroelectricity (Lakes 2001; Wang & Lakes 2001, 2004, 2006).

Violation of positive-definiteness of the elastic moduli raises stability issues. It is well known that a homogeneous, unconstrained, linear elastic body must have positive-definite elastic moduli to be stable overall (Kelvin 1888). However, Drugan (2007) and Kochmann & Drugan (2009) recently proved that the geometrical constraint imposed by an encapsulating positive-definite phase in a composite can provide sufficient stabilization to permit violation of positive-definiteness in the composite's encapsulated phase, while the composite remains stable overall. The derivation of the stability conditions for composites is in general quite complex, requiring numerical evaluation except for very special cases.

Thus, the main goal of this paper is to introduce methods for deriving *closed-form analytical* sufficient conditions for the stability of two-phase solids and general composite materials, and to employ these to derive new specific results for several important basic composite geometries.

## 2. Incremental stability conditions for elastic solids

The state of strain in a body experiencing infinitesimal displacement gradients is characterized by the symmetric infinitesimal strain tensor , where ** u** is the displacement vector field and a comma in a subscript denotes differentiation with respect to the ensuing coordinate(s). For a linear elastic solid, the symmetric stress tensor

**is related to the strain tensor via Hooke's law, , being the fourth-order tensor of elastic moduli. If the solid is also isotropic, the modulus tensor contains only two independent scalar moduli, e.g. the shear and bulk moduli,**

*σ**μ*and

*κ*, respectively; or Young's modulus

*E*and Poisson's ratio

*ν*; or the Lamé moduli

*λ*and

*μ*. For the last case, so that

*σ*

_{ij}=2

*με*

_{ij}+

*λε*

_{kk}

*δ*

_{ij}, where

*δ*

_{ij}is the Kronecker delta. The strain energy density is .

The stability of an elastic body can be investigated by several different methods yielding different conditions in general. First, it is essential to differentiate between material and structural conditions of stability of a macroscopic solid. The *conditions of material stability* (Mandel 1966) are local in nature and derived from the governing (pointwise) differential equations of the elastic medium by requiring that acceleration waves can propagate at finite speeds. Therefore, the material conditions of stability ensure real-valuedness of the speed of travelling waves, which requires strong ellipticity of the stiffness tensor (Kelvin 1888) and positive-definiteness of the acoustic tensor ** T**(

**) with components in every direction**

*n***(Mandel 1962). They also imply rank-one convexity of the strain energy density. For an isotropic, linear elastic material, the conditions of material stability reduce to**

*n**μ*>0 and

*λ*+2

*μ*>0.

The *conditions of structural stability*, in contrast, are non-local and ensure stability of the overall body. Satisfying these guarantees that an arbitrary infinitesimal perturbation of the displacement field from an equilibrium state remains infinitesimal for all time. (It is important to note that this only applies to infinitesimal perturbations; therefore, one may speak of incremental stability in this context. For a sufficiently large imposed perturbation, the system may become unstable and deform into the next stable state. However, as we focus on linear elasticity, the condition of incremental stability is sufficient for our purposes.) The conditions of structural stability are derived from either energetic considerations, enforcing uniqueness of solutions or quasi-convexity of the energy (Kirchhoff 1859; Morrey 1952; Pearson 1956; Hill 1957), or from a dynamic approach that seeks to analyse the eigenmodes of a free vibration with respect to stability in the sense of Lyapunov (1966). The energy method yields the conditions of so-called *static stability*, whereas the dynamic approach determines *kinematic stability* (Ziegler 1953). While both methods provide identical results for conservative systems (Koiter 1965) (such as the linear elastic composites studied in the sequel), they differ for non-conservative systems (e.g. if circulatory or gyroscopic forces are present; cf. Kochmann & Drugan 2011), where only the dynamic approach yields the correct structural stability conditions (Ziegler 1953). For an externally unconstrained homogeneous solid, the conditions of structural stability require positive-definiteness of the fourth-order stiffness tensor, which, for an isotropic, linear elastic medium reduces to *μ*>0 and *κ*>0. Owing to the different forms of the effective bulk modulus *κ*, these conditions are interpreted differently in terms of the Lamé moduli in two and three dimensions: *μ*>0 and must hold in three dimensions, while *μ*>0 and *λ*+*μ*>0 are the corresponding conditions in two dimensions (i.e. plane strain).

For a macroscopic solid without microstructure, the local conditions of material stability are necessary, and the non-local conditions of structural stability are sufficient for overall stability. The situation becomes more complex when studying stability in microstructured solids with randomized or periodic microstructure where cooperative phenomena may arise. Here, the two conditions guarantee periodic (local or short wavelength) stability and long-range (or long wavelength) stability, respectively, and both conditions must be checked independently to ensure overall stability. In particular, in a composite with a distinct separation of macro- and micro-length scales, structural stability on the microscale corresponds to material stability on the macroscale.

Sufficient stability conditions for complex composite geometries are generally difficult to determine analytically, necessitating numerical techniques. Kochmann (in press) showed how sufficient (static) conditions of stability for composites with negative-stiffness phases can be obtained in a straightforward manner from a modal decomposition and eigenvalue analysis based on a finite-element formulation. As shown there, it is a simple exercise to verify that the sufficient condition of stability for an elastic composite (i.e. the condition of structural stability) translates into the requirement of positive-definiteness of the overall stiffness matrix ** K** of the finite-element model. For the separable displacement field of the form

**(**

*u***,**

*x**t*)=

**(**

*v***)**

*x**e*

^{iω t}, this can be verified by computing Rayleigh's quotient, which provides an upper bound on the (squared) lowest eigenfrequency

*ω*

_{0}of the body: 2.1where

**is the mass matrix and**

*M***the vector of all nodal displacements. Equality holds if and only if**

*v***is the eigenvector corresponding to the lowest mode; otherwise, the quotient provides an upper bound for arbitrary ,**

*v***≠**

*v***0**. (Depending on the number of modes of rigid body motion, the finite-element model will contain that number of zero-energy modes and hence zero eigenfrequencies; these do not affect the stability analysis.) For a linear elastic solid, stability requires all eigenfrequencies to be real-valued, i.e. for all

*i*. A complex eigenfrequency (or its complex conjugate which is also a solution) would give rise to exponential growth of displacements, hence instability. (For general solids,

*ω*may be complex and stability requires that Im(

*ω*)>0.) Therefore, the stability of an elastic solid requires that or . As

**is positive-definite in general, it follows that the numerator in (2.1) must be positive for stability; ergo, the stiffness matrix**

*M***must be positive-definite for overall stability. Thus, the requirements of overall stability coincide with the conditions that guarantee positive-definiteness of the stiffness matrix. The numerical solution is mesh-dependent, in general, but shows rapid convergence towards the exact stability limit (Kochmann in press). For all results presented in subsequent sections (although not shown explicitly for conciseness), sufficient accuracy is ensured by repeating simulations with increasing uniform mesh refinement and confirming convergence.**

*K*## 3. Kinematic stability conditions for simple composites

### (a) Modal decomposition

The stability of elastic composites with phases having non-positive-definite moduli has thus far been investigated for geometrically simple two-phase composites. Using an energy method, Drugan (2007) analysed composites consisting of a coated, infinitely extended circular cylinder (two-dimensional plane strain), and a coated sphere (three dimensions); in both cases, the inclusion was assumed to violate positive-definiteness and the (positive-definite) coating was assumed thin compared with the inclusion radius. Kochmann & Drugan (2009) performed a dynamic analysis of the same two-dimension-coated cylinder composite for the full range of coating thicknesses. In that analysis, the displacement field ** u** of the circular cylinder was computed for free vibrations in plane strain, yielding the solution in the modal form (in polar coordinates

*r*and

*φ*; taking the real part of the right side being implied) 3.1where

*t*is time, ,

*ω*

_{mn}the eigenfrequencies and

*v*_{mn}the corresponding eigenmodes. The results show that each mode

*m*provides independent stability conditions, so that overall stability is determined by the most restrictive of all

*m*-modes. Figure 1 illustrates sample numerical results for the stability limits corresponding to different modes

*m*for a coated cylinder (figure 3

*a*, ratio of coating thickness to inclusion radius

*t*/

*a*=3%, ratio of inclusion/coating shear moduli

*μ*

^{I}/

*μ*

^{II}=1/15). The stable regime of positive-definiteness of both phases is considerably expanded owing to the geometric constraint provided by the coating, thus permitting a violation of positive-definiteness in the inclusion while the composite remains stable overall. The strongest restriction on the elastic moduli for overall composite stability is given by the rotationally symmetric

*m*=0 deformation mode; this was found to be true in all cases of positive-definite coating and non-positive-definite inclusion.

A similar full dynamic analysis has not been performed for the three-dimensional-coated sphere. However, it is reasonable to assume, based on the coated cylinder full analysis by Kochmann & Drugan (2009), that the rotationally symmetric mode is also the critical one in causing first instability of the coated sphere, so that only this deformation mode must be analysed to derive the sufficient conditions of stability. We will make this assumption to perform a stability analysis for the full range of coating thicknesses of the coated sphere; the veracity of the assumption will be confirmed in multiple ways. We will first derive the general displacement field representation and then determine the stability limit for a homogeneous sphere (to illustrate and confirm the approach), and then treat the coated sphere. An analogous procedure for the coated cylinder shown in figure 3*a* is summarized in appendix A, as it represents a great simplification to the full analysis of Kochmann & Drugan (2009), and it provides closed-form analytical stability conditions.

### (b) Rotationally symmetric stability analysis for a homogeneous sphere

For the dynamic stability analysis, we consider a homogeneous, isotropic, linear elastic sphere of radius *b*, whose governing differential equations when no body forces act are the equations of motion together with Hooke's constitutive law:
3.2where *ρ* is mass density and superposed dots denote time differentiation. As justified above, we assume the critical deformation mode leading to first instability to be rotationally symmetric, so that the displacement components (*u*_{R},*u*_{θ},*u*_{φ}) in spherical coordinates (*R*,*θ*,*φ*) can be written in the separable form (with the real part of the right side being implied):
3.3Applying (3.3) to (3.2), we arrive at the Navier equations
3.4a
3.4b
and
3.4cThe prime denotes differentiation with respect to radius *R*. As the left-hand sides of (3.4) are independent of *θ* and the equations must hold for all *θ*, the right-hand sides must be zero, requiring *v*_{θ}(*R*)=*v*_{φ}(*R*)=0. This reduces (3.4) to a single ordinary differential equation of Bessel type ((3.4a) with zero right-hand side). Introducing dimensionless radius *r*=*R*/*b* ∈[0,1] with *b* the outer radius, and the abbreviation
3.5we obtain the general solution
3.6Here, *j*_{m}(*x*) and *y*_{m}(*x*) are the order-*m* spherical Bessel functions of the first and second kind, respectively, and *c*_{i} are complex constants. The radial parts of the relevant stress components follow from these via Hooke's law specialized to rotational symmetry:
3.7so that in particular (using dimensionless radius *r*)
3.8It will be useful to employ the following trigonometric representations (Olver 1970):
3.9To avoid a displacement singularity at *r*=0, we must set *c*_{2}=0.

Stability conditions are obtained by applying specific boundary conditions to determine the corresponding infinite set of eigenfrequencies. We enforce zero tractions on the entire sphere surface: *σ*_{rr}(1)=*σ*_{rθ}(1)=*σ*_{rφ}(1)=0. The latter two conditions are satisfied automatically; the eigenfrequencies are determined from *σ*_{rr}(1)=0. Rearranging this (assuming *χ*≠0 and ), we obtain
3.10This transcendental equation yields the infinite set of eigenfrequencies *ω*_{i} from the roots *χ*_{i}. The infinite number of real solutions to (3.10) is illustrated in figure 2*a*.

Instability requires the existence of at least one complex eigenfrequency to result in exponential growth of displacements. Owing to the absence of dissipation in the linear elastic solid, all solutions *ω*_{i} are either real or pure imaginary. If there exists at least one imaginary *ω*_{i}, this will give rise to an imaginary *χ*_{i} (given the necessary conditions of stability are satisfied). Thus, we introduce imaginary roots *χ*=*iy* with and so equation (3.10) becomes
3.11Function *f*(*y*) is illustrated in figure 2*b*. With , it is easily verified that there exist either two symmetric solutions or no solution at all (except *y*=0, which is stable). As the slope of the curves in figure 2*b* vanishes at the origin, the existence of unstable roots is determined from the first non-zero higher-order derivative of *f*(*y*) at the origin. Consequently, unstable roots only exist if
3.12which is equivalent to . Enforcing also the requirement of strong ellipticity (which results from the requirement that for real eigenfrequencies), we obtain the expected stability conditions of positive-definiteness in three dimensions:
3.13If the second condition (positive bulk modulus) is violated, the above analysis shows that the critical deformation mode to cause instability is the lowest, rotationally symmetric mode. This provides one confirmation of our approach. The same conclusion can be drawn for the case of the coated cylinder, presented in appendix A.

### (c) Rotationally symmetric stability analysis for a coated sphere

Next, we derive stability conditions for the composite of figure 3*b*: a homogeneous, isotropic, linear elastic spherical inclusion (radius *a*, moduli *λ*^{I}, *μ*^{I}), embedded in a concentric, homogeneous, isotropic, linear elastic coating of uniform thickness (outer radius *b*, moduli *λ*^{II}, *μ*^{II}). The general displacement field representation in the inclusion is the same as for the homogeneous sphere (employing (3.6) with *c*_{2}=0), while the coating requires employment of the full representation (3.6). Just as for the homogeneous sphere (and the two-dimensional cylinder, cf. appendix A), the angular components of the displacement and traction conditions will not contribute to the stability conditions. Thus, we again need to analyse only the radial components. The radial displacement components in inclusion and coating are (employing (3.9)), respectively,
3.14aand
3.14bwith
3.15and material properties of inclusion (*i*=*I*) and matrix (*i*=*II*). We again define *r*=*R*/*b*∈[0,1] with *b* the outer coating radius, *x*=*a*/*b*, and stresses are again obtained from the displacement fields via Hooke's law. We impose boundary conditions of zero tractions on the entire outer surface and continuity of displacements and tractions across the inclusion/coating interface:
3.16Introducing the dimensionless quantities
3.17the system (3.16) becomes
3.18a
3.18b
and
3.18cWe assume the same mass density for both materials; different densities alter the exact solutions but not the stability results (Kochmann & Drugan 2009; Kochmann in press). Equations (3.18) are written in matrix form as ** M**⋅(

*c*

_{1},

*d*

_{1},

*d*

_{2})

^{T}=

**0**, where

**is the corresponding coefficient matrix; non-trivial solutions require . This is the characteristic equation to be solved for the infinite set of eigenfrequencies**

*M**ω*

_{i}(i.e. one solves for the infinite set of solutions

*j*

_{i}). Employing a numerical search algorithm for the roots of , we numerically determine whether or not imaginary roots

*ω*exist, and we can determine those moduli and radii combinations that allow for imaginary and hence unstable eigenfrequencies, in order to obtain the complete stability map (for positive-definite coating moduli).

A closed-form analytical representation of the stability conditions can be obtained for the case when the inclusion violates positive-definiteness, while the coating remains positive-definite. Stability is lost when at least one imaginary solution *j* exists that solves the characteristic equation . In direct analogy to the homogeneous cylinder discussed in the appendix, such an imaginary solution will exist only if the slope of the characteristic equation is non-zero at the origin; thus, the boundary between stable and unstable regimes is found from the requirement that the characteristic equation's slope vanishes at the origin, i.e. with ,
3.19Let us introduce the parametrization , where *β*∈[0,1] is a dimensionless parameter (Drugan 2007). *β*=0 corresponds to positive-definiteness and *β*=1 to strong ellipticity, so *β*∈[0,1] covers the full range of interest. One can show that (3.19) is satisfied when
3.20which is the analytical expression of the stability limiting curve (i.e. replacing = by > ensures stability).

This result is valid for arbitrary choices of the radii *a* and *b*, i.e. it generalizes to arbitrary coating thicknesses the results of Drugan (2007), which relied on the thin-coating assumption. Figure 4 summarizes the stability conditions for the example case of equal shear moduli and for various relative coating thicknesses. The results obtained from solving the characteristic equation numerically and from (3.20) are identical. Further confirmation of the veracity of these new results for all coating thicknesses is provided by a comparison with results from a finite-element analysis. To this end, figure 4 also includes results obtained from a numerical eigenvalue analysis to determine the static overall stability conditions of the composite sphere (Kochmann in press). It is clear that, within the tolerance of numerical errors, the derived results are confirmed by the finite-element procedure.

## 4. Static stability conditions for composites

### (a) Stability conditions from overall elastic moduli

Thus far, we have shown for two specific composite geometries, via detailed dynamic analyses, that an inclusion material violating positive-definiteness can be stabilized by the surrounding matrix material. We were fortunately able, in these basic geometries, to extract closed-form analytical statements of the stability requirements. For geometrically more complex composites, the stability analysis method of choice would appear to be the numerical finite-element approach summarized at the end of §2. This alone would, of course, provide only numerical stability results. We now introduce a method to determine *closed-form analytical* composite stability requirements for more geometrically complex composite materials, to be used in conjunction with the numerical finite-element approach. This is accomplished by deriving the relationship between overall stability and overall stiffness.

The bulk modulus *κ* of an elastic medium, a scalar measure of the body's resistance to volumetric changes, is defined in general by
4.1where *p* denotes the outward hydrostatic pressure acting on the surface of a body of volume *V* . For a composite of volume *V* with stress ** σ** and strain

**that vary in its interior, one introduces the volumetric averages (**

*ε***denoting position) 4.2so that the bulk modulus of an infinitesimally strained composite can be written as (with**

*X**d*the dimension, and tr

**=**

*T**T*

_{ii}the trace of a tensor) 4.3Note that for anisotropic solids this defines a structural bulk modulus, whereas for overall isotropic solids, this definition yields the effective linear elastic bulk modulus of the composite. Recall that, with Δ

*V*being the change in volume, 4.4With Gauss' theorem and equilibrium without body forces (div

**=**

*σ***0**), we transform the integral for the average stress into a surface integral over ∂

*V*: 4.5with

**=**

*t*

*σ***denoting the surface traction. When a uniform outward hydrostatic pressure**

*n**p*acts on the entire surface,

**=(**

*t**p*

**)**

*I***=**

*n**p*

**on the surface. Then, 4.6Hence, the overall effective bulk modulus of the body subject to a constant**

*n**p*is 4.7

Next, we link the bulk modulus to the stability of the body. The elastic energy stored in a deformed, linear elastic solid is
4.8As we assume no body forces act, the total work *P* performed on the body is due solely to the constant, outward external pressure *p*:
4.9When the deformed body under hydrostatic pressure is in equilibrium, then and consequently . Here, we assume the unstrained body (** ε**=

**0**) to be in equilibrium, so that the constant must vanish. Therefore, we conclude that 4.10Stability of a deformed linear elastic solid subject to infinitesimal strain and stress requires that the stored energy be positive for any non-zero kinematically admissible strain field

**(Pearson 1956; Hill 1957). Applied to the equilibrium configuration, stability requires 4.11For a constant pressure**

*ε**p*, stability hence implies that

*p*and the volume change must always be of the same sign. Therefore, it follows from (4.7) that a linear elastic solid in stable equilibrium must exhibit a positive overall (effective) bulk modulus: 4.12

Our reasoning is sufficiently general for any linear elastic solid: no assumptions were made about isotropy or any other specific form of the elastic modulus tensor, nor about the geometry of the body (except that it be continuous and simply connected). However, this result is not bijective: a positive bulk modulus does not generally guarantee stability for the following reason. The first inequality in (4.11) must hold for arbitrary non-zero kinematically admissible ** ε** for overall stability; a positive bulk modulus, however, only ensures that the integral is positive for the equilibrium configuration and not for any perturbation in general. What this means is that, if a rigorous (e.g. numerical) stability analysis shows a heterogeneous linear elastic body to be in stable equilibrium,

*analytical*requirements for this stability can be obtained from the condition that the body's overall bulk modulus be positive. We will demonstrate this concept in §5, where we derive analytical sufficient stability conditions for example composites and confirm their validity. Furthermore, we conclude that (4.12) is

*necessary*for overall composite stability; our proposed approach for geometrically complex composites is to verify its sufficiency by application of a numerical finite-element approach. We will discuss the implications for general composites in §4

*b*and illustrate the concept explicitly for a periodic matrix/inclusion composite in §5

*c*.

We note in passing that application of (4.12) to the simple specific cases of negative-stiffness-phase composites analysed here shows that, in these cases, as the inclusion bulk modulus is tuned to increasingly negative values, the overall composite will go unstable before infinite composite bulk modulus can be attained. Whether this is true for general composites is still an open question; our research in progress shows that other material properties *can* attain positive extreme values owing to negative inclusion bulk modulus while retaining overall (composite) stability.

For completeness, we also investigate the effect of pre-stress on composite stability. For any infinitesimal pre-stress *σ*_{0}, the above energy method can be modified to yield
and
so that (4.10) remains unaltered and the final conclusions are the same.

### (b) Stability conditions from effective elastic moduli

The geometrical and topological details of actual composite materials are of enormous complexity and often span multiple length scales. For most composites in engineering applications, there exists a considerable distinction between the body's macroscale (i.e. the length scale of the body's outer extensions or its geometrical shape) and its microscale (i.e. the length scale of all microstructural characteristics, e.g. the average diameter of inclusions in a particle-matrix composite). As a consequence, it is legitimate to assume a clear separation of those scales and to characterize the macroscopic material behaviour by effective properties and constitutive laws to be determined from the microscale (e.g. via assuming the existence of a representative volume element of a periodic microstructure and applying the theory of homogenization to obtain the effective mechanical response). We now extend the concept outlined above to general composites and link the sufficient conditions of overall stability to the effective properties of elastic media with microstructure.

We follow the classical theory of composite materials and assume a clear separation of scales between the solid's macroscale and its microstructure, cf. figure 5*a*. Assuming furthermore the same (representative) microstructure to exist throughout the composite, the macroscale composite can be treated as a homogeneous solid with effective properties to be determined from the microstructure via appropriate methods (Hill 1972; Hashin 1972). A common technique to obtain the effective response is the application of analytical or computational homogenization methods to a representative volume element. For simplicity, we restrict ourselves to linear elastic media as in previous sections, for which effective elastic moduli are available in analytical form for many special cases.

Based on the separation of scales, we treat the macroscopic composite as homogeneous with constant effective elastic moduli. For a homogeneous medium subjected to pure traction or mixed traction/displacement boundary conditions, overall stability requires positive-definiteness of the (now effective) elasticity tensor, as was shown for homogeneous solids in the previous section. This condition can be enforced pointwise on the constant effective moduli, so that the condition of structural stability requires the positive-definiteness of the effective elasticity tensor . For an overall isotropic linear elastic solid with microstructure, this gives
4.13where and are the effective shear and bulk moduli, respectively. Therefore, knowledge of the effective elastic properties (for overall isotropic composites of arbitrary geometry) is, in principle, sufficient to obtain the conditions of overall stability (but see provisos in the next paragraph). For any composite model, the effective shear modulus will depend on the constituent shear moduli, and the necessary condition of stability requires (pointwise) strong ellipticity of each constituent. Hence, the individual shear moduli must be positive for stability, so the first condition in (4.13) is generally satisfied automatically. The second condition thus becomes the critical stability condition which, again, requires the effective, overall bulk modulus to be positive for stability. This was already found in the previous section for macroscopic composites of arbitrary geometry, cf. equation (4.12). Here, we conclude that the same condition of stability also applies to general composites with microstructure when replacing the overall elastic properties by the effective (homogenized) ones. In multi-cell materials, cooperative phenomena can affect the overall stability of the solid. Therefore, it is imperative to check stability by two criteria: (i) the local (pointwise) material stability (which only detects periodic phenomena on the level of the representative volume element) and (ii) the structural (global) stability of the structure or solid (which detects long-range instabilities or localization phenomena). Therefore, we must also check macroscopic stability, i.e. we interpret (4.13) as *necessary* conditions of stability that can be employed to derive closed-form analytical stability conditions whose *sufficiency* is to be verified by, e.g. a finite-element study.

Note that this approach only applies to composites whose entire surface phase has positive-definite elastic moduli. In the cases of a non-positive-definite coating, this coating material (comprising the outer free surface of the composite) will always cause global instability (this can be seen explicitly in the coated cylinder composite via the general analysis of Kochmann & Drugan 2009).

## 5. Analytical stability conditions from effective bulk moduli

The conclusion that composite stability requires a positive overall bulk modulus is perhaps unsurprising in light of the known *homogeneous* material stability requirements, yet it facilitates an elegant and alternative derivation of the static composite stability conditions *in a closed form*. We now confirm the veracity of this approach via three specific examples. First, we employ (4.12) to derive closed-form conditions of static stability for the examples of coated cylinder and sphere directly from the overall bulk modulus. Then, we use (4.13) to derive an analytical sufficient condition of static stability for a geometrically more complex matrix/inclusion composite.

### (a) Coated cylinder

The overall bulk modulus of the coated cylinder shown in figure 3*a* is
5.1with *x*=*a*/*b* and *μ*^{i}, *λ*^{i} the elastic moduli of inclusion (*I*) and coating (*II*). Stability conditions follow from requiring , which admits two possible solutions; the one giving the correct stability requirement is
5.2The other solution is discarded because, when violating positive-definiteness, the overall bulk modulus first becomes negative, which is the stability limit sought; it later returns to positive values, to which the discarded inequality corresponds.

The necessary conditions of pointwise stability require *μ*^{I},*μ*^{II}≥0. The composite is stable overall if both phases are positive-definite, i.e. if *λ*^{i}+*μ*^{i}>0 for both materials. The range of interest hence reduces to those cases where the inclusion violates positive-definiteness but still satisfies the necessary conditions of strong ellipticity. Therefore, it is convenient to make use of the dimensionless parameter *β* to define *λ*^{I}=−(1+*β*)*μ*^{I} (Drugan 2007): *β*=0 corresponds to positive-definiteness and *β*=1 to strong ellipticity, so *β*∈[0,1] covers the full range of interest. Inserting *λ*^{I}=−(1+*β*)*μ*^{I} into (5.2) along with the necessary conditions of stability and *x*≤1, we find the following sufficient conditions of overall stability:
5.3Equations (5.3) give new sufficient conditions for stability of the coated cylinder for all coating thicknesses in a closed form; they agree exactly with those determined numerically by Kochmann & Drugan (2009). Furthermore, they reduce exactly to those derived by Drugan (2007) for a thin coating by defining the coating thickness *t* as *b*=*a*+*t*, so that *x*=*a*/(*a*+*t*) and assuming *t*/*a*≪1. The complete set of sufficient conditions in dimensionless form both for the full range of coating thicknesses and for the thin-coating approximation is
Appendix A shows the dynamic approach yields identical conditions.

The present analysis also permits derivation of the sufficient stability conditions for a dilute particle/matrix composite by investigating the limit of large coating thicknesses *b*: taking the limit of in the above solution yields
5.4The fourth condition requires the coating merely to be positive-definite, whereas the third condition allows the inclusion to violate positive-definiteness, with the stable regime for the inclusion moduli increasing with increasing shear stiffness of the coating material. In particular, for *μ*^{II}≥*μ*^{I}, the stability conditions for the inclusion reduce to those of strong ellipticity.

### (b) Coated sphere

The same method is applied to a coated sphere, whose overall bulk modulus is The complete set of sufficient conditions for the full range of coating thicknesses and for the thin-coating approximation is

These results give the sufficient conditions for static stability of the coated sphere for all coating thicknesses in a closed form. They agree exactly with (3.20) determined from the dynamic analysis outlined in §3(c). They also coincide with those of Drugan (2007) in the thin coating limit (*t*/*a*≪1), except for the lower bound on *μ*^{II}/*μ*^{I}, for which Drugan (2007) derived the stronger restriction *μ*^{II}/*μ*^{I}≥*β* *a*/*t*. As the present approach is valid for the full range of coating thicknesses, the above results generalize those of Drugan (2007).

The sufficient stability conditions for the dilute composite in three dimensions are obtained by taking the limit of in the above solution:
5.5Just as for the two-dimensional problem, the fourth condition requires the coating to be merely positive-definite, while the third condition allows the inclusion to violate positive-definiteness, with the stable regime for the inclusion moduli increasing with increasing shear stiffness of the coating material. For *μ*^{II}≥*μ*^{I}, the stability conditions for the inclusion reduce to those of strong ellipticity.

### (c) Solid with periodic microstructure

As a concrete example to confirm the conclusions of §4*b* for more general composites, we apply our approach to an elastic two-phase composite with equal constituent shear moduli *μ*_{1}=*μ*_{2}=*μ*, bulk moduli *κ*_{1} and *κ*_{2} and volume fractions *v*_{1} and *v*_{2}, for which Hill (1963) derived the following exact results for the effective moduli for an isotropic composite (or one with cubic material symmetry) with arbitrarily shaped constituents:
5.6where *κ*_{R} and *κ*_{V} are the Reuss and Voigt bounds on the bulk modulus, respectively:
5.7Importantly, this result does not depend on a positive-definite strain energy density: it also applies to negative constituent bulk moduli. Application of conditions (4.13) directly gives the prospective (to be confirmed by a numerical finite-element stability analysis) sufficient stability conditions for otherwise arbitrary overall isotropic or cubic two-phase composites whose entire surface borders the positive-definite phase:
5.8

First, it is interesting to compare these sufficient stability conditions to those derived in the previous section for a single coated-sphere composite which, rewritten in terms of shear and bulk moduli with equal sphere and coating bulk modulus *μ*>0 (and *κ*_{2}>0), are:
5.9Since *v*_{1}=(*a*/*b*)^{3}, conditions (5.9) and (5.8) exactly coincide.

We now confirm that (5.8) are indeed sufficient stability conditions for a specific sample of a matrix/distributed-inclusions composite material under zero-traction boundary conditions via application of the numerical finite-element stability analysis described earlier. The specific linear elastic composite consists of a square matrix of positive-definite material containing a regular array of *n*×*n* identical square non-positive-definite inclusions with equal spacing (to ensure cubic symmetry) that are fully encapsulated by the matrix material; see the large figure in figure 5*b*. Numerical results were obtained by investigating the eigenvalues of the full overall stiffness matrix (which corresponds to zero-traction boundary conditions on the entire outer surface of the constructed composite). The sought stability limit is found when the first negative eigenvalue arises. This procedure is repeated for composite samples with increasing numbers of inclusions (specifically, we compute results for *n*=1,2,3,6). To permit a comparison with the analytical conditions (5.8), the shear moduli of inclusion and matrix material are chosen to be equal. Furthermore, note that, for comparison with the two-dimensional (plane strain) finite-element model, Hill's result must be modified by using the plane strain bulk moduli of the composite phases. Figure 6 gives an overview of numerical results obtained from the finite-element method. It shows the lower bound on the (normalized) negative inclusion bulk modulus *κ*_{1} versus the (normalized) positive matrix bulk modulus *κ*_{2} for overall composite stability. With an increasing number of inclusions, the finite-element solutions rapidly approach result (5.8) (shown as a solid line). It is clear that the 6×6-composite already shows convincing agreement with the analytical stability condition in both cases shown (for two different volume fractions). Therefore, the numerical results confirm the applicability and sufficiency of the analytical stability conditions (5.8) for a square matrix/inclusion composite having at least 36 inclusions (so that the Hill formula, derived for a representative volume element away from specimen boundaries, gives a sufficiently accurate expression for the bulk modulus for the finite square composite) and of the analytical approach in general, allowing for the stability limit of the elastic solid with microstructure to be determined from the effective elastic moduli of the macroscopic solid.

## 6. Conclusions

We have presented methods to derive closed-form analytical conditions of static and kinematic stability for elastic composites and multi-phase solids having a negative-stiffness phase, which demonstrate the expanded regime of stability owing to the geometric constraints enforced by appropriately placed positive-definite phases. We have applied these methods to the fundamental composites of coated cylinder and coated sphere, for all coating thicknesses, and thereby confirmed and generalized those stability conditions previously available. Furthermore, an analogous method has been outlined for solids with microstructures and its applicability was illustrated by deriving closed-form analytical stability requirements for a specific matrix/distributed-inclusions composite, which were confirmed by a numerical finite-element stability analysis. Our results provide closed-form analytical stability requirements crucial to the exploration and development of novel, stable negative-stiffness-phase composite materials tuned to exhibit extreme properties.

## Acknowledgements

This research was supported by the National Science Foundation (NSF) under grant DMR-0949254, and the Army Research Office/Defense Advanced Research Projects Agency (ARO/DARPA) under grant 57492-EG-DRP.

## Appendix A. Dynamic stability analysis in plane strain

(*a*) *Rotationally symmetric stability analysis for a homogeneous cylinder*

Analogously to the derivation in §3, the displacement components for rotational symmetry in plane strain can be written in the separable form (where the real part of the right sides is implied)
A1using polar coordinates (*R*,*θ*). Applying this to the dynamic governing equations absent body forces, we arrive at the uncoupled Navier equations for plane strain
A2and
A3Using dimensionless radius *r*=*R*/*b*∈[0,1] with *b* the outer radius of the circular cylinder, the general solution representations for the radial parts of the displacement components are
A4with
A5Here, *J*_{m}(*x*) and *Y* _{m}(*x*) are the order-*m* Bessel functions of the first and second kind, respectively, and *c*_{i} are complex constants. The corresponding stresses are obtained from Hooke's law, which in plane strain reads for rotational symmetry
A6Using the above displacement field and dimensionless radius *r*, the radial parts of the stresses are
A7and
A8To avoid a displacement singularity at the origin in the analysis of the homogeneous cylinder, we must have *c*_{2}=*c*_{4}=0. To determine the eigenmodes, we assume zero tractions on the outer boundary, *σ*_{rr}(1)=*σ*_{rθ}(1)=0, giving the two uncoupled conditions to be solved for the infinite set of eigenfrequencies (for *J*_{0}(*χ*_{p})≠0):
A9The definition of the Bessel function of the first kind shows that (A9)_{2} has an infinite number of real solutions, but no imaginary solutions (i.e. it cannot cause instability). Figure 7*a* illustrates the infinite number of real solutions of (A9)_{1}.

Instability arises if at least one of the eigenfrequencies becomes imaginary to result in exponential growth of displacements. Therefore, rewrite *χ*_{p}=*i* *y* and , so that the critical condition (A9)_{1} becomes
A10where *I*_{m}(*x*) is the order-*m* modified Bessel function of the first kind. This condition can have real roots *y*, as illustrated in figure 7*b*. With , it becomes apparent that there exist either two symmetric roots or no root at all (except *y*=0, which is stable). To determine the condition of stability, we must identify those combinations of the elastic moduli for which roots exist. From the limits at , we infer that roots exist only if *f*′(0)>0, giving the stability condition
A11Along with the necessary conditions, *μ*>0 and *λ*+2*μ*>0 (which are apparent from the definitions of *χ*_{p} and *χ*_{s}) we hence deduce the stability conditions
A12the expected conditions of elastic moduli positive-definiteness in two dimensions.

(*b*) *Rotationally symmetric stability analysis for a coated cylinder*

Based on the previous analysis, we now derive the stability conditions for the coated cylinder of figure 3*a*, which consists of a homogeneous, isotropic, linear elastic inclusion (radius *a*, elastic moduli *λ*^{I} and *μ*^{I}) and a homogeneous, isotropic, linear elastic coating (outer radius *b*, moduli *λ*^{II} and *μ*^{II}). The radial part of the general solution for the inclusion is the same as above (i.e. with the Bessel-Y terms being omitted)
A13while the solution in the coating material requires the full representation
A14with
A15with the respective material properties of inclusion (superscript *I*) and matrix (*II*). As before, we employ the dimensionless radius *r*=*R*/*b*∈[0,1] with *b* the outer radius of the cylindrical body, and *a* the radius of the inclusion. We introduce the dimensionless radius ratio *x*=*a*/*b*. The stresses are obtained from the displacement fields in the same manner as for the homogeneous solid.

We investigate the case of zero tractions on the outer boundary, and we enforce continuity of displacements and of tractions across the interface between inclusion and coating. This gives the following system of six equations:
Note that, as for the homogeneous body, the radial and angular components of the displacements and tractions are independent, i.e. the above set of six equations can be reduced to two independent sets of three equations each (one system of equations involving constants *c*_{1},*d*_{1},*d*_{2} and one set involving *c*_{2},*d*_{3},*d*_{4}). This will be beneficial for the determination of the unstable eigenfrequencies.

For brevity, we use abbreviations (3.17). With these dimensionless variables, the two systems of equations assume the following forms. The first system of equations stems from continuity of *v*_{r} and *σ*_{rr}, and zero tractions *σ*_{rr} on the boundary:
The second system of equations results from continuity of *v*_{θ} and *σ*_{rθ}, and zero tractions *σ*_{rθ} on the boundary:
In matrix form, where *M*_{i} are the corresponding coefficient matrices, we have
A16which indicates that non-trivial solutions exist if
A17This is the characteristic equation to be solved for the infinite set of eigenfrequencies *ω*. One finds that only yields real-valued eigenfrequencies. We numerically determine whether or not imaginary roots *j*=*iy* () exist and construct the complete stability map (for positive-definite matrix moduli). The results are identical to those obtained from the full analysis (Kochmann & Drugan 2009) (figure 8).

Just as for the coated sphere, we can obtain a closed-form analytical expression for the stability limit by defining *λ*^{I}=−(1+*β*)*μ*^{I} and evaluating equation (3.19) for matrix *M*_{1} of the present formulation. This yields the stability limit as
A18which agrees exactly with the numerical solution and in the limit of thin coatings agrees with the thin-coating solution of Drugan (2007) as illustrated in figure 8.

- Received September 9, 2011.
- Accepted February 15, 2012.

- This journal is © 2012 The Royal Society