Oscar Lopez-Pamies

Pedro Ponte Castañeda

Second-Order Estimates for the Macroscopic Response and Loss of Ellipticity in Porous Rubbers at Large Deformations

Keywords: homogenization, finite deformations, elastomeric foams, microstructure evolution

. The reasons for this result have been linked to the evolution of the microstructure, which, under appropriate loading conditions, can induce geometric softening leading to overall loss of ellipticity. Furthermore, the "second-order" homogenization method has the merit that it recovers the exact evolution of the porosity under a finite-deformation history in the limit of incompressible behavior for the matrix.

Introduction

This article is concerned with the use of recently developed homogenization methods to determine the macroscopic behavior of porous elastomers, as well as the associated evolution of their microstructure (e.g., porosity) and the possible development of instabilities, when these materials are subjected to finite deformations. In an early contribution, Blatz and Ko [START_REF] Blatz | Application of finite elastic theory to the deformation of rubbery materials[END_REF] performed various experiments on a polyurethane rubber with a random distribution of voids about 40 μ in diameter and an approximate volume fraction of 50%. The experimental results generated allowed these authors to propose a phenomenological model, known as the Blatz-Ko model, which turns out to give adequate predictions at least for certain types of loadings. A physically appealing property of this constitutive model is that it is found to lose ellipticity at finite deformation histories [START_REF] Knowles | On the ellipticity of the equations of nonlinear elastostatics for a special material[END_REF]. This property is in agreement with experimental evidence suggesting that such materials can develop macroscopic shear band instabilities at sufficiently high deformations, which correspond to buckling of the ligaments between the pores at the micro scale. Motivated by this earlier work, Abeyaratne and Triantafyllidis [START_REF] Abeyaratne | An investigation of localization in a porous elastic material using homogenization theory[END_REF] attempted a numerical study of the overall behavior of a nearly incompressible Neo-Hookean matrix with a periodic distribution of cylindrical pores. They made use of the results of homogenization theory for periodic media [START_REF] Sanchez-Palencia | Comportements local et macroscopique d'un type de milieux physiques heterogenes[END_REF], which allow the reduction of the problem of determining the effective behavior of a composite to a numerical calculation on a unit cell. An interesting and important finding in this work was that the homogenized constitutive model for the porous material loses ellipticity, even when the matrix material itself does not. This pioneering work has been generalized and developed further by Triantafyllidis and coworkers [START_REF] Triantafyllidis | On the comparison between microscopic and macroscopic instability mechanisms in a class of fiber-reinforced composites[END_REF][START_REF] Geymonat | Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity[END_REF], always in the context of periodic hyperelastic composites. One of the main conclusions of this work is that loss of strong ellipticity for the homogenized material can be shown rigorously to correspond to the bifurcation of the composite at wavelengths much larger than the size of the unit cell. Furthermore, the overall loss of strong ellipticity (i.e., the possible emergence of shear bands) in the homogenized composite, which is due to the buckling of the actual material at the microstructural level, provides an upper bound for the stable domain of the composite.

Although a proper homogenization framework has been available for some time for hyperelastic composites with random microstructures [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF], the application of this framework to random porous elastomers has apparently not yet been attempted, presumably because of the technical difficulties associated with this problem. There are, however, some rigorous estimates for special loadings [START_REF] Hashin | Large isotropic elastic deformation of composites and porous media[END_REF], as well as other estimates based on various types of ad hoc approximations, mostly for low-density foams (see, e.g., [START_REF] Gent | The deformation of foamed elastic materials[END_REF][START_REF] Feng | Nonlinear deformation of elastomeric foams[END_REF][START_REF] Gibson | Cellular Solids[END_REF]). Our proposal here is to make use of the "second-order" homogenization method, originally developed by Ponte Castañeda [START_REF] Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations[END_REF] for viscoplastic materials, and extended recently for hyperelastic composites by Lopez-Pamies and Ponte Castañeda [START_REF] Lopez-Pamies | Second-order estimates for the large-deformation response of particle-reinforced rubbers[END_REF][START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF]. For comparison purposes, we will also make use of an earlier version of the method due to Ponte Castañeda and Tiberio [START_REF] Castañeda | A second-order homogenization procedure in finite elasticity and applications to black-filled elastomers[END_REF] (see also [START_REF] Castañeda | Exact second-order estimates for the effective mechanical properties of nonlinear composite materials[END_REF]). The advantage of these methods is that they can be used for any type of composite system, and that they make use of standard estimates for a suitably optimized "linear comparison composites" to generate corresponding estimates for the nonlinear hyperelastic composite. These methods have already been used to estimate the behavior of particle-reinforced elastomers [START_REF] Castañeda | A second-order homogenization procedure in finite elasticity and applications to black-filled elastomers[END_REF][START_REF] Lopez-Pamies | Second-order estimates for the large-deformation response of particle-reinforced rubbers[END_REF][START_REF] Lahellec | Second-order estimate of the macroscopic behavior of periodic hyperelastic composites: Theory and experimental validation[END_REF][START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF], and have been shown to be able to handle the strongly nonlinear constraint of material incompressibility (a constraint on the determinant of the deformation) for these material systems. These encouraging results for particle-reinforced elastomers are suggestive that the methods can also be used successfully for porous elastomers. This being our first attempt to handle porous elastomers in the context of finite deformations, explicit results will be generated only for a model problem involving aligned two-dimensional pores distributed randomly and isotropically in an (in)compressible, isotropic elastomer. It is important to emphasize, however, that the method can be applied for general classes of constitutive models in finite elasticity, as well as for very general classes of microstructures. The aim of this first work is to explore the capabilities of the methods in the context of a simple example, albeit one that incorporates the essential features of the problem including strongly nonlinear behavior for the matrix phase, the possible evolution of the microstructure and its implications for the overall stability of the material, as determined by the strong ellipticity condition.

Preliminaries on Hyperelastic Composites

Consider a material made up of N different (homogeneous) phases distributed randomly in a specimen occupying a volume 0 in the reference configuration.

Here, the characteristic length of the inhomogeneities (e.g., voids) is assumed to be much smaller than the size of the specimen and the scale of variation of the applied loading.

The constitutive behavior of the phases is characterized by stored-energy functions W (r) (r = 1, . . . , N) that are nonconvex functions of the deformation gradient F. The local energy function of the composite may thus be written as

W (X, F) = N r=1 χ (r) (X)W (r) (F), (1)
where the functions χ (r) are equal to 1 if the position vector X is inside phase r (i.e., X ∈ (r) 0 ) and zero otherwise. The stored-energy functions of the phases are, of course, assumed to be objective in the sense that W (r) (QF) = W (r) (F) for all proper orthogonal Q and arbitrary deformation gradients F. Making use of the polar decomposition F = RU, where U is the right stretch tensor and R is the rotation tensor, it follows, in particular, that W (r) (F) = W (r)

(U).

The local or microscopic constitutive relation for the material is then given by

S = ∂W ∂F (X, F), (2)
where S denotes the first Piola-Kirchhoff stress tensor. Note that sufficient smoothness has been assumed for the W (r) on F. Furthermore, the stored-energy functions W (r) will be assumed to be such that W (r) (F) → ∞ as det F → 0+, to ensure the material impenetrability condition: det F(X) > 0 for X in 0 . This condition would be automatically satisfied for incompressible materials, where det F is required to be exactly 1.

Following Hill [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF] and Hill and Rice [START_REF] Hill | Elastic potentials and the structure of inelastic constitutive laws[END_REF], the global or macroscopic constitutive relation for the composite is defined by

S = ∂ W ∂F , (3)
where S = S , F = F are the average stress and average deformation gradient, respectively, and

W (F) = min F∈K(F) W (X, F) = min F∈K(F) N r=1 c (r) W (r) (F) (r) (4)
is the effective stored-energy function of the composite. In the above expressions, the symbols • and • (r) denote volume averages over the composite ( 0 ) and over the phase r ( (r) 0 ), respectively, so that the scalars c (r) = χ (r) represent the (initial) volume fractions of the given phases. Furthermore, K denotes the set of admissible deformation gradients:

K(F) = F | x = χ (X) with F = Gradχ in 0 , x = FX on ∂ 0 . ( 5
)
Note that W physically represents the average elastic energy stored in the composite when subjected to an affine displacement boundary condition. Moreover, from the definition (4) and the objectivity of W (r) , it can be shown that W is objective, namely, W (F) = W (U). Here, U represents the macroscopic right stretch tensor associated with the macroscopic polar decomposition F = R U, with R denoting the macroscopic rotation tensor. In turn, the objectivity of W implies that the macroscopic rotational balance equation S F T = F S T is automatically satisfied [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF][START_REF] Ogden | Extremun principles in non-linear elasticity and their application to composites -I Theory[END_REF].

It is further recalled that since W cannot be convex suitable hypothesis are needed to ensure the existence of minimizers in (4). Ball [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF] has provided sufficient conditions for the existence of such minimizers, including the hypothesis of polyconvexity of W , together with suitable growth conditions for W . More mathematically precise definitions of the effective energy W are available at least for periodic microstructures [START_REF] Braides | Homogenization of some almost periodic coercive functionals[END_REF][START_REF] Müller | Homogenization of nonconvex integral functionals and cellular elastic materials[END_REF]. Such definitions generalize the classical definition of the effective energy for periodic media with convex energies [START_REF] Marcellini | Periodic solutions and homogenization of nonlinear variational problems[END_REF], by allowing for possible interactions between unit cells, essentially by taking an infimum over the set of all possible combinations of units cells. Physically, this corresponds to accounting for the possible development of instabilities in the composite at sufficiently high deformation. In this connection, it is important to remark that Geymonat et al. [START_REF] Geymonat | Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity[END_REF], following earlier work by Triantafyllidis and Maker [START_REF] Triantafyllidis | On the comparison between microscopic and macroscopic instability mechanisms in a class of fiber-reinforced composites[END_REF], have shown rigorously that loss of strong ellipticity in the homogenized behavior of the composite corresponds to the development of long wavelength instabilities in the form of localized shear bands. Furthermore, the "failure surfaces" defined by the macroscopic loss of strong ellipticity condition are actually upper bounds for the onset of other types of instabilities.

The objective of this work is to obtain estimates for the effective stored-energy function W of the above-defined hyperelastic composites subjected to finite deformations, with particular interest in the special case of porous elastomers, where the second phase is vacuous. Motivated by the earlier work of Abeyaratne and Triantafyllidis [START_REF] Abeyaratne | An investigation of localization in a porous elastic material using homogenization theory[END_REF], a second objective will be to investigate whether or not the homogenized behavior of the porous material can lose strong ellipticity, even when the local behavior of the matrix phase is assumed to be strongly elliptic itself. The determination of the effective stored-energy function of a porous elastomer is a very difficult problem, because it amounts to solving a set of highly nonlinear partial differential equations with random coefficients. As a consequence, there are very few analytical estimates for W . Ogden [START_REF] Ogden | Extremun principles in non-linear elasticity and their application to composites -I Theory[END_REF] noted that use of the trial field F = F in the definition (4) for W leads to an upper bound analogous to the wellknown Voigt upper bound in linear elasticity. Due to the well-known difficulties associated with the definition of a complementary energy principle in finite elasticity, the equivalent of a Reuss-type bound in linear elasticity is not straightforward. A non-trivial lower bound that used only information on the phase volume fractions was proposed by Ponte Castañeda [START_REF] Castañeda | The overall constitutive behaviour of nonlinearly elastic composites[END_REF], exploiting the polyconvexity hypothesis. Our proposal will be to use the "second-order" homogenization method for hyperelastic composites developed by Lopez-Pamies and Ponte Castañeda [START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF]. This is a general homogenization technique, exact to second order in the heterogeneity contrast, which has the capability to incorporate statistical information beyond the phase volume fraction and that can be applied to large classes of hyperelastic composites including reinforced and porous rubbers. For completeness, a brief outline of the second-order method is included in the following section. For a more detailed description of the theory, the reader is referred to [START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF].

Outline of the Second-Order Variational Methods

The key concept behind the second-order homogenization method for hyperelastic composites developed by Lopez-Pamies and Ponte Castañeda [START_REF] Lopez-Pamies | Second-order estimates for the large-deformation response of particle-reinforced rubbers[END_REF][START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF], as well as the earlier version of the method [START_REF] Castañeda | A second-order homogenization procedure in finite elasticity and applications to black-filled elastomers[END_REF] is the introduction of a fictitious "linear comparison composite" (LCC) with the same microstructure as the nonlinear composite (i.e., the same χ (r) ). Thus, the stored-energy function of the LCC can be formally expressed as

W T (X, F) = N r=1 χ (r) (X) W (r) T (F), (6)
where the stored-energy functions of the phases W (r) T are given by the second-order Taylor approximations of the nonlinear stored-energy functions W (r) about certain uniform reference deformation gradients F (r) :

W (r) T (F) = W (r) (F (r) ) + S (r) (F (r) ) • (F -F (r) ) + 1 2 (F -F (r) ) • L (r) 0 (F -F (r) ). ( 7
)
In this relation, the L (r) 0 are fourth-order constant tensors, which together with the F (r) , are left to be specified later, and use has been made of the notation:

S (r) (F) . = ∂W (r) (F) ∂F . ( 8
)
Note further that the constitutive relation of the phases in the LCC are given by the expressions:

S = L (r) 0 F + S (r) , (9)
where S (r) = S (r) (F (r) ) -L (r) 0 F (r) is a fixed polarization stress for fixed L (r) 0 and F (r) . Thus, the LCC can be thought of as a linear "thermoelastic" composite, but in a generalized sense since the relevant "stress" and "strain" measures are not the standard linearized measures of stress and strain [START_REF] Castañeda | A second-order homogenization procedure in finite elasticity and applications to black-filled elastomers[END_REF]. While it is well known that such a material model is not suitable to describe the constitutive behavior of actual elastomers at finite strain, the LCC is only an intermediate construction that will allow the simplification of the original fully nonlinear problem, as described by ( 4) with ( 1) and ( 2). Next, "corrector" functions measuring the error made in the approximation of the stored-energy functions W (r) of the nonlinear composite by the corresponding stored-energy functions W (r) T of the LCC are introduced such that

V (r) (F (r) , L (r) 0 ) = stat F (r) W (r) ( F (r) ) -W (r) T ( F (r) ) , ( 10
)
where the optimization operation stat with respect to a variable means differentiation with respect to that variable and setting the result equal to zero to generate an expression for the optimal value of the variable. Then, it follows that the local stored-energy functions of the phases of the nonlinear composite may be approximated as r) , L (r) 0 ), [START_REF] Hashin | Large isotropic elastic deformation of composites and porous media[END_REF] and therefore that the effective stored-energy function W of the nonlinear composite may be correspondingly approximated as

W (r) (F) ≈ W (r) T (F) + V (r) (F (
W (F) ≈ W T (F; F (s) , L (s) 0 ) + N r=1 c (r) V (r) (F (r) , L (r) 0 ), (12)
where

W T (F; F (s) , L (s) 0 ) = min F∈K W T (X, F) = min F∈K N r=1 c (r) W (r) T (F) (r) ( 13
)
is the effective energy function associated with the LCC. In this last expression, it is implicitly assumed that the variables L (s) 0 remain strongly elliptic at least up to the point where macroscopic instabilities could develop. These macroscopic instabilities correspond to the loss of strong ellipticity of the effective constitutive relation for the nonlinear composite; the use of the approximation (12) could become questionable beyond the onset of these instabilities.

The approximation [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF] is valid for any reference deformations F (r) and moduli L (r) 0 , which naturally suggests its optimization with respect to these variables. In fact, the solution of this optimization problem for the variables F (r) and L (r) 0 in the estimate [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF] for W depends on the solution of the optimization problems [START_REF] Gibson | Cellular Solids[END_REF] defining the "corrector" functions V (r) . Thus, stationarity with respect to the variables F (r) in [START_REF] Gibson | Cellular Solids[END_REF] leads to the conditions: S (r) ( F (r) ) -S (r) (F (r) ) = L (r) 0 ( F (r) -F (r) ), [START_REF] Hill | Elastic potentials and the structure of inelastic constitutive laws[END_REF] which correspond to the linearization of the nonlinear constitutive relation for the hyperelastic material in the phases. Now, it is observed that relations [START_REF] Hill | Elastic potentials and the structure of inelastic constitutive laws[END_REF] have several possible solutions. For example, these conditions could be identically satisfied by taking the F (r) equal to the F (r) , which would lead to the so-called "tangent" approximation. On the other hand, solutions are also possible such that F (r) = F (r) and F (r) = 0, which leads to a new type of linearization which has been referred to as a "generalized secant" approximation [START_REF] Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations[END_REF]. Depending on the type of linearization chosen, the general expression [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF] can deliver, as discussed in the next sections, various types of second-order estimates. It is recalled that the reason for calling these estimates "second-order" is related to the fact that they are exact to second-order in the limit of small heterogeneity contrast among the phases [START_REF] Suquet | Small-contrast perturbation expansions for the effective properties of nonlinear composites[END_REF].

TANGENT SECOND-ORDER ESTIMATES (VERSION 1)

Considering the limit as F (r) tends to F (r) in ( 14) makes the functions V (r) vanish identically. Then, optimality of the reference deformations F (r) in ( 12) leads to the prescriptions:

F (r) = F (r) . = F (r) , ( 15
)
where F (r) has been defined as the average deformation field over phase r in the LCC, defined by relations ( 6) and [START_REF] Gent | A new constitutive relation for rubber[END_REF]. As shown by Ponte Castañeda and Tiberio [START_REF] Castañeda | A second-order homogenization procedure in finite elasticity and applications to black-filled elastomers[END_REF], under condition [START_REF] Knowles | On the ellipticity of the equations of nonlinear elastostatics for a special material[END_REF], the general second-order estimate (12) simplifies to

W (F) = N r=1 c (r) W (r) (F (r) ) + 1 2 S (r) (F (r) ) • (F -F (r) ) . ( 16
)
A key disadvantage of the estimates ( 16) is that by setting F (r) = F (r) , the optimality conditions for the moduli L (r) 0 in expression [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF], which for this case specialize to

(F -F (r) ) ⊗ (F -F (r) ) (r) = 0, ( 17
)
where it is recalled that F = F(X) is the solution of the LCC problem (13), cannot be satisfied in general [START_REF] Castañeda | Variational second-order estimates for nonlinear composites[END_REF]. Instead, we implement the physically motivated prescription:

L (r) 0 = L (r) t . = ∂ 2 W (r) ∂F 2 (F (r) ), (18)
which, as already mentioned, is consistent with the limit F (r) → F (r) in [START_REF] Hill | Elastic potentials and the structure of inelastic constitutive laws[END_REF]. Note that the estimate ( 16) then depends exclusively on the average fields F (r) over the phases of the linear comparison composite.

SECOND-ORDER ESTIMATES WITH FLUCTUATIONS (VERSION 2)

Considering the more elaborate "generalized secant" linearization scheme, where F (r) = F (r) and F (r) = 0, and optimizing the resulting expression [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF] with respect to the moduli L (r) 0 , leads to the following conditions:

(F -F (r) ) ⊗ (F -F (r) ) (r) = ( F (r) -F (r) ) ⊗ ( F (r) -F (r) ). ( 19
)
Unlike, the corresponding conditions (17) associated with the tangent secondorder estimate (Version 1), the new conditions [START_REF] Lopez-Pamies | Second-order estimates for the large-deformation response of particle-reinforced rubbers[END_REF] are more flexible in that they allow fluctuations of the fields in the phases of the linear comparison composite. The question then becomes what is the optimal choice of the reference deformations F (r) . Unfortunately, this is a difficult question that has not yet been completely resolved [START_REF] Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations[END_REF]. For this reason, two different approximate choices -neither of which is expected to be optimal -will be considered here. The first is to identify the reference deformations with the phase averages of the deformation gradients in the LCC:

F (r) = F (r) . ( 20
)
This prescription was used by Lopez-Pamies and Ponte Castañeda [START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF] and has the advantage that it makes stationary (with respect to F (r) ) the stored energy W T of the LCC. Furthermore, this prescription can be shown [START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF] to lead to the following second-order estimate:

W (F) = N r=1 c (r) W (r) ( F (r) ) -S (r) (F (r) ) • ( F (r) -F (r) ) , (21)
where the phase moduli tensors L (r) 0 in the LCC are determined by the conditions:

S (r) ( F (r) ) -S (r) (F (r) ) = L (r) 0 ( F (r) -F (r) ), ( 22
)
and the variables F (r) by the conditions:

( F (r) -F (r) ) ⊗ ( F (r) -F (r) ) = (F -F (r) ) ⊗ (F -F (r) ) (r) . = C (r) F . ( 23
)
In this last relation, it is useful to note that the phase fluctuation covariance tensors C (r) F may be estimated via

C (r) F = 2 c (r) ∂ W T ∂L (r) 0 F (r) =F (r) . ( 24
)
It should be emphasized that, because the phase fluctuation tensors C (r) F are not of rank 1, it is not possible to satisfy conditions [START_REF] Lopez-Pamies | Second-order estimates for the large-deformation response of particle-reinforced rubbers[END_REF] in full generality. Instead, as explained later in more detail, only appropriate traces of these expressions should be enforced [START_REF] Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations[END_REF].

SECOND-ORDER ESTIMATES WITH FLUCTUATIONS: A SIMPLIFIED VERSION (VERSION 3)

An alternative choice for the reference deformation F (r) , which is also probably not optimal, is provided by

F (r) = F. ( 25
)
This prescription has the advantage of being simpler than [START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF], while still keeping dependence on the field fluctuations.

Using condition [START_REF] Castañeda | The overall constitutive behaviour of nonlinearly elastic composites[END_REF], together with the appropriate specialization of ( 19), the general second-order estimate [START_REF] Hill | On constitutive macro-variables for heterogeneous solids at finite strain[END_REF] can be shown to reduce to

W (F) = N r=1 c (r) W (r) ( F (r) ) -S (r) (F) • ( F (r) -F (r) ) , ( 26
)
where now the phase moduli tensors L (r) 0 in the LCC are determined by the conditions:

S (r) ( F (r) ) -S (r) (F) = L (r) 0 ( F (r) -F), (27)
and the variables F (r) by (appropriate traces of) the conditions:

( F (r) -F) ⊗ ( F (r) -F) = C (r) F + (F (r) -F) ⊗ (F (r) -F). (28)
It is seen that the second-order estimate [START_REF] Castañeda | Exact second-order estimates for the effective mechanical properties of nonlinear composite materials[END_REF] depends explicitly on the variables F (r) . In addition, as opposed to [START_REF] Knowles | On the failure of ellipticity of the equations for finite elastostatic plane strain[END_REF] and like [START_REF] Marcellini | Periodic solutions and homogenization of nonlinear variational problems[END_REF], the estimate (26) also depends directly on the variables F (r) , which are associated with the fluctuations of the deformation fields in the phases, as specified by relation [START_REF] Castañeda | Nonlinear composites[END_REF]. All three estimates, ( 16), [START_REF] Marcellini | Periodic solutions and homogenization of nonlinear variational problems[END_REF], and [START_REF] Castañeda | Exact second-order estimates for the effective mechanical properties of nonlinear composite materials[END_REF], can be shown to be exact to second order in the heterogeneity contrast, provided that the corresponding estimates for the LCC are also taken to be exact to second order in the contrast. For instance, the fact that both the Hashin-Shtrikman (HS) and Self-Consistent (SC) estimates are exact to second order in the contrast for linear composites implies that the corresponding hyperelastic HS and SC estimates for W obtained from ( 16), [START_REF] Marcellini | Periodic solutions and homogenization of nonlinear variational problems[END_REF], and (26) will be also exact to second order in the contrast. However, it should also be recalled [START_REF] Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations[END_REF] that the second-order methods exhibit a small "duality gap", which has the implication that the overall stress-strain relation for the nonlinear hyperelastic composite, as generated from equation (3), is not exactly the same as that for the LCC.

As a final remark it is noted that Lopez-Pamies and Ponte Castañeda [START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF] have shown that the estimates delivered by [START_REF] Marcellini | Periodic solutions and homogenization of nonlinear variational problems[END_REF] were superior to those delivered by [START_REF] Knowles | On the failure of ellipticity of the equations for finite elastostatic plane strain[END_REF], which (as previously stated) do not take into account field fluctuations, in the context of incompressible elastomers reinforced with rigid particles. In this particular context, the incorporation of fluctuations by the version [START_REF] Marcellini | Periodic solutions and homogenization of nonlinear variational problems[END_REF] proved to be crucial in order to recover the correct overall incompressibility constraint. In the present context of porous elastomers, there is not such an incompressibility constraint as the overall behavior of a porous elastomer remains compressible even when the matrix phase is incompressible. Instead, the challenge for porous systems is to predict correctly the evolution of the relevant microstructural variable (the porosity).

Effective Behavior of Porous Elastomers

In this section, the second-order estimates ( 16) and ( 26) for the effective storedenergy function W are applied to the special case of two-phase composites consisting of vacuous (i.e., W (2) = 0) inclusions with given initial volume fraction c (2) = f o in a compressible elastomeric matrix with stored-energy function W (1) = W . Note that for this particular case the average stress in phase 2 is identically zero, so that the average stress in phase 1 is given by S (1) = (1/c (1) )S. It is also emphasized that by exploiting the objectivity of W , only macroscopic pure stretch loading histories (i.e., F = U; R = I) need to be considered. Here, general expressions are derived for two types of geometry and distribution of the pores: (i) initially spherical pores distributed randomly and isotropically in the reference configuration, and (ii) aligned cylindrical pores with initially circular cross section distributed randomly and isotropically in the transverse plane of the reference configuration. It is remarked that the former situation corresponds to a statistically isotropic composite, whereas the latter corresponds to a statistically transversely isotropic one. The loading will be assumed to be general three-dimensional in the first case, and in-plane two-dimensional in the second. Before proceeding with the various second-order estimates, it is useful to recall the classical Voigt upper bound, which depends only on the phase concentrations and follows from the principle of minimum potential energy [START_REF] Ogden | Extremun principles in non-linear elasticity and their application to composites -I Theory[END_REF]. Thus, the specialization of this Voigt bound to porous elastomers with hyperelastic matrix phase W leads to

W V (U) = (1 -f o )W (U). ( 29
)
Note that in the limit when the matrix phase is made incompressible, the Voigt upper bound becomes infinite for all loadings except for those with macroscopically isochoric deformations (i.e., det F = J = 1). Although rigorously an upper bound, the Voigt estimate is physically unrealistic, because it would suggest that a porous elastomer with an incompressible matrix phase would be itself incompressible, which is in contradiction with experimental evidence. This spectacular failure of the Voigt bound can be used as motivation for generating the new types of estimates that we propose to develop in this paper. Although they are less rigorous in the sense that they will not be bounds, they will be much more accurate and will give realistic predictions, at least for non-isochoric overall deformations. We conclude this section by noting that the available lower bounds [START_REF] Castañeda | The overall constitutive behaviour of nonlinearly elastic composites[END_REF] vanish identically for the case of porous elastomers, and are therefore also of little practical value in this context.

THE LINEAR COMPARISON COMPOSITE

In order to compute the second-order estimates ( 16) and ( 26) for porous rubbers it is necessary to determine the effective stored-energy function ( 13) associated with a fictitious linear porous "thermoelastic" composite (LCC) with the same microstructure as the original elastomer, as well as the corresponding phase averages

F (r)
and fluctuations C (r) F . Note that this fictitious linear thermoelastic problem involves measures of stress and deformation that are not symmetric, and hence suitable generalizations of the classical thermoelastic problem are required. These generalizations are straightforward and were provided by Ponte Castañeda and Tiberio [START_REF] Castañeda | A second-order homogenization procedure in finite elasticity and applications to black-filled elastomers[END_REF] in the broader context of N-phase thermoelastic composites. For conciseness, the general expressions will not be repeated here, instead, only the relevant results specialized to two-phase porous systems will be considered. To this end, it is first noted that for the special class of two-phase composites, the work of Levin [START_REF] Levin | Thermal expansion coefficients of heterogeneous materials[END_REF] allows great simplification of the general relations of linear thermoelastic composites. In fact, the effective potential energy (13) for W T for two-phase composites may be simply written in the form:

W T (F) = f + T • H + 1 2 H • L 0 H + 1 2 H + ( L 0 ) -1 ( T) • ( L 0 -L 0 ) H + ( L 0 ) -1 ( T) , ( 30
)
where the notations H = F -I and H (r) = F (r) -I have been introduced for convenience. Also here,

f (r) = W (r) (F (r) )-T (r) •H (r) -1 2 H (r) •L (r) 0 H (r) with T (r) = S (r) (F (r) ) -L (r)
0 H (r) , and L 0 = L (1) 0 -L (2) 0 , T 0 = T (1) 0 -T (2) 0 . Furthermore, in this relation, f and L 0 are the volume averages of f and L 0 , while L 0 is the tensor of effective modulus of the two-phase, linear-elastic comparison composite with moduli L (1) 0 and L (2) 0 , and the same microstructure, in its undeformed configuration, as the nonlinear hyperelastic composite. Now, by letting L (2) 0 → 0 and defining L 0 = L (1) 0 and M 0 = L -1 0 , it is straightforward to show that relation [START_REF] Castañeda | Variational second-order estimates for nonlinear composites[END_REF] in the case of porous systems specializes to (1) + M 0 S(F (1) ) • L 0 F -F (1) + M 0 S(F (1) ) . [START_REF] Sanchez-Palencia | Comportements local et macroscopique d'un type de milieux physiques heterogenes[END_REF] Next, the average deformation F (1) and the fluctuations C (1) F in the matrix phase can be determined from the stored-energy function (31) (see, for example, [START_REF] Castañeda | Nonlinear composites[END_REF]), and may be simply written as (1) + M 0 S(F (1) ) [START_REF] Suquet | Small-contrast perturbation expansions for the effective properties of nonlinear composites[END_REF] and

W T (F) = (1 -f o )W (F (1) ) - 1 -f o 2 S(F (1) ) • M 0 S(F (1) ) + 1 2 F -F
F (1) = F + 1 1 -f o M 0 L 0 -(1 -f o )L 0 F -F
C (1) F = 2 1 -f o ∂ W T ∂L 0 , ( 33
)
respectively. Note that expressions (31) through [START_REF] Triantafyllidis | On the comparison between microscopic and macroscopic instability mechanisms in a class of fiber-reinforced composites[END_REF] are valid for any reference deformation tensor F (1) and modulus L 0 . Furthermore, these expressions are valid for any effective modulus tensor L 0 . For example, for the case of a random and isotropic distribution of initially spherical pores, use can be made of the isotropic Hashin-Shtrikman-type estimate [START_REF] Willis | Variational and related methods for the overall properties of composites[END_REF]:

L 0 = L -1 0 + f o 1 -f o Q -1 -1 , ( 34
)
where Q = L 0 -L 0 PL 0 , with P being obtained by setting L (0) equal to L 0 in the expression:

P = 1 4π |ξ |=1 H (0) (ξ ) dS, (35)
with

K (0) ik = L (0) ij kl ξ j ξ l , N (0) = K (0) -1
, and H (0) ij kl (ξ ) = N (0) ik ξ j ξ l . Similarly, for the case of a random and isotropic distribution of aligned cylindrical pores with initially circular cross section, the corresponding HS type estimate is given by [START_REF] Wang | Stability and vibrations of elastic thick-walled cylindrical and spherical shells subjected to pressure[END_REF] but with P being obtained by setting L (0) equal to L 0 in the expression [START_REF] Willis | Variational and related methods for the overall properties of composites[END_REF]:

P = 1 2π ξ 2 1 +ξ 2 2 =1 H (0) (ξ 1 , ξ 2 , ξ 3 = 0) dS, (36)
where the cylindrical pores have been aligned in the x 3 direction. From a computational point of view, it is seen that the microstructural tensor P depends on the anisotropy of the modulus L 0 , which in turn depends on the functional form of the potential W , as well as the particular type of loading, as determined by F = U.

TANGENT SECOND-ORDER ESTIMATES (VERSION 1)

The specialization of the second-order estimate [START_REF] Knowles | On the failure of ellipticity of the equations for finite elastostatic plane strain[END_REF] to the case of elastomeric porous composites leads to

W (F) = (1 -f o ) W (F (1) ) + 1 2 S(F (1) ) • (F -F (1) ) , ( 37
)
where the variable F (1) needs to be determined. Now, by recognizing that within the framework of the tangent second-order estimates F (1) = F (1) , equation ( 32) can be readily shown to reduce to 1) + M 0 S(F (1) ) .

F (1) = F + 1 1 -f o M 0 L 0 -(1 -f o )L 0 F -F (
(38)

Next, recalling that in this context L 0 = ∂ 2 W (F (1) )/∂F 2 , it is seen that (38) constitutes a system of nine nonlinear algebraic equations for the components of the average deformation F (1) . The solution of these equations can then be used to compute the effective stored-energy function W for the porous elastomer using (37).

SECOND-ORDER ESTIMATES WITH FLUCTUATIONS (VERSION 3)

The second-order estimate [START_REF] Castañeda | Exact second-order estimates for the effective mechanical properties of nonlinear composite materials[END_REF] for the case of elastomeric porous composites specializes to (1) ) -S(F) • ( F (1) -F (1) ) .

W (F) = (1 -f o ) W ( F
(39)

Here, the variables F (1) , F (1) , as well as the modulus tensor L 0 of the matrix phase of the linear comparison composite, need to be determined. Now, recognizing that in connection with the second-order estimate (26) the reference field F (1) = F, equation ( 32) can be shown to yield the following expression for the average deformation gradient F (1) associated with the estimate (39):

F (1) = F + 1 1 -f o M 0 L 0 -(1 -f o )L 0 M 0 S(F). ( 40
)
With regard to the above equations, it is necessary to remark that (40) provides an explicit expression for F (1) in terms of the modulus L 0 , which can be determined, along with F (1) , making use of the generalized secant condition:

S( F (1) ) -S(F) = L 0 ( F (1) -F), (41)
together with suitably chosen traces of the expression:

( F (1) -F) ⊗ ( F (1) -F) = 2 (1 -f o ) ∂ W T ∂L 0 . ( 42
)
More specifically, the traces to be taken depend on the choice of the form of L 0 , as will be explained later in more detail. Moreover, in this last expression, W T is the stored-energy function for the relevant LCC given by

W T (F) = (1 -f o )W (F) + 1 2 S(F) • M 0 L 0 -(1 -f o )L 0 M 0 S(F). ( 43
)

THE STRONG ELLIPTICITY CONDITION

A complete study of the stability of porous elastomers with random microstructures is an extremely difficult problem, and well beyond the scope of this work. However, following Geymonat et al. [START_REF] Geymonat | Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity[END_REF], the onset of macroscopic instabilities can be identified from the loss of strong ellipticity of the effective constitutive behavior of the porous elastomers, which, as has already been seen in the prior subsections, can be estimated easily and efficiently by means of the second-order variational procedure.

In this subsection, it is quickly recalled that the condition of strong ellipticity for a given constitutive relation is that the corresponding acoustic tensor must be positive definite. More precisely, the condition of strong ellipticity for the homogenized porous elastomer characterized by the effective stored-energy function W can be written in the form:

L ij kl N j N l m i m k > 0 (44)
for all m ⊗ N = 0. Here, L = ∂ 2 W /∂F 2 is the fourth-order tensor of first-order elastic moduli characterizing the effective incremental behavior of the porous composite. When condition (44) is satisfied, the associated macroscopic equilibrium equations form a strongly elliptic system of partial differential equations.

In connection with condition (44) it should be emphasized that L = L 0 . In other words, the modulus associated with the effective stored-energy function of the porous elastomer does not correspond exactly to the effective modulus associated with the auxiliary linear thermoelastic composite. This is related to the abovementioned existence of a (small) "duality gap" in the second-order variational estimates.

Parenthetically, it is recalled that the condition of ordinary ellipticity requires the acoustic tensor to be merely nonsingular and not necessarily positive definite. Hence, strong ellipticity implies ordinary ellipticity but the converse is not true in general. However, the interest here is in the determination of the loss of strong ellipticity of homogenized materials that are strictly convex, and therefore strongly elliptic, under infinitesimal deformations. Then, elliptic regions which contain the point λ1 = λ2 = λ3 = 1, by continuity, will be necessarily strongly elliptic as well. That is, for the cases studied here, ellipticity and strong ellipticity are equivalent.

In summary, use can be made of the second-order expressions (37) and (39) to produce estimates for the effective stored-energy function W in order to determine the strongly elliptic domain of porous elastomers through condition (44).

Plane Strain Loading of Transversely Isotropic Porous Elastomers

The results presented in the previous section are specific for porous elastomers with (3-D and 2-D) isotropic microstructures, hence the use of the Hashin-Shtrikman (HS)-type estimate [START_REF] Wang | Stability and vibrations of elastic thick-walled cylindrical and spherical shells subjected to pressure[END_REF] for the LCC. However, the results are general as far as the material behavior and loading conditions are concerned. The aim of this paper is to make use of these results for the first time, and therefore some additional hypotheses will be made in this section in order to simplify the calculations involved, in an attempt to generate results that are as explicit as possible. Thus, the problem of plane-strain loading of porous elastomers consisting of cylindrical voids (with initially circular cross-section), perpendicular to the plane of deformation and aligned in the x 3 direction, will be considered. Note that the applied deformation F = U in this context is entirely characterized by the macroscopic principal stretches λ1 and λ2 , with the out-of-plane principal stretch λ3 being identically 1. Furthermore, attention will be restricted to porous elastomers with isotropic matrix phases. This restriction, along with objectivity and plane strain conditions, implies that the stored-energy function of the matrix phase can be expressed as a function of the principal invariants of the right Cauchy-Green deformation tensor C = F T F:

I = tr C = λ 2 1 + λ 2 2 , (45) J = √ detC = λ 1 λ 2 ,
or equivalently, as a symmetric function of the principal stretches λ 1 , λ 2 of F. Therefore, W may be written as

W (F) = ϕ(I, J ) = (λ 1 , λ 2 ). ( 46
)
However, actual rubber being nearly incompressible, it will suffice to consider isotropic matrix phases characterized by the stored-energy function:

W (F) = g(I ) + h(J ) + μ 2 (J -1) 2 , ( 47
)
where the Lamé parameter μ will be taken to tend to infinity in order to recover incompressible behavior (J → 1). Here, g and h are assumed to be twice continuously differentiable and to satisfy the conditions: g(2) = h(1) = 0, g I (2) = μ/2, h J (1) = -μ, and 4g I I (2) + h J J (1) = μ, where the subscripts I and J denote partial differentiation with respect to these invariants. Note that when these conditions are satisfied W (F) → (1/2)μ (trε) 2 + μtrε 2 , where ε is the infinitesimal strain tensor, as F → I, so that the stored-energy function (47) linearizes properly.

A relatively simple model of the general type (47), which captures the limiting chain extensibility of elastomers, is the Gent model [START_REF] Gent | A new constitutive relation for rubber[END_REF]:

W (F) = - μJ m 2 ln 1 - I -2 J m -μ ln J + μ 2 - μ J m (J -1) 2 , ( 48
)
where the parameter J m is the limiting value for I -2 at which the material locks up. Note that (48) reduces to a Neo-Hookean material on taking the limit J m → ∞.

5.1. TANGENT SECOND-ORDER ESTIMATES (VERSION 1)

Compressible Matrix

Within the framework of the second-order estimate (37), it is seen that under plane strain loading and given matrix material behavior (47), expression (38) reduces simply to two nonlinear equations for the components F

(1)

11 = λ(1)
1 and F

(1)

22 = λ(1)
2 (the other in-plane components being identically zero) of the average deformation F (1) in the matrix phase of the linear comparison composite. They are:

f o ( λ(
11 = 0, (49)
where

L 1111 = 2 ḡ(1) I + 4 λ(1) 1 2 ḡ(1) I I + h(1) J J + μ λ(1) 2 2 , L 2222 = 2 ḡ(1) I + 4 λ(1) 2 2 ḡ(1) I I + h(1) J J + μ λ(1) 1 2 , L 1122 = h(1) J -μ + 4 ḡ(1) I I + h(1) J J + 2μ λ(1) 1 λ (1) 2 and

S

(1)

11 = 2 ḡ(1) I λ(1) 1 + h(1) J λ(1) 2 + μ λ(1) 1 λ(1) 2 -1 λ(1) 2 , S (1) 22 = 2 ḡ(1) I λ(1) 2 + h(1) J λ(1) 1 + μ λ(1) 1 λ(1) 2 -1 λ(1) 1 .
Here, for convenience, the subscript '0' has been dropped out for the components of L 0 . Also, use has been made of the notation ḡ(1) = g( Ī (1) ) and h(1) = h(J (1) ),

where

Ī (1) = ( λ(1) 1 ) 2 + ( λ(1)
2 ) 2 and J (1) = λ(1)

1 λ(1)
2 denote the principal invariants associated with F (1) . In passing, it is noted that for the tangent modulus tensor of (47) and the applied loading, the relevant tensor P can be computed analytically, but for brevity, the explicit expressions will not be included here. In summary, upon solving numerically the two nonlinear algebraic equations (49) for the nonzero components of F (1) , the effective stored-energy function W can be readily computed using (37).

Incompressible Matrix

The limit when the matrix phase is made incompressible (i.e., μ → ∞) can be shown to simplify the above expressions considerably. As already stated, this limit is interesting from a practical perspective, given that actual elastomers exhibit a nearly incompressible behavior (i.e., they usually exhibit a ratio between the bulk and shear moduli of the order of 10 4 ). Due to the cumbersomeness of the final expressions, the asymptotic analysis of the incompressible limit involving the general stored-energy function (47) will not be included here. Instead, for illustrative purposes, only the particular case of a Neo-Hookean matrix phase will be discussed.

The details of the relevant asymptotic analysis are given in Appendix A, but the final expression for the tangent second-order estimate (37) for the effective storedenergy function of a porous elastomer with incompressible Neo-Hookean matrix phase may be written as

W I (F) = I ( λ1 , λ2 ) = (1 -f o ) 2u μ (u 2 -1)( λ1 -λ2 ) + (1 + u 2 )( λ2 u 2 -2u + λ1 )((1 + f o )u 2 + ( λ2 -λ1 )u -1 -f o ) u( λ1 -λ2 u 2 ) . ( 50
)
In this expression, u is the solution of the equation

( λ2 2 + f 2 o -1)u 4 + 2( λ1 + (f o -1) λ2 )u 3 + ( λ2 2 -λ2 1 )u 2 -2((f o -1) λ1 + λ2 )u + 1 -f 2 o -λ2 1 = 0, (51)
which can be determined explicitly as a function of f o , λ1 , and λ2 . Note that (51) is a quartic polynomial equation and hence it may be solved in closed-form. However, for practical purposes, it proves helpful to solve (51) numerically. In this regard, it is emphasized that only one of the 4 roots of (51) gives the correct linearized behavior for the effective response of porous materials; this is indeed the root that should be chosen at least initially.

For the special case of in-plane hydrostatic finite expansion or contraction, i.e., λ1 = λ2 = λ, the second-order estimate (50) can be shown to further reduce to

I ( λ, λ) = 2μ (1 -f o )( λ -1) 2 f o + λ -1 . ( 52
)
For later use, it is also noted that the average deformation field in the matrix phase associated with the stored-energy function ( 52) is given by

F (1) = I. (53)
Finally, it is noted that the result (53) holds true, not only for a porous elastomer with incompressible Neo-Hookean matrix, but for a porous elastomer with a general incompressible isotropic matrix phase.

5.2. SECOND-ORDER ESTIMATES WITH FLUCTUATIONS (VERSION 3)

Compressible Matrix

Given the transverse isotropy of the microstructure and the orthogonal symmetry of the prescribed loading, it is reasonable to assume that the linear comparison problem of relevance here will also exhibit orthotropic symmetry, with the symmetry axes aligned with the applied loading F = U. In order to carry out the computation for the HS-type second-order estimates (39) for porous elastomers with matrix phase (47) under plane strain conditions, it suffices to consider the inplane components of a general deformation field F relative to the symmetry axes.

It is convenient to denote these components as a vector in 4 [

F 11 F 22 F 12 F 21 ] T . ( 54
)
The modulus tensor L 0 of the matrix phase of the linear comparison composite is consequently denoted as a matrix in 4×4 ⎡ ⎢ ⎢ ⎣

L 1111 L 1122 0 0 L 1122 L 2222 0 0 0 0 L 1212 L 1221 0 0 L 1221 L 2121 ⎤ ⎥ ⎥ ⎦ , ( 55
)
where, for convenience, the subscript '0' has been suppressed for the components of L 0 , and use has been made of the major symmetry of the tensor L 0 (i.e., L ij kl = L klij ). Next, due to the above-stated assumptions, it is seen that the tensor F (1) has at most 4 independent components ( F (1) 11 , F (1) 22 , F (1) 12 , F (1) 21 ) which should be determined from equations (42). This implies that the modulus L 0 must be chosen to have at most four independent components, with respect to which W T will be differentiated to generate 4 consistent equations for the components of F (1) . Note that these four conditions correspond to 4 different "traces" of the equations (42). For this reason, the components of L 0 will be required to satisfy the prescriptions:

L 1212 = L 2121 and L 1221 + L 1122 = (L 1111 -L 1212 )(L 2222 -L 1212 ). ( 56
)
These relations reduce the components of the L 0 to only 4 independent ones, namely, L 1111 , L 2222 , L 1122 , and L 1212 . It is recalled [START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF] that the motivation for the choices (56) was twofold: (i) the components of the tangent modulus of a Neo-Hookean material, expressed relative to the symmetry axes, satisfy (56); and (ii) the conditions (56) simplify significantly the computation of the microstructural tensor P (see [20, Appendix A]). Now, making use of the equations [START_REF] Wang | Stability and vibrations of elastic thick-walled cylindrical and spherical shells subjected to pressure[END_REF] for the HS estimate for L 0 , the equations (42), together with (43), can be used to generate 4 equations for the 4 components of F (1) , which are of the form:

F (1)
11 -λ1 2 + 2f 1 F (1) 12 F (1) 21 = k 1 , F (1) 22 -λ2 2 + 2f 2 F (1) 12 F (1) 21 = k 2 , (57) F (1) 12 2 + F (1) 21 2 + 2f 3 F (1) 12 F (1) 21 = k 3 , F (1) 11 -λ1 F (1) 22 -λ2 -F (1) 12 F (1)

21 = k 4 ,
where

f 1 , f 2 , f 3 , k 1 , k 2 , k 3 , k 4 are
functions of the components of L 0 , i.e., L 1111 , L 2222 , L 1122 , and L 1212 , as well as of the deformation F, the initial volume fraction of voids f o , and the constitutive functions g, h, and μ of the matrix phase. Equations (57) can be shown to yield two distinct solutions for F (1) 11 and F (1) 22 , in terms of which F (1) 12 and F (1) 21 may be determined. Note, however, that the variables F (1) 12 and F (1) 21 only enter the equations through the combinations F (1) 12 F (1) 21 and ( F (1) 12 ) 2 + ( F (1) 21 ) 2 , and hence, only these combinations can be determined uniquely from (57). The two solutions for F (1) 11 and F (1) 22 are as follows:

F (1) 11 -λ1 = ± 2f 1 k 4 + k 1 4f 2 1 k 2 + 4f 1 k 4 + k 1 , ( 58
) F (1) 22 -λ2 = ±

2f 1 k 2 + k 4 4f 2 1 k 2 + 4f 1 k 4 + k 1
, where it must be emphasized that the positive (and negative) signs in the roots for F (1) 11 and F (1) 22 go together. Next, each of the two distinct roots of (57) can be substituted into the generalized secant condition (41) to obtain a system of 4 equations for the 4 unknowns L 1111 , L 2222 , L 1122 , and L 1212 , which must be solved numerically. Having computed the values of all the components of L 0 for a given initial porosity (f o ), given material behavior (g, h, and μ ), and given loading ( λ1 and λ2 ), the values of the components of F (1) and F (1) can be readily determined using relations (40) and (58), respectively. In turn, the second-order estimate (39) for the effective stored-energy function W for porous isotropic elastomers can be computed using these results. Finally, it is noted that the two above roots lead to very similar results for the effective behavior of the porous elastomer when small values of μ are considered (i.e., for μ of the order of μ). However, for larger values of μ , the estimates produced by the two distinct roots are very different. In fact, as explained in more detail in the following subsection, it can be shown that for large values of μ only one of the two roots provides physically meaningful results. Consequently, this is the root that should be chosen to compute the effective behavior of the porous elastomer.

Incompressible Matrix

The above expressions can be simplified in the limit of incompressibility of the matrix phase (i.e., μ → ∞). In this context, it is important to realize that the asymptotic behavior of one of the two above roots leads to nonphysical predictions in the limit as μ becomes unbounded. More specifically, for deformations satisfying ē1 + ē2 0 only the "positive" (+) root provides physically reasonable predictions, whereas for deformations with ē1 + ē2 0, only the "negative" (-) root has the physically plausible asymptotic behavior; here, the logarithmic strains ē1 = ln( λ1 ) and ē2 = ln( λ2 ) have been introduced for convenience. Taking this observation into account, it can be shown that the second-order estimate (39) for the the effective stored-energy function of a porous elastomer with an incompressible isotropic matrix phase may be written as

W I (U) = I ( λ1 , λ2 ) = (1 -f o )g( Î (1) ), ( 59
)
where Î (1) = Î (1) (α, β, γ ), and α, β, γ are the solution to three nonlinear algebraic equations, not shown here for their bulkiness, which can be solved for in terms of the initial porosity f o , given material behavior g, and given loading λ1 and λ2 . In general, it is not possible to solve these equations in closed-form. However, for the particular case of a porous elastomer with an incompressible Neo-Hookean matrix phase, the general estimate (59) can be shown (see Appendix B) to reduce to

I ( λ1 , λ2 ) = (1 -f o )μ 2 p 4 v 4 + p 3 v 3 + p 2 v 2 + p 1 v + p 0 (q 2 v 2 + q 1 v + q 0 ) 2 -2 , ( 60
)
where v is the solution of the quartic polynomial:

r 4 v 4 + r 3 v 3 + r 2 v 2 + r 1 v + r 0 = 0. (61)
Here, the coefficients p 0 , p 1 , p 2 , p 3 , p 4 , q 0 , q 1 , q 2 , r 0 , r 1 , r 2 , r 3 and r 4 , which depend on f o , λ1 and λ2 , are given in explicit form in Appendix C. Since the es-timate (60) depends effectively on the solution of the quartic polynomial equation (61), it may be written in closed-form. However, for all practical purposes, it is simpler to solve (61) numerically. In this regard, it is emphasized that only one of the 4 roots of (61) gives the correct linearized behavior for the effective response of porous materials; this is indeed the root that should be chosen at least initially. It is useful, for comparison purposes, to spell out the simplification of the second-order estimate (60) for the case of in-plane hydrostatic loading, i.e., λ1 = λ2 = λ. The result reads as follows:

I ( λ, λ) = 2μ 1 -f o (1 + f o ) λ2 + f o -1 -2 λ f o ( λ2 + f o -1) . ( 62
)
For later use, it is noted that the average deformation field in the matrix phase associated with the stored-energy function ( 62) is given by

F (1) = λI I, (63)
where

λI = f o ( λ2 + f o -1) - λ f o -1 . ( 64
)
Finally, it is emphasized that the result (64) holds true, not only for a porous elastomer with incompressible Neo-Hookean matrix, but in fact for a porous elastomer with general incompressible isotropic matrix phase.

COMPARISON OF THE SECOND-ORDER ESTIMATES WITH EXACT RESULTS

Hydrostatic Loading

With regard to porous elastomers subjected to finite deformations, there are very few exact results available. For the special case of hydrostatic loading, Hashin [START_REF] Hashin | Large isotropic elastic deformation of composites and porous media[END_REF] obtained the exact equilibrium solution by making use of the idea of the composite spheres assemblage. Following that work, it is straightforward to show that the exact stored-energy function for the in-plane hydrostatic deformation of a porous rubber with incompressible isotropic matrix W (F) = (λ 1 , λ 2 ) may be written as

W I (U) = 2 1 √ f o (λ, λ -1 )R dR, (65)
where

λ = 1 + λ2 -1 R 2 . ( 66
)
The corresponding exact average deformation in the matrix phase reads as

F (1) = λI I, (67)
with

λI = f o ( λ2 + f o -1) - λ f o -1 , ( 68
)
where λ must be greater than √ 1f o . In general, the integral in (65) cannot be computed analytically; however, for the particular case of a porous elastomer with incompressible Neo-Hookean matrix phase, the exact function may be expressed as

I ( λ, λ) = μ 2 ( λ2 -1) ln λ2 + f o -1 f o -ln( λ2 ) . ( 69
)
It can thus be seen that the two versions of the second-order estimates, as defined by ( 52) and (62), do not recover the exact result (69) for the effective storedenergy function of porous elastomers with incompressible Neo-Hookean matrix phase subjected to general in-plane hydrostatic finite deformations. Nonetheless, both estimates can be shown to be exact to third order in the infinitesimal strain (i.e., up to O( λ-1) 3 ). For larger finite deformations, however, the behavior of both estimates is very different. As it will be seen in the following section, whereas the expression (62), which takes into account the field fluctuations, provides estimates that are in very good agreement with the exact result, the corresponding tangent expression (52) delivers estimates that deviate drastically from (69). This disparity is mainly due to the difference in the prediction of the evolution of the microstructure. In fact, while the average deformation gradient in the matrix phase (53) predicted by the tangent second-order theory is exactly equal to the identity, the corresponding F (1) (63) predicted by the second-order method with fluctuations is consistent with the exact result (67). This is a remarkable result. Indeed, by taking into account the fluctuations, the second-order estimate ( 26) is able to improve on the tangent second-order estimate [START_REF] Knowles | On the failure of ellipticity of the equations for finite elastostatic plane strain[END_REF] in that it recovers the exact average deformation fields in a porous elastomer with an incompressible isotropic matrix under finite in-plane hydrostatic loading.

General Loading

For more general loadings, there are no other known exact solutions to the problem for porous elastomers. However, for incompressible matrix phase materials, a simple conservation of mass argument (for the matrix phase) allows for the determination of the evolution of the porosity f as a function of deformation. The result is

f = 1 - 1 -f o J . ( 70
)
Now, it can be shown that the estimate for the porosity associated with the secondorder estimate with fluctuations (59) for porous elastomers with incompressible isotropic matrix phases is in exact agreement with the exact result (70). The proof of this is sketched out for the particular case of a Neo-Hookean porous material at the end of Appendix B. On the other hand, the corresponding estimate for the porosity arising from the tangent second-order estimate (37) deviates radically from (70) for large deformations, as will be shown explicitly in the next section. In summary, Version 3 (with fluctuations), unlike Version 1 (tangent), of the second-order method is able to predict the exact evolution of the porosity for general finite loading, and consequently the exact average deformation fields for hydrostatic loading, for a porous elastomer with incompressible isotropic matrix phase. This is a nontrivial result, since the actual fields in a deformed porous elastomer are highly heterogeneous. However, it appears that the "generalized secant" linearization together with the incorporation of field fluctuations leads to improved approximations capable of capturing better the heterogeneity of these fields in order to deliver the exact results mentioned above.

LOSS OF STRONG ELLIPTICITY

For the particular case of plane strain deformations, the strong ellipticity condition (44) can be written more explicitly. In fact, under plane strain deformations, necessary and sufficient conditions for strong ellipticity of the effective constitutive behavior of the type of porous systems considered here have been shown [START_REF] Knowles | On the failure of ellipticity of the equations for finite elastostatic plane strain[END_REF][START_REF] Hill | On the theory of plane strain in finitely deformed compressible materials[END_REF] to reduce to:

L 1111 > 0, L 2222 > 0, L 1212 > 0, L 2121 > 0, and (71)
L 1111 L 2222 + L 1212 L 2121 -( L 1122 + L 1221 ) 2 > -2 L 1111 L 2222 L 1212 L 2121 ,
where

L iijj = ∂ 2 W ∂ λi ∂ λj , L ij ij = 1 λ2 i -λ2 j λi ∂ W ∂ λi -λj ∂ W ∂ λj , i = j, ( 72
)
L ijj i = 1 λ2 i -λ2 j λj ∂ W ∂ λi -λi ∂ W ∂ λj , i = j,
(i, j = 1, 2) are the components of the modulus L written in the Lagrangean principal axes. Note that the third and fourth conditions in (71) are equivalent and that for loadings with λi = λj (i = j ), suitable limits must be taken for some of the components in (72). In particular, equations (72) 2 and (72) 3 transform into:

L ij ij = 1 2 L iiii -L iijj + 1 λi ∂ W ∂ λi , i = j, ( 73
)
L ijj i = 1 2 L iiii -L iijj - 1 λi ∂ W ∂ λi , i = j, respectively.
It is important to remark here that most stored-energy functions of the form (47) describe best the actual behavior of elastomers when "calibrated" to be strongly elliptic. For example, a compressible Gent material, characterized by the storedenergy function (48), is strongly elliptic under plane-strain deformations if (but not only if) μ > 0, J m > 0, and μ > 2μ/J m . Note that for a Neo-Hookean elastomer, these sufficient conditions simplify to μ > 0 and μ > 0. In fact, for realistic elastomers, μ > 0, J m > 0, and μ is not only positive but is several orders of magnitude larger than μ, namely, μ /μ ≈ 10 4 . Consequently, the Gent elastomers utilized in this work to model the matrix behavior of the porous elastomers are strongly elliptic for all deformations. As a final remark, it is recalled that having strict convexity in the linearized region, i.e., μ > 0 and κ = μ + μ > 0 (where κ denotes the plane-strain bulk modulus) in the limit λ1 → 1 and λ2 → 1, does not ensure strong ellipticity of a Gent material for all deformations.

In the next section, the second-order methods described in this section will be used to generate estimates for the strongly elliptic domains of random porous elastomers with incompressible Gent and Neo-Hookean matrix phases subjected to plane-strain deformations. It will be shown that even though the behavior of the matrix phase is strongly elliptic, the homogenized behavior of the porous elastomer can lose strong ellipticity. This is consistent with the observations of Abeyaratne and Triantafyllidis [START_REF] Abeyaratne | An investigation of localization in a porous elastic material using homogenization theory[END_REF] for porous elastomers with periodic microstructures.

Results for General Plane-Strain Loading

This section presents results associated with the (Versions 1 and 3) second-order HS estimates for in-plane hydrostatic, uniaxial, and pure shear loading of porous elastomers with incompressible Gent and Neo-Hookean matrix phases. Results are given for μ = 1 and various levels of initial porosity f o , and were computed up to the point at which the effective incremental moduli were found to lose strong ellipticity, or truncated at some sufficiently large strain if no such loss was found. For clarity, the points at which loss of strong ellipticity occurred are denoted with the symbols and • for Versions 1 and 3, respectively. The characterization of the strongly elliptic domains is given in the last subsection. It is further noted that exact results and bounds are presented when available.

HYDROSTATIC LOADING

Figure 1 presents the comparison between the effective behavior as predicted by Versions 1 and 3 of the second-order method and the exact result, for a Neo-Hookean porous elastomer with incompressible matrix phase (μ → ∞) under Recall that closed-form expressions for the effective stored-energy functions shown in Figure 1(a) are given by ( 52), (62), and (69) for Version 1, Version 3, and the exact result, respectively.

The main observation that can be made from Figure 1 is that Version 3 of the second-order variational procedure provides estimates for the effective constitutive behavior which are in excellent agreement with the exact result. Version 1 also delivers estimates that compare reasonably well with the exact result for compressive loadings. However, the predictions by Version 1 deviate significantly from the exact behavior for large tensile deformations. It is also seen that both versions of the second-order method predict loss of strong ellipticity of the homogenized porous elastomer under in-plane hydrostatic compression, while no such behavior is observed under in-plane hydrostatic tension. Moreover, both second-order estimates for the effective behavior exhibit a better agreement with the exact result for higher values of f o . Finally, it is interesting to note from Figure 1(b) that the overall constitutive behavior of the composite consistently exhibits hardening under compression and softening under tension. This feature will be shown shortly to be due mainly to a geometric effect caused by the evolution of the porosity.

Figure 2 provides plots associated with the results shown in Figure 1 for: (a) the porosity as a function of the logarithmic strain ē = ln( λ); and (b) the critical stretch λcrit at which the loss of strong ellipticity of the homogenized elastomer takes place, as a function of initial porosity f o . First, a key point to be drawn from Figure 2(a) is that the porosity decreases for compressive deformations and increases for tensile ones. This entails a geometric hardening/softening on the overall response of the porous elastomer which is entirely consistent with the hardening/softening exhibited by the effective constitutive behavior observed in Figure 1(b). Moreover, the porosity predicted by Version 3 of the second-order method reduces to the exact result, as already pointed out in the previous section (see expressions (63) and ( 67)). On the other hand, the porosity delivered by Version 1 deviates from the exact evolution for large finite deformations, especially for tensile hydrostatic loading. In fact, under hydrostatic tension, Version 1, which does not take into account information about the field fluctuations, predicts unrealistic values for the porosity (i.e., greater than unity). This explains the extremely soft effective constitutive behavior observed in Figure 1 for Version 1 of the second-order method under hydrostatic tension. Also, in accordance with the trend discerned from Figure 1, it appears that the porosity evolution predicted by Version 1 gets worse with decreasing initial porosity.

The main observation with regard to Figure 2(b) is the somewhat counterintuitive result that the porous material becomes more stable ( λcrit smaller) with increasing initial values of the porosity. In this connection, it is relevant to remark that while exact results are available for the effective stored-energy function and porosity evolution in in-plane hydrostatic loading of composite cylinders (with incompressible matrix phase), the loss of strong ellipticity of these structures has not been studied. However, Wang and Ertepinar [START_REF] Wang | Stability and vibrations of elastic thick-walled cylindrical and spherical shells subjected to pressure[END_REF] did study the stability of an isolated cylindrical Neo-Hookean shell under in-plane hydrostatic loading. Results of that work comprising the buckling flexural modes n = 2, which corresponds to the collapse to an oval shape, and n = 3 have been included in Figure 2(b), for reference purposes. Although, it should be emphasized that the buckling behavior of an isolated shell cannot be rigorously identified with the buckling instabilities that would take place in an actual composite system, even for a composite with the Hashin composite-shell microstructure, it is interesting to remark that these results appear to be consistent with the results derived from the second-order theory in that the overall stability is enhanced with increasing initial porosity, at least initially. Moreover, it is noted that the critical stretches characterizing the loss of strong ellipticity predicted by both versions of the second-order variational procedure are smaller than the corresponding critical stretches associated with the buckling modes for an isolated cylindrical shell. Figure 2(b) also shows that for the interval 0 < f o < 0.4 the loss of strong ellipticity predicted by Version 3 is slightly smaller than the one obtained from Version 1. In contrast, for initial porosities higher than 0.4, the prediction of loss of strong ellipticity by Version 1 becomes smaller than the one computed from Version 3. The difference between the results of these two versions becomes more pronounced in the limit f o → 1, where λcrit → 0.73 and λcrit → 0 for Versions 1 and 3, respectively. We expect the estimate for the critical stretch associated with Version 3 to be more accurate, but we do not have an explanation for its relatively low values at high porosities. However, it should be kept in mind that it is expected that other (smaller wavelength) instabilities would take place before reaching the loss of ellipticity condition.

UNIAXIAL LOADING

In Figure 3, plots for the effective stress-strain behavior associated with Versions 1 and 3 are presented for a Neo-Hookean porous elastomer with incompressible matrix phase (μ → ∞) under uniaxial loading with λ2 = 1, λ1 = λ. The results for the stress components S 11 and S 22 are presented in parts (a) and (b), respectively, for values of f o = 30, 50 and 70%, as functions of the logarithmic strain ē = ln( λ). Similar to the case of in-plane hydrostatic loading, the results for compressive (tensile) deformations shown in Figure 3 exhibit a clear hardening (softening) behavior with increasing deformation, but less pronounced than the corresponding results for in-plane hydrostatic loading. From Figure 3(a) it is seen that the effective constitutive behavior for S 11 obtained from Version 1 is significantly softer than the one obtained from Version 3. This is also the case for the component S 22 as shown by Figure 3(b). In fact, the prediction for this component of the stress by Version 1 is not only much softer than the corresponding stress predicted by Version 3, but it even decreases for tensile loadings reaching negative values, which is physically unrealistic. This suggests that the predictions of Version 1 could be too soft for large finite deformations, especially for tensile loadings. Furthermore, as it was the case for hydrostatic loadings, loss of ellipticity was found for compressive loadings but not for tensile ones.

Figure 4 provides corresponding results for: (a) the porosity; and (b) the average aspect ratio of the pores ω, as function of the logarithmic strain ē = ln( λ). Note that the aspect ratio has been defined as ω = λ(2)

1 / λ(2) 2 , with λ(2) 1 and λ(2) 2 denoting the principal stretches associated with the average deformation gradient tensor of the vacuous phase F (2) , so that ω > (<) 1 correspond to an oblate (prolate) average shape of the pores. As it was the case for hydrostatic loading, it is seen from Figure 4(a) that the porosity decreases for compressive deformations and increases for tensile ones. In turn, this can be related to the aforementioned hardening/softening trend exhibited by the effective stress-strain behavior in Figure 3.

As already anticipated in Section 5.3.2, Figure 4(a) also shows that the prediction for the evolution of the porosity by Version 3 of the second-order method is in agreement with the exact result, whereas the prediction by Version 1 deviates from the correct behavior for large deformations. This deviation, which is much more drastic for tensile loadings and lower values of f o , helps explain the unphysical behavior observed in Figure 3(b) for S [START_REF] Müller | Homogenization of nonconvex integral functionals and cellular elastic materials[END_REF] . In particular, it is seen that S 22 tends to negative values whenever f approaches one. Figure 4(b) shows that both versions of the second-order method give similar predictions for the average aspect ratio of the pores. Note that in compression the changes in aspect ratio are more rapid for smaller f o . It is concluded from the observations made in the context of these figures for uniaxial stretch, as well as the earlier figures for hydrostatic deformation, that Version 3 of the second-order method leads to more consistent predictions for the overall behavior and microstructure evolution of the porous elastomers, and therefore it should be preferred over Version 1. For this reason, in the following sections only results associated with Version 3 will be presented.

PURE SHEAR LOADING

Figure 5 provides Version 3 second-order estimates for a Gent porous elastomer with incompressible matrix phase under pure shear ( λ1 = 1/ λ2 = λ). Results are shown for an initial porosity of 10% and values of the lock-up parameter J m = 42, 100, and J m → ∞ as functions of the logarithmic strain ē = ln( λ) for: (a) the effective stored-energy function compared with the Voigt upper bound; and (b) the evolution of the aspect ratio ω. First, it is observed from Figure 5(a) that the Version 3 estimates for the effective stored-energy function satisfy the rigorous Voigt upper bound. It is emphasized again that this bound is only helpful when considering isochoric loadings, like the one considered in this section, since it becomes unbounded otherwise. Note that no loss of ellipticity was detected at any level of deformation. In connection with the evolution of the microstructure, it is remarked that the porosity does not evolve under pure shear deformations. On the other hand, as clearly shown by Figure 5(b), the aspect ratio of the pores does increase fairly rapidly with increasing strains. Furthermore, note that ω appears to be insensitive to the value of the material parameter J m . Figure 6 presents plots of the corresponding results for the stress components: (a) S 11 ; and (b) S 22 . One of the main points that can be drawn from Figure 6 is the strong dependence of the effective stress-strain relation of the porous rubber on the lock-up parameter J m of the matrix phase. Interestingly, it can also be deduced from these figures that the evolution of the aspect ratio appears to have little effect on the effective constitutive behavior of the porous elastomer under pure shear. In order to aid the discussion of the results, the boundary at which the porosity vanishes has also been included in Figure 7. Note that once the zero-porosity boundary is reached, further compressive deformation of the material is not possible.

LOSS OF STRONG ELLIPTICITY

An interesting observation that can be made from Figure 7 is that the loci of points describing loss of strong ellipticity satisfy ē2 + ē1 < 0, which implies that a necessary condition for loss of strong ellipticity to occur is the existence of a compressive component in the state of deformation. Also, note that the predictions from both versions of the second-order method have roughly the same qualitative behavior; however, the results of Version 1 appear to be more restrictive than those of Version 3 for all initial values of porosity and loadings, with the exception of cases satisfying f o < 0.3, ē1 < 0, and ē2 < 0, for which the onset of loss of strong ellipticity of Version 3 precedes that one of Version 1. Furthermore, it is interesting (a) (b) to note that the strongly elliptic (and non-elliptic) domains shown in Figure 7 are similar to the results obtained by Abeyaratne and Triantafyllidis [START_REF] Abeyaratne | An investigation of localization in a porous elastic material using homogenization theory[END_REF] for the loss of strong ellipticity of periodic porous elastomers with a nearly incompressible Neo-Hookean matrix phase. However, their results appear to be more restrictive than the ones obtained here. In particular, these investigators did find loss of (strong) ellipticity for deformations with ē2 + ē1 > 0, which includes pure shear loading. These discrepancies seem to be consistent with their periodic microstructure, as it is more susceptible to instabilities than the random microstructure utilized in this work. Another important point that deserves further comment is the trend followed by the onset of loss of ellipticity as a function of initial porosity. In effect, Figure 7 suggests that a Neo-Hookean porous elastomer with random and isotropic microstructure becomes more stable with increasing value of initial porosity. This behavior is counterintuitive as one might expect an elastomer to be more unstable with increasing porosity. However, this is a complex and difficult problem, which will be pursued in future work. Finally it is interesting to remark that it was through the failure of the third (and equivalently fourth) condition of (71) that the porous elastomer with incompressible Neo-Hookean matrix phase lost strong ellipticity systematically. Indeed, whereas the evolution of the microstructure for compressive loadings led to the already-mentioned hardening of some of the components of the effective incremental modulus, it also led to the softening of the shear component L 1212 , which resulted in the overall loss of ellipticity of the porous elastomer.

For completeness, it is noted that the corresponding domains of strong ellipticity for porous elastomers with incompressible Gent matrix phases are essentially identical to those shown in Figure 7. Indeed, the results predicted by the secondorder theory indicate that the value of the lock-up parameter J m does not play a major role in estimating the onset of loss of ellipticity of porous elastomers with random and isotropic microstructures.

In summary, the second-order estimates for the homogenized constitutive behavior of porous elastomers with isotropic, strongly elliptic, matrix phases have been found to admit loss of strong ellipticity at reasonable levels of deformation. This behavior has been linked to the evolution of the microstructure under finite deformations, which, depending on the specific loading conditions, was found to induce hardening or softening behavior resulting in the loss of strong ellipticity for the porous elastomer.

Concluding Remarks

In this work, analytical estimates have been derived for the effective behavior of porous elastomers with random microstructure subjected to finite deformation, by means of an implementation of the second-order procedure of Lopez-Pamies and Ponte Castañeda [START_REF] Lopez-Pamies | Second-order estimates for the large-deformation response of particle-reinforced rubbers[END_REF][START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF]. It is emphasized that this homogenization technique, which is an extension of the variational method developed by Ponte Castañeda [START_REF] Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations[END_REF] in the context of viscoplastic materials, is applicable to a large class of hyperelastic composites including reinforced and porous rubbers.

A key issue in the general framework of the second-order variational procedure is the scheme employed for the linearization of the constitutive relation of the hyperelastic phases in the composite. In this regard, it has been seen that the earlier tangent linearization proposed by Ponte Castañeda and Tiberio [START_REF] Castañeda | A second-order homogenization procedure in finite elasticity and applications to black-filled elastomers[END_REF] results in estimates for the effective stored-energy function that depend exclusively on the average fields of the constituent phases. On the other hand, the estimates associated with the generalized secant linearization scheme not only depend on the average fields, but also exhibit a direct dependence on the field fluctuations. The difference between these two approaches has already been shown to be significant in the context of reinforced incompressible elastomers, where the incorporation of field fluctuations proved necessary to obtain the correct overall incompressibility constraint for these materials (see [START_REF] Lopez-Pamies | Second-order homogenization estimates incorporating field fluctuations in finite elasticity[END_REF]). Within the richer class of porous elastomers, the direct incorporation of field fluctuations into the computation of the effective behavior has turned out to be essential as well.

Thus, by incorporating field fluctuations, Version 3 of the second-order method has been shown to lead to the exact evolution of the porosity in porous elastomers with incompressible, isotropic, matrix phases, under general plane strain loading. This is a remarkable result in view of the strong nonlinearity of the problem. Furthermore, for the particular case of hydrostatic loading, the effective constitutive estimates delivered by Version 3 exhibit excellent agreement with the available exact result. This can be related to the correct prediction of the porosity evolution. Unfortunately, no other exact results are available for the effective constitutive behavior of porous elastomers. However, based on the comparisons presented, it seems plausible that Version 3 of the second-order variational procedure should be also able to deliver accurate estimates for the homogenized behavior of porous elastomers for more general loading conditions.

On the contrary, Version 1 of the second-order method, which only makes use of the average fields, delivers predictions for the evolution of the microstructure that deviate rapidly from the expected behavior for finite deformations, especially for tensile loadings. The negative consequences of this deviation were put in evidence by the comparisons with the exact result for hydrostatic loading, where the Version 1 estimates, even though exact to third order in the infinitesimal strain, break down under large tensile deformations.

A major result of this work is the strong influence of the microstructure evolution on the overall behavior of porous elastomers, in particular, through geometric hardening/softening mechanisms arising as a consequence of the evolution of the pore microstructure during a finite-deformation history. Indeed, it was seen that the decrease of the porosity during compressive deformations results in a significant hardening of the effective constitutive behavior of the porous elastomer. On the other hand, the increase of the porosity associated with tensile deformations leads to a pronounced softening.

Finally, it has been shown that loss of strong ellipticity, corresponding to the possible development of shear-band instabilities, can take place in porous elastomers with random microstructures at physically realistic levels of compressive deformation. This is consistent with earlier findings by Abeyaratne and Triantafyllidis [START_REF] Abeyaratne | An investigation of localization in a porous elastic material using homogenization theory[END_REF] for porous systems with periodic microstructures. Indeed, in this work, we have been able to relate softening mechanisms associated with the evolution of the microstructure under finite deformations with the possible onset of macroscopic instabilities, even for materials with strongly elliptic matrix phases. These encouraging results for two-dimensional microstructures should provide ample motivation to carry out corresponding analyses for porous and other types of elastomeric composites with more general three-dimensional, random microstructures, where comparisons with appropriate experimental results should be feasible.

Next, under condition (A.2), the equations of order O( 0 ) yield the relationship:

α 2 + α 2 1 β 2 = (1 + α 2 1 )(-1 -f o + (1 + f o )α 2 1 + ( λ2 -λ1 )α 1 )μ λ1 -α 2 1 λ2 , (A.3)
which determines the combination α 2 + α 2 1 β 2 in terms of α 1 . Finally, the equations of order O( ) derived from (49) are considered. Making use of relations (A.2) and (A.3) in these equations can be shown to lead to the following expressions:

α 3 + G 1 (α 1 , α 2 )β 3 = G 2 (α 1 , α 2 ) (A.4)
and

( λ2 2 + f 2 o -1)α 4 1 + 2( λ1 + (f o -1) λ2 )α 3 1 + ( λ2 2 -λ2 1 )α 2 1 -2((f o -1) λ1 + λ2 )α 1 + 1 -f 2 o -λ2 1 = 0, (A.5)
where G 1 and G 2 are (known functions of their arguments) too cumbersome to be included here.

It is noted that (A.4) establishes a linear relationship between β 3 and α 3 analogous to the one established by equation (A.3) between β 2 and α 2 . More importantly, (A.5) provides a fourth-order polynomial equation for the coefficient α 1 in terms of the initial concentration of pores f o , and the applied loading as determined by λ1 and λ2 . This equation is precisely the equation (51) given in the main body of the text, where for clarity of notation α 1 was denoted as u. It turns out that the leading order term of the effective energy (37) in the limit of incompressibility may eventually be characterized entirely in terms of the coefficient α 1 . The final result is given by expression (50) in the text, where, as already pointed out, u must be identified with α 1 .

It is noted that for the particular case of hydrostatic loading, i.e., λ2 = λ1 = λ, a suitable limit must be taken in the above expressions. For this type of deformation, it is straightforward to show that λ(1) 2 = λ(1) 1 , and hence that, β 1 = α 1 , β 2 = α 2 , and β 3 = α 3 . Now, making use of these relations together with the equation of order O( -1 ) given by (A.2) leads to α 1 = 1. In turn, this result for α 1 makes the equation (of order O( 0 )) (A.3) be satisfied trivially, whereas the one of order O( 1 ) can be shown to render the following identities:

α 2 = λ -1 λ + f o -1 μ, (A.6) α 3 = ( λ -1)(3 -5f o + 2f 2 o + (7f o -6) λ + 3 λ2 ) 2( λ + f o -1) 3 μ 2 .
Recognizing now that under expression (A.1), hydrostatic loading, and α 1 = 1, the expansion of the second-order estimate (37) in the incompressibility limit can be written, to first order, as

W I (F) = I ( λ, λ) = 2(1 -f o )( λ -1)α 2 + O( ), (A.7)
together with (A.6) 1 , leads to the final result (52).

Appendix B. Incompressible Limit for a Neo-Hookean Porous Elastomer (Version 3)

In this appendix, a brief outline of the asymptotic analysis corresponding to the incompressibility limit associated with the second-order estimate (39) for a porous elastomer with a Neo-Hookean matrix phase is presented. As discussed in the main body of the text, only one of the roots derivable from this version of the second-order method has a physically consistent asymptotic behavior in the limit of incompressibility. The limit associated with this root is the one presented here. It is noted that the results obtained from the following asymptotic analysis have been checked to be in agreement with the full numerical solution.

Based on numerical evidence from the results for general μ , an expansion is attempted in the limit as μ → ∞ of the following form: First, it is remarked that for the particular case of a Neo-Hookean matrix phase one of the generalized secant equations (41) can be solved exactly for the variable L 1212 in terms of the other components of the modulus L 0 . This, together with the constraints (56), can be shown to result into the following simplifications: Next, introducing relations (B.1) and (B.2) in the general expression (58) for the components of F (1) -F can be shown to lead to the following expansions: Prescriptions (B.4) through (B.8) can be shown to be sufficient to fully determine the first-order terms of all of the components of F (1) -F and F (1) . The final expressions may be written as 21 = 0. At this point, it is important to remark that relations (B.9) and (B.10), by means of (B.6), ultimately depend on the variable n 2 , which can be determined in closed-form by solving the fourth-order polynomial equation (B.7). This is precisely the same equation as (61) given in the main body of text, where for clarity of notation n 2 was relabelled as v. Under the above development, it is then straightforward to show that the leading order term of the expansion of the second-order estimate (39) in the limit of incompressibility may be expressed in closed-form, as it ultimately depends on the coefficient n 2 . The final explicit expression (in terms of the variable n 2 = v) is given by (60) in the text.

Next, it is shown that the porosity associated with the second-order estimate (60) for a porous elastomer with an incompressible Neo-Hookean matrix phase reduces to the exact result (70). Given that a HS-type approximation is utilized in the homogenization process, the fields in the porous phase are assumed constant. This implies that the average change in volume of the porous phase is simply given by J (2) = det(F) (2) = det( F (2) )

= (f o λ1 -(1 -f o )x 1 )(f o λ2 -(1 -f o )y 1 ) f 2 o , (B.11)
where use has been made of the relation F = (1f o )F (1) + f o F (2) . Expression (B.11) can now be used to compute the porosity associated with the secondorder estimate (60) through the relation

f = J (2) J f o , (B.12)
which, after some simplification, can be shown to reduce to the exact result (70). Finally it should be emphasized that this result has been proven to hold not only for Neo-Hookean porous elastomers, but more generally, for porous elastomers with incompressible isotropic matrix phases. Note that the factors of μ in the above expressions have been included for consistency with the results from Appendix B. However, these factors cancel out in equations ( 60) and (61), and therefore μ may be dropped from the above expressions.

r 1 =

Figure 1 .

 1 Figure 1. Comparisons of the second-order estimates (Versions 1 and 3) with the exact results for the effective response of a porous rubber subjected to in-plane hydrostatic loading ( λ2 = λ1 = λ). The results correspond to an incompressible Neo-Hookean matrix phase with various initial concentrations f o of aligned cylindrical voids, and are shown as a function of the logarithmic strain ē = ln( λ). (a) The stored-energy function W ; and (b) the corresponding stresses S 11 = S 22 = S.

Figure 2 .

 2 Figure 2. In-plane hydrostatic loading ( λ2 = λ1 = λ) of a porous rubber with an incompressible, Neo-Hookean matrix phase with various initial concentrations f o of aligned cylindrical voids. (a) The evolution of porosity f as predicted by Versions 1 and 3 of the second-order method compared with the exact result as a function of the logarithmic strain ē = ln( λ). (b) The critical stretches λcrit at which the loss of strong ellipticity of the homogenized elastomer takes place as a function of initial porosity f o . (This last plot also includes the critical loads for the first two buckling modes (n = 2 and 3) of a cylindrical shell [34].)

Figure 3 .

 3 Figure 3. Versions 1 and 3 estimates of the second-order method for the effective response of a porous rubber subjected to uniaxial loading ( λ2 = 1 and λ1 = λ). The results correspond to an incompressible Neo-Hookean matrix phase with various initial concentrations f o of aligned cylindrical voids, and are shown as a function of the logarithmic strain ē = ln( λ). (a) The stress component S 11 . (b) The stress component S 22 .

Figure 4 .

 4 Figure 4. Versions 1 and 3 estimates of the second-order method for the effective response of a porous rubber subjected to uniaxial loading ( λ2 = 1 and λ1 = λ). The results correspond to an incompressible Neo-Hookean matrix phase with various initial concentrations f o of aligned cylindrical voids, and are shown as a function of the logarithmic strain ē = ln( λ). (a) The evolution of the porosity f compared with the exact result. (b) The evolution of the average aspect ratio of the voids ω.

Figure 5 .

 5 Figure 5. Version 3 estimates of the second-order method for the effective response of a porous rubber subjected to pure shear loading ( λ1 = 1/ λ2 = λ). The results correspond to an incompressible Gent matrix phase with given initial porosity f o = 0.1 and various values of the material parameter J m , and are shown as a function of the logarithmic strain ē = ln( λ). (a) The effective stored-energy function W compared with the Voigt upper bound. (b) The evolution of the aspect ratio ω.

Figure 6 .

 6 Figure 6. Version 3 estimates of the second-order method for the effective response of a porous rubber subjected to pure shear loading ( λ1 = 1/ λ2 = λ). The results correspond to an incompressible Gent matrix phase with given initial porosity f o = 0.1 and various values of the material parameter J m , and are shown as a function of the logarithmic strain ē = ln( λ). (a) The stress component S 11 . (b) The stress component S 22 .

Figure 7

 7 Figure 7 displays the strongly elliptic (and non-elliptic) domains for the 2-D porous elastomer with incompressible Neo-Hookean matrix phase, subjected to in-plane deformations. The results are shown in the plane ( ē1 -ē2 ) for: (a) Versions 1 and 3 estimates for initial porosities f o = 30, 50 and 70%; and for (b) Version 3 estimates for initial porosities f o = 10, 20 and 30%.In order to aid the discussion of the results, the boundary at which the porosity vanishes has also been included in Figure7. Note that once the zero-porosity boundary is reached, further compressive deformation of the material is not possible.An interesting observation that can be made from Figure7is that the loci of points describing loss of strong ellipticity satisfy ē2 + ē1 < 0, which implies that a necessary condition for loss of strong ellipticity to occur is the existence of a compressive component in the state of deformation. Also, note that the predictions from both versions of the second-order method have roughly the same qualitative behavior; however, the results of Version 1 appear to be more restrictive than those of Version 3 for all initial values of porosity and loadings, with the exception of cases satisfying f o < 0.3, ē1 < 0, and ē2 < 0, for which the onset of loss of strong ellipticity of Version 3 precedes that one of Version 1. Furthermore, it is interesting

Figure 7 .

 7 Figure 7. Domains of strong ellipticity on the ( ē1 -ē2 )-plane for a porous elastomer with incompressible, Neo-Hookean, matrix phase and various levels of initial concentrations f o of aligned cylindrical voids, as determined by Versions 1 and 3 of the second-order variational procedure. The dotted lines denote the boundary at which the level of zero porosity has been reached upon compressive deformation. (a) Comparisons between the Versions 1 and 3 estimates; and (b) Version 3 estimates for low initial porosity.

L 1111 = a 1 +

 1 a 2 + a 3 + O( 2 ), L 2222 = b 1 + b 2 + b 3 + O( 2 ), L 1122 = c 1 + c 2 + c 3 + O( 2 ), (B.1)L 1212 = d 1 + d 2 + O( ), L 1221 = e 1 + e 2 + O( ),where . = 1/μ is a small parameter and a 1 ,a 2 , a 3 , b 1 , b 2 , b 3 , c 1 , c 2 , c 3 , d 1 , d 2 , e 1, and e 2 are unknown coefficients that ultimately depend on the applied loading F, the initial concentration of voids f o , and the material parameter μ.

d 1

 1 = 0, d 2 = μ, e 1 = a 1 b 1c 1 , and (B.2) e 2 = a 2 b 1 + a 1 b 2 -(a 1 + b 1 )μ 2 √ a 1 b 1 c 2 .

1 )

 1 + λ2 ( λ1 λ2 -1)((f o -1)n 2 + 2μ( λ1 + λ2 ))] a 1 [2(f o -1)n 2 λ1μ( λ1 + λ2 )((f o -3) λ1 -(1 + f o ) λ2 )],p 1 = x 1 y 1 + λ2 x 1 + λ1 y 1 + λ1 λ2 -1, (B.9)s 1 = μ( λ1 -λ2 ) 2 ( λ1 + λ2 ) 2 a 2 1 f o λ1 [2(f o -1)n 2 λ1μ( λ1 + λ2 )((f o -3) λ1 -(1 + f o ) λ2 )] 2 × a 1 f o μ( λ1 + λ2 ) 2 + (f o -1) λ2 ( λ1 λ2 -1) n 2μ( λ1 + λ2 ) × 2a 1 f o λ1 + λ2 (1 + f o -(1 + f o ) λ1 λ2 )

  1) 2 -λ2 )(L 1122 P 1111 + L 2222 P 1122 ) + ( λ(1) 1 -λ1 ) (L 1122 P 1122 -1)(1 + (f o -1)L 1122 P 1122 ) -(f o -1)(L 2222 + L 2 1122 P 1111 -L 1111 L 2222 P 1111 )P 2222 + L 1111 (P 1111 -(f o -1)L 2222 P 2 1122 ) f o P 1111 -(f o -1)L 2222 P 2 1122 + (f o -1)L 2222 P 1111 P 2222 S -(f o -1)(L 1111 + L 2 1122 P 2222 -L 2222 L 1111 P 2222 )P 1111 + L 2222 (P 2222 -(f o -1)L 1111 P 2 1122 ) f o P 2222 -(f o -1)L 1111 P 2 1122 + (f o -1)L 1111 P 2222 P 1111 S

	(1)
	11
	(1) 22 = 0
	and
	f o ( λ(1)
	(1)
	22

f o P 1122 + (f o -1)L 1122 P 2 1122 -(f o -1)L 1122 P 1111 P 2222 S 1 -λ1 )(L 1122 P 2222 + L 1111 P 1122 ) + ( λ(1) 2 -λ2 ) (L 1122 P 1122 -1)(1 + (f o -1)L 1122 P 1122 ) f o P 1122 + (f o -1)L 1122 P 2 1122 -(f o -1)L 1122 P 2222 P 1111 S (1)

  4μ 3 λ1 ( λ1 + λ2 ) (-1 + f o ) 3 λ4 + f o ) 4 -(1 + f o (-4 + f o + f 2 + λ2 2f o (4 -4f o + f 3 o + (-3 + f o ) λ2 2 )) , r 0 = -μ 4 ( λ1 + λ2 ) 2 (-1 + f o ) 4 λ4 1 -(-1 + f o ) 2 λ3 1 (2f o (1 + f o ) -(-1 + f o ) λ2 1 ) λ2 + λ2 1 (-2 + f o (2 + 3f o + f 3 o -2(-3 + f 2

		1 λ2 + λ1	λ2 2 (-1 + f 2 o
	+ (1 + f o + f 2 o ) λ2 2 ) + λ3 1 ((-1 o )) λ2 2 )
	+ λ3 2 (-1 + f o (f o + λ2 2 ))
	+ λ2 1 λ2 (1 o ) λ2 1 )) λ2 2
	+ λ1 (2f o (1 + f o ) 2 + (2 + f o (4 + f o + f 2 o )) λ2 1 ) λ3 2
	+ (1 + f o )(1 + f o + 2f o	λ2 1 ) λ4 2 + (-1 + f o ) λ1	λ5 2 .

Acknowledgement

This work was supported by NSF grant DMS-0204617.

Appendix A. Incompressible Limit for a Neo-Hookean Porous Elastomer (Version 1)

In this appendix some details are presented concerning the incompressibility limit associated with the tangent second-order estimate (37) for a porous elastomer with Neo-Hookean matrix. The asymptotic solution resulting from this heuristic derivation has been checked to be in agreement with the full numerical results.

Motivated by the observed properties of the numerical solution for general μ , an expansion is attempted in the limit as μ → ∞ of the following form:

where . = 1/μ is a small parameter, and α 1 , α 2 , α 3 , β 1 , β 2 , and β 3 are unknown coefficients which ultimately depend on the applied loading F, the initial concentration of voids f o , and the material parameter μ.

By making use of expressions (A.1) in relation (49), a hierarchical system of equations is obtained for the coefficients α 1 , α 2 , α 3 , β 1 , β 2 , and β 3 . The leading order terms O( -1 ) of these equations can be shown to lead to the following relationship:

which implies that the determinant of F (1) , denoted by J (1) , is exactly equal to one in the incompressible limit.

F (1) 11 -λ1 = x 1 + x 2 + O( 2 ), F (1) 22 -λ2 = y 1 + y 2 + O( 2), (B.3) F (1) 12 F (1) 21 = p 1 + p 2 + O( 2), ( F (1) 12 ) 2 + ( F (1) 21

The explicit expressions for the coefficients of these expansions have not been included here for their bulkiness; however, it is useful to spell out their dependence on the variables introduced in (B.1). Thus, the coefficients of first order x 1 , y In connection with relations (B.3), it is necessary to clarify that the asymptotic expressions for the combinations F (1) 12 F (1) 12 and ( F (1) 12 ) 2 + ( F (1) 12 ) 2 have been specified in (B.3), rather than those for the independent components F (1) 12 and F (1) 12 , since, as discussed previously, they are the relevant variables in this problem. Now, by introducing expressions (B.1)-(B.3) into the three reduced (recall that L 1212 = μ) generalized secant equations (41), a hierarchical system of equations is obtained for the remaining unknown coefficients introduced in (B.1). Thus, the equations of first-order O( -1 ) lead to the following results:

whereas the equations of second-order O( 0 ), by making use of (B.4), can be shown to render the following relations:

where n 2 = λ1 a 2 -λ2 c 2 ,

and q 2 , q 1 , q 0 , r 4 , r 3 , r 2 , r 1 , and r 0 have been given in explicit form in Appendix C.

Appendix C. Coefficients Associated with the Incompressible Limit for a Neo-Hookean Porous Elastomer (Version 3)

In this appendix, the expression for the coefficients introduced in relations ( 60) and (61) are given in explicit form in terms of λ1 , λ2 , f o , and μ: