A Simple Explicit Model Constructed from a Homogenization Solution for the Large-Strain Mechanical Response of Elastomeric Syntactic Foams

This paper introduces a homogenization-based constitutive model for the large-strain mechanical response of elastomeric syntactic foams subject to arbitrary quasistatic loading and unloading conditions. Based on direct observations from experiments, this class of emerging foams are considered to be made of a nonlinear elastic matrix filled with a random isotropic distribution of hollow thin-walled spherical shells — commonly termed microspheres or microballoons — each having the same mean diameter that are made of an elastic brittle material that is much stiffer than the elastomeric matrix, typically, either glass or a hard polymer. Accordingly, such underlying microballoons behave effectively as rigid particles initially. Along a given loading path, however, they may fracture (in the case of glass microballoons) or buckle (in the case of polymer microballoons) at which point they abruptly transition to behave effectively as vacuous pores. On that account, the proposed constitutive model corresponds to a homogenization solution for the nonlinear elastic response of particle-filled porous elastomers — precisely, elastomers embedding both rigid spherical particles and vacuous spherical pores of equal monodisperse size — wherein the volume fraction of pores corresponds to the volume fraction of fractured or buckled microballoons and hence is not a fixed parameter but rather an internal variable of state that evolves as a function of the loading history. After the general presentation of the model, where its theoretical and practical features are discussed, its descriptive and predictive capabilities are showcased via comparisons with experimental data for silicone syntactic foams filled with glass microballoons.


Introduction
The origin of syntactic foams, namely, materials comprised of a metal, polymer, or ceramic matrix filled with hollow spheres -typically made of a ceramic or a hard polymer and commonly termed microspheres or microballoons -of micron or sub-micron size, dates back to the 1950s. Since their discovery, they have increasingly proven very effective in advancing a wide range of technologies and as a result have occupied the interest of numerous investigators in a plurality of fields; see, e.g., Anderson (1970), Rand (1978), Shutov (1986), John and Nair (2010), Gupta et al. (2014).
Despite early indications of their potential to enable new technologies (Rand, 1978), investigations of soft elastomeric syntactic foams have not been pursued to nearly the same extent as those of stiffer counterparts, presumably because of the technical challenges associated with finite deformations. Indeed, although the linear elastic response of syntactic foams has been studied quite extensively previously -see, e.g., Bardella and Genna (2001), Gupta and Woldesenbet (2004), Porfiri and Gupta (2009), and Tagliavia et al. (2009) -only recently have experiments been reported on the macroscopic mechanical response of elastomeric syntactic foams at finite deformations (Brown et al., 2018;Yousaf et al., 2020) and on the associated fracturing and buckling of the underlying microballoons (Shorter et al., 2008;Croom et al., 2019a;2019b;Yousaf et al., 2020). Theoretical models to describe the buckling of plastic shells were described in Fok and Allwright (2001), Jones (2008), and Thorpe (2016). The only prior model associated with finite deformations is the micromechanical model introduced by De Pascalis et al. (2014), which deals exclusively with the case of the prediction of the nonlinear loading curve for hydrostatic deformations.
In this context, the objective of this paper is to introduce a constitutive model aimed at explaining, describing, and predicting the mechanical response of elastomeric syntactic foams under arbitrarily large quasistatic deformations directly in terms of their microscopic properties and, in particular, of the fracturing or buckling of their underlying microballoons. The focus is on the prominent class of isotropic elastomeric syntactic foams filled with glass microballoons or microballoons made of a hard polymer.
As elaborated in Section 3, based on direct observations from experiments, one of two central ideas of the proposed model is to consider that the microballoons in elastomeric syntactic foams behave initially as rigid particles. Along a given loading path, however, they may fracture (in the case of glass microballoons) or buckle (in the case of polymer microballoons). At soon as they do so, they are assumed to behave effectively as vacuous pores. The second central idea, also based on direct observations from experiments, is to consider that at any given state of loading, elastomeric syntactic foams contain a statistically uniform distribution of a mixture of intact and fractured/buckled microballoons. In other words, the fracturing or buckling of microballoons occurs uniformly throughout the foams, rather than occurring in isolated localized regions.
Granted the above-two premises, it follows that the mechanical response of elastomeric syntactic foams can be identified as the homogenized response of a class of particle-filled porous elastomers -precisely, elastomers embedding both rigid spherical particles and vacuous spherical pores of equal monodisperse size 1 -wherein the volume fraction of pores corresponds to the volume fraction of fractured or buckled microballoons and hence is not a fixed parameter but rather an internal variable of state that evolves as a function of the loading history. For clarity of exposition, the relevant homogenization solution is presented in Section 2, prior to its summons in Section 3.
Aimed at illustrating the use of the proposed model as well as at showcasing its descriptive and predictive capabilities, Sections 4 and 5 are devoted to presenting comparisons with the experimental results of Brown et al. (2016; for silicone syntactic foams filled with glass microballoons and to placing on record some concluding remarks.

A homogenization solution for the nonlinear elastic response of particle-filled porous elastomers
As anticipated in the Introduction and as will become apparent in the next section, an essential result needed in the construction of the model proposed in this work for elastomeric syntactic foams is the homogenization solution for the nonlinear elastic response at finite deformations of a specific class of particle-filled porous elastomer.
The microscopic description of the particle-filled porous elastomers. Precisely, the class of particle-filled porous elastomers of interest here consists of a random isotropic distribution of rigid spherical particles firmly bonded to a homogeneous nonlinear elastic matrix that also contains a random isotropic distribution of vacuous spherical pores; such microstructures are reminiscent of those found in a variety of other material systems, for instance, concrete, cortical bone, and bamboo (Giorgio et al., 2019;Tal and Fish, 2018;Mannan et al., 2017;Contrafatto et al., 2016). The domain occupied by any such three-phase composite material in its undeformed and stress-free ground state is denoted by Ω 0 and its boundary by ∂Ω 0 . All particles and pores are taken to have an identical initial diameter that is much smaller than the length scale of Ω 0 ; see Fig. 1 for an schematic.

• •
Figure 1: Schematic of the particle-filled porous elastomer of interest in its ground configuration Ω 0 depicting its underlying microstructure comprised of a random isotropic distribution of rigid spherical particles and vacuous spherical pores. All particles and pores have identical diameter. The nonlinear elastic behaviors of the elastomeric matrix, the rigid particles, and the vacuous pores are characterized by stored-energy functions Wm, Wr, and Wp. The macroscopic or homogenized nonlinear elastic behavior of the particle-filled porous elastomer is characterized by the effective stored-energy function W .
The nonlinear elastic behavior of the matrix (m) is characterized by the isotropic and incompressible storedenergy function where I 1 = F · F = tr C, J = det F = √ det C, with C = F T F, and where Ψ m stands for any non-negative function of choice satisfying the linearization conditions Ψ m (3) = 0 and Ψ m (3) = dΨ(3)/dI 1 = µ m /2, with µ m denoting the initial shear modulus of the matrix, and the strong ellipticity conditions Ψ m (I 1 ) > 0 and Ψ m (I 1 ) + 2(I 1 − λ 2 i − 2λ −1 i )Ψ m (I 1 ) > 0 (i = 1, 2, 3), with λ 1 , λ 2 , λ 3 denoting the singular values of F; see, e.g., Section 4.2.7 in the monograph by Ogden (1997) for an in-depth discussion of the algebraic manipulation of stored-energy functions subject to constraints.
Examples for the function Ψ m include, for instance, the Neo-Hookean model introduced by Treloar (1943), the model introduced by Lopez-Pamies (2010), as well as the Arruda and Boyce (1993) and Gent (1996) models among others. In this last expression, µ 1 > 0, µ 2 ≥ 0, α 1 , α 2 are realvalued material parameters, the first two of which satisfy the condition µ 1 + µ 2 = µ m . The rigid (r) and vacuous constitutive behaviors of the particles and the pores (p) are characterized by the storedenergy functions and W p (F) = 0; see, e.g., Chi et al. (2016) for an in-depth discussion of the algebraic manipulation of stored-energy functions subject to the rigidity constraint (4). At any material point X ∈ Ω 0 , it follows that the first Piola-Kirchhoff stress tensor S is given by where θ r and θ p are the characteristic functions describing the spatial locations occupied by the particles and the pores in Ω 0 , i.e., θ r takes the value of 1 if the position vector X lies within a particle and zero otherwise, while θ p takes the value of 1 if the position vector X lies within a pore and zero otherwise. For later use, we introduce the notation and c = c r + c p to express the volume fractions of the particles and the pores, as well as the combined total volume fraction of particles and pores in the initial domain Ω 0 occupied by the particle-filled porous elastomer.
The macroscopic or homogenized response. Granted the isotropic (and hence statistically uniform) distribution of the particles and the pores and their vanishingly small size relative to the macroscopic size of Ω 0 -taken henceforth to denote a representative volume element -the macroscopic or homogenized constitutive response of the particle-filled porous elastomer can be formally defined as the relation between the volume average of the first Piola-Kirchhoff stress S := |Ω 0 | −1 Ω0 S(X)dX and the volume average of the deformation gradient F := |Ω 0 | −1 Ω0 F(X)dX under affine boundary conditions on ∂Ω 0 (Hill, 1972).
It follows from standard arguments that such a relation can be expediently cast in the variational form where W , termed the effective stored-energy function, represents physically the total elastic energy (per unit undeformed volume) stored in the material; in this last expression, K stands for a sufficiently large set of admissible deformation gradients F(X) with prescribed volume average F. It further follows from the overall local geometric and constitutive isotropy of the material at hand that the effective stored-energy function W defined by the variational problem (6) admits the representation in terms of the standard principal invariants In a string of recent contributions, Lopez-Pamies et al. (2013a; 2013b), Goudarzi et al. (2015), Lefèvre and Lopez-Pamies (2017a;2017b), and Shrimali et al. (2019) have derived results that include as special cases the solution for the homogenization problem (6) when the solid is filled only with rigid spherical particles (i.e., when θ p = 0) or only with vacuous spherical pores (i.e., when θ r = 0). In the next two subsections, we leverage the same computational and analytical techniques used in those works to generate the more general solution for the homogenization problem (6) when rigid spherical particles and vacuous spherical pores are both present (i.e., when θ r = 0 and θ p = 0).

Computational homogenization
In the footstep of a now well-established practice in infinitesimal and finite elasticity (Gusev, 1997;Michel et al., 1999;Segurado and Llorca, 2002;Ghossein and Levesque, 2012;Lopez-Pamies et al., 2013b;Lefèvre and Lopez-Pamies, 2017a;Shrimali et al., 2019), the particle-filled porous elastomers studied in this work can be accurately approximated as infinite media made out of the periodic repetition of a unit cell containing a random distribution of a sufficiently large but finite number N = N r + N p of N r rigid spherical particles and N p vacuous spherical pores.
For expediency in the implementation of periodicity, we select here the defining unit cell to be a unit cube Y 0 = (0, L) 3 = (0, L) × (0, L) × (0, L) with edges of length L. To be able to span the large range of combined total volume fractions of particles and pores c = c r +c p ∈ [0, 0.5] of relevance for elastomeric syntactic foams, we make use of the algorithm introduced by Lubachevsky and Stillinger (1990); see also Lubachevsky et al. (1991). Roughly speaking, the idea behind this algorithm is to randomly seed at once in the unit cell Y 0 the desired total number N of particles and pores as points endowed with random velocities and a uniform radial growth rate. As the points move and grow into spheres, their collision with one another are described by conservation of momentum, while their crossings through the boundaries of the unit cell are described by periodicity. When the desired volume fraction c is reached, the algorithm is stopped. The final step consists in randomly selecting which N r of the N spheres are rigid spherical particles and which N p of them are vacuous spherical pores. For a desired volume fraction c r of rigid spherical particles and a desired volume fraction c p of vacuous spherical pores, where R = 3L 3 c 4πN 1/3 stands for the common radius of the particles and pores. Clearly, for a given pair N and c, relations (7) may not lead to integer values for N r and N p in terms of the desired values of c r and c p , in which case the latter target pair may only be approximated. For definiteness, all the computational homogenization results presented in this work correspond to unit cells that contain a total of N = 120 spheres. Accordingly, in terms of the combined volume fraction c of particles and pores, their radius is given by R = L(3c/4πN ) 1/3 = 0.1258Lc 1/3 . We emphasize that the number N = 120 of spheres in the unit cells was carefully selected based on a parametric study that confirmed that it is sufficiently large to render overall nonlinear elastic behaviors for the materials of interest here that are indeed approximately isotropic; see Remark 1 in Leonard et al. (2020). By way of an example, Fig. 2 shows two representative unit cells at the same volume fraction c = c r + c p = 0.2 of particles and pores with two different combinations of volume fractions c r and c p .

FE discretization
The computations for the homogenized response (5)-(6) of the particle-filled porous elastomers with the aboveoutlined microstructures can be expediently carried out with the finite-element (FE) method.
In this work, we make use of the open-source mesh generator code Netgen (Schöberl, 1997) to discretize the constructed unit cells with non-overlapping 10-node tetrahedral elements. Meshes with about 1.75 million elements were checked to be sufficiently refined to deliver accurate solutions.
Because of the incompressibility of the matrix material, we make use of a hybrid re-formulation of the variational problem (6) in which both the displacement field u and a pressure field p are the independent fields in the problem. That re-formulation also restricts attention to solutions that are periodic of period the unit cell Y 0 ; see Section 5 in Lefèvre and Lopez-Pamies (2017b) for the relevant details. Within the hybrid re-formulation, we make use of tetrahedral elements which feature approximations that are quadratic in the displacement field and linear in the pressure field. The numerical implementation of the relevant discretized equations and the computation of their solution are carried out in the commercial FE code ABAQUS.
As a final note in this subsection, we remark that all the FE results that are presented in the sequel correspond to the average of two realizations. For each type of particlefilled porous elastomer considered, the responses of all two corresponding realizations showed very small differences (less than 3%) between one another. Figure 3 presents FE solutions (solid circles) for the macroscopic elastic response in the limit of small deformations of particle-filled porous elastomers with total volume fractions of particles and pores c = c r +c p = 0.1, 0.2, 0.3, 0.4, and 0.5; details of the FE calculations involved in the asymptotic limit of small deformations can be found in Appendix A of Spinelli et al. (2015).

FE solutions for small deformations
Specifically, plots are presented for the normalized values of the effective shear modulus µ/µ m in Fig. 3(a) and of the effective bulk modulus κ/µ m in Fig. 3(b) as functions of the volume fraction of pores c p . The plots also include the corresponding results (solid lines) from the explicit approximation proposed below. Those are given by the formulae (10) and are clearly in good agreement with the FE results for any combination of volume fractions of particles and pores in the entire range c = c r + c p ∈ [0, 0.5] considered here.
Remark 1. Particle/pore size polydispersity. The effective moduli (10) are exact results for particle-filled porous elastomers with two opposite limiting classes of microstructures: i) porous elastomers wherein the matrix material contains rigid particles that are much smaller than the pores and ii) particle-filled elastomers wherein the matrix material contains vacuous pores that are much smaller than the rigid particles. The fact that the formulae (10) fare so well in Fig. 3 with the FE solutions for particlefilled porous elastomers wherein the rigid particles and the vacuous pores are of the same size indicates that the size disparity between particles and pores has little effect on the macroscopic response of these material systems all the way up to c = 0.5, which is already remarkably close to the percolation threshold c ≈ 0.64 (Scott, 1960;Lubachevsky et al., 1991;Torquato and Stillinger, 2010).

Sample FE solutions for finite deformations
Figures 4 and 5 present representative FE solutions (solid circles) for the macroscopic elastic response of particlefilled porous elastomers at finite deformations. The results pertain to materials wherein the underlying elastomeric matrix is characterized by the basic Neo-Hookean storedenergy function (2) and the total volume fraction of particles and pores is c = 0.2 corresponding to two different combinations of particles and pores, c r = 0.1, c p = 0.1 and c r = 0.05, c p = 0.15.
Specifically, Fig. 4 presents plots of the normalized effective stored-energy function W /µ m in terms of the invariant I 1 in part (a) and of the invariant I 2 in part (b) for fixed values of the remaining two invariants, I 2 = 3.63, J = 1.08 and I 1 = 3.59, J = 1.17, respectively. Figure  5 presents plots of W /µ m in terms of the invariant J for the fixed values of the remaining two invariants I 1 = 4.5, I 2 = 5.55 and I 1 = 4.25, I 2 = 5.05 and the volume fraction of pores c p = 0.1 in part (a) and c p = 0.15 in part (b). Here, it is fitting to recall that fixing the values of any pair of invariants I 1 , I 2 , J restricts the range of physically allowable values of the remaining invariant. The results displayed in Figs. 4 and 5 pertain to the entire range of allowable values for each of the cases that is presented. Now, it is plain from Figs. 4 and 5 that the computed effective stored-energy function (6) for the particle-filled porous elastomers of interest in this work: • depends roughly linearly on I 1 , • is essentially independent of I 2 , • depends nonlinearly on J, and • is approximately of the separable form when the underlying matrix material is Neo-Hookean. These are precisely the same four key features identified by Lopez-Pamies et al. (2013a;2013b), Goudarzi et al. (2015), Lefèvre and Lopez-Pamies (2017a;2017b), and Shrimali et al. (2019) for filled and porous Neo-Hookean elastomers. Following then the same two-step approach put forth in those works 2 , one can readily construct a simple explicit approximate homogenization solution for the effective storedenergy function W of particle-filled porous elastomers with any non-Gaussian -not just Neo-Hookean -matrix material. We present such an approximate solution next.

An explicit approximate homogenization solution
As alluded to above, in view of the functional features of W revealed by the FE solutions, it is not difficult to deduce that the computational homogenization results presented in the preceding subsection for the effective storedenergy function of the particle-filled porous elastomers can be described by the explicit approximation Before proceeding with the use of the result (8) in the next section to model the mechanical response of elastomeric syntactic foams, it is appropriate to spell out a number of its features.
Remark 2. The macroscopic constitutive response. The macroscopic constitutive relation (5) implied by the effective stored-energy function (8) reads as Remark 3. The limit of small deformations. In the limit of small deformations as F → I, the effective stored-energy function (8) reduces asymptotically to 10) stand for the effective shear and bulk moduli of the particlefilled porous elastomer.
In the further limit of dilute volume fraction of particles as c r 0, the effective moduli (10) reduce to leading order to which are nothing more than the classical Hashin-Shtrikman bounds for the shear and bulk moduli of an isotropic porous elastic material with incompressible matrix (Hashin and Shtrikman, 1963); it is important to recall here that the Hashin-Shtrikman bounds (10) are realizable by a multitude of microstructures. On the other hand, in the opposite limit of dilute volume fraction of pores as c p 0, the effective moduli (10) reduce to leading order to the classical exact result due to Roscoe (1973) for an isotropic incompressible linear elastic material reinforced by an isotropic distribution of rigid spherical particles of infinitely many diverse sizes, namely, µ = 1 (1 − c r ) 5/2 µ m and κ = +∞.
Remark 4. The dilute limit of particles. In the limit as the volume fraction of rigid spherical particles c r 0, the effective stored-energy function (8) for the effective stored-energy function of porous non-Gaussian elastomers, cf. equation (5) in that reference.
Remark 5. The dilute limit of pores. In the limit as the volume fraction of vacuous spherical pores c p 0, the effective stored-energy function (8) reduces asymptotically to For isochoric deformations when J = 1, the leading order term in (11) agrees identically with the result of Lopez-Pamies et al. (2013b) for the effective stored-energy function of particle-filled non-Gaussian elastomers, cf. equation (2) in that reference. For non-isochoric deformations when J > 1, on the other hand, the leading order term in (11) differs from the result of Lopez-Pamies et al. (2013b) in that it is finite. This is because under loading conditions with sufficiently large ratios of hydrostatic-to-shear stresses, the initially zero-volume pores can "cavitate" and grow elastically to occupy a finite volume within the material; see, e.g., Lopez-Pamies et al. (2011a;2011b).
Remark 6. Closure of the pores. The effective storedenergy function (8) features the following asymptotic behavior: This is the macroscopic manifestation of the purely kinematical fact that the current porosity in the particle-filled porous elastomer vanishes, that is, complete closure of the underlying pores ensues, whenever the determinant of the macroscopic deformation gradient F reaches the critical value J = 1 − c p , at which point the particle-filled porous elastomer (because of the incompressibility of the underlying matrix and rigid particles) behaves as an incompressible solid for loadings with any further volumetric compression.
Remark 7. The case of Gaussian or Neo-Hookean matrix. For the basic case when the matrix is Neo-Hookean, namely, when the function Ψ m is given by (2), the effective stored-energy function (8) specializes to In turn, the associated macroscopic response (9) specializes to Remark 8. Accuracy at small and finite deformations.
In spite of its simple and explicit form, as illustrated by the direct comparisons with computational homogenization results presented in Figs. 3 through 5 above and as demonstrated by a plurality of other comparisons not included here, the approximation (8) describes fairly accurately the macroscopic nonlinear elastic response of the particle-filled porous elastomers of interest here at small as well as at finite deformations, this for Gaussian and non-Gaussian behaviors of the underlying matrix and for any combination of volume fractions of particles and pores in the entire range c = c r + c p ∈ [0, 0.5] of practical relevance for elastomeric syntactic foams.

The proposed model for elastomeric syntactic foams
As outlined in the Introduction, the microballoons in elastomeric syntactic foams are typically made of either glass or a hard polymer, both of which are much stiffer than the embedding elastomeric matrix. As a result, they behave effectively as rigid spherical particles initially. Along a given loading path, however, the experiments show that the glass microballoons may fracture while the polymer microballoons may buckle. The experiments also show that the fracturing or the buckling of the microballoons tend to occur uniformly throughout the entire specimen being tested, as opposed to this occurring in a single localized region; see, e.g., Mae et al. (2008), Croom et al. (2019a;2019b), Lu et al. (2019), and Yousaf et al. (2020). This evidence prompts the following two idealizations: • after fracturing or buckling, the microballoons in elastomeric syntactic foams transition to behave effectively as vacuous pores and • at any given state of loading, elastomeric syntactic foams contain a statistically uniform distribution of a mixture of intact and fractured/buckled microballoons.
Granted these two idealizations, it is plain that the mechanical response of elastomeric syntactic foams can be described by the homogenization solution worked out in the preceding section for particle-filled porous elastomers with the caveat that the volume fraction of pores c p corresponds now to the volume fraction of fractured/buckled microballoons and hence is not a fixed parameter, but, instead, a variable that evolves as a function of the loading history.
In the sequel, we make the above statement precise. We do so by considering c p to be an internal variable within the classical constitutive framework based on internal variables of state; see, e.g., Coleman and Gurtin (1967), Germain et al. (1983), Maugin (2015), and Kumar and Lopez-Pamies (2016). Within that framework, absent changes in temperature, the constitutive relation describing the macroscopic mechanical response of elastomeric syntactic foams can be formally written in the compact form where S once again denotes the first Piola-Kirchhoff stress tensor, ψ = ψ(F, c p ) stands for the free-energy function, and E(F, c p ) is the evolution equation for the internal variable c p , taken henceforth, again, to correspond to the volume fraction of fractured/buckled microballoons. Given its physical meaning, c p is bounded from below and above according to 0 ≤ c p ≤ c, where c stands for the volume fraction of microballons in the syntactic foam at hand. We proceed in Subsections 3.1 and 3.2 by presenting specific prescriptions for the free-energy function ψ(F, c p ) and the evolution equation E(F, c p ) = 0 and then spell out the constitutive response (13) that they imply in Subsection 3.3 together with some of its salient features.

The free-energy function ψ(F, c p )
Upon transcribing the homogenization solution (8) to the current setting, we prescribe the free-energy function in (13) to be given by Here, we emphasize that c stands for the volume fraction of microballoons and hence is a fixed parameter for any given elastomeric syntactic foam. We also recall that the function Ψ m denotes the finite branch of the storedenergy function (1) describing the isotropic and incompressible elasticity of the underlying elastomeric matrix material, while, again, I 1 = F · F and J = det F.

The evolution equation E(F, c p ) = 0 for c p
In principle, the evolution equation for the volume fraction c p of fractured/buckled microballoons in elastomeric syntactic foams along a given loading path can be determined by carrying out the appropriate homogenization calculations wherein the presence of the microballoons is resolved explicitly and the proper bifurcation and postbifurcation analyses are accounted for from the outset; see Geymonat et al. (1993) for the relevant theoretical framework and Michel et al. (2007;2010) for two pertinent examples.
Alternatively, the evolution equation can be prescribed phenomenologically on the basis of macroscopic observations. In this work, we follow this latter top-down approach.
While still relatively scarce and almost exclusively limited to uniaxial compression loading paths, all the currently available experiments on elastomeric syntactic foams consistently indicate that initially intact microballoons remain roughly so provided that the applied loads are relatively small. As the applied loads increase, there is a distinct point at which the microballoons start fracturing or buckling in significant numbers. Upon continuing loading, more and more microballoons may fracture or buckle until all of them had done so, that is, c p = c, and the process is exhausted.
To model the above-outlined behavior along any given loading path, say parameterized by the time t ∈ [0, T ], we propose the following basic class of rate-independent evolution equations: with initial condition for c p . In these expressions, we have made use of the standard "dot" notation to indicate derivatives with respect to the time t and f (I 1 , J) is any (suitably well-behaved) non-negative material function subject to the constraints c 0 p ≤ f (I 1 , J) ≤ c for all I 1 , J > 0. Remark 9. The initial condition (16). For newly fabricated elastomeric syntactic foams, the initial volume fraction of fractured/buckled microballoons is expected to be vanishingly small and thus c 0 p = 0+. In practice, it suffices to set c 0 p = 10 −5 . Remark 10. Irreversibility. The class of evolution equations (15) describes the processes of fracturing and buck-ling of the microballoons as irreversible. While the fracture of glass microballoons is indeed irreversible,  have shown that buckled polymer microballoons may unbuckle back to their original spherical form upon unloading. In this work, we restrict attention to evolution equations of the form (15) and will not consider generalizations that account for the possible reversible buckling of polymer microballoons.
Remark 11. The function f (I 1 , J). In an ideal world, the function f (I 1 , J) would be determined for a given elastomeric syntactic foam of interest by fitting its mechanical response for a variety of loading conditions that span the entire space of uniform deformations. At present, however, most of the experimental data available in the literature is narrowly restricted to uniaxial compression. So one must resort to the use of a model to extrapolate that limited available data to the entire deformation space.
In this work, for demonstration purposes, we shall consider functions of the form f (I 1 , J) = g(J). (17) It is expected from long-standing knowledge on other types of foams (syntactic or otherwise) that volume-decreasing loading conditions when J < 1 should favor the fracturing or buckling of the underlying microballoons in elastomeric syntactic foams, while volume-increasing loading conditions when J > 1 should hinder them; see, e.g., the experimental results of Gupta et al. (2010) and Landauer et al. (2019). The constitutive choice (17) is among the simplest that allow to describe such a "tension"/"compression" asymmetry.
As illustrated below in Section 4, fits to uniaxial compression data indicate that the particular form where w r ≥ 0 with R r=1 w r = 1, γ r > 0, β r > 0, and δ r > 0 are c-dependent material parameters, is descriptive of a variety of elastomeric syntactic foams.

The constitutive response
Having introduced the free-energy function (14) and the evolution equation (15) with (17), it is now a simple matter to spell out the constitutive relation (13) that they imply. Indeed, the first Piola-Kirchhoff stress S reads as where I 1 is given explicitly by (14) 2 and where c p is defined, in terms of the initial volume fraction of fractured/ buckled microballoons c 0 p and the deformation history, by the ordinary differential equatioṅ Numerical solution of the evolution equation (20). In general, because of its two distinct branches, the evolution equation (20) for c p needs to be solved numerically. All the same, the construction of such numerical solutions is trivial. Indeed, upon partitioning the time interval of interest [0, T ] into the discrete times 0 = t 0 , t 1 , ..., t m , t m+1 , ..., t M = T , the solution for c m+1 p = c p (t m+1 ) at the discrete time t m+1 is given explicitly in terms of the solution c m p = c p (t m ) at the previous discrete time t m and the associated values of the invariant J m+1 = J(t m+1 ) at t m+1 and J k = J(t k ) at all previous times t k by the rule (21) Clearly, for functions g(J) that are non-decreasing, such as (18), and loading paths along which J increases monotonically in time, the discrete solution (21) can be written continuously in time as c p (t) = g(J(t)).
The limit of small deformations. In the limit of small deformations as F → I, recalling the linearization properties Ψ m (3) = 0 and Ψ m (3) = µ m /2 of the stored-energy function Ψ m (I 1 ) characterizing the elastic response of the underlying elastomeric matrix, where, again, µ m denotes its initial shear modulus, the constitutive relation (19) reduces asymptotically to where = (F + F T − 2I)/2 stands for the infinitesimal strain tensor, and denote the shear and bulk moduli of the elastomeric syntactic foam, and where c p is, of course, still defined by the evolution equation (20).
Specialization to uniaxial compression. As mentioned above, the archetypical experiment to probe the mechanical behavior of elastomeric syntactic foams has been uniaxial compression. Customarily, a cylindrical specimen of the elastomeric syntactic foam of interest is stretched uniaxially along its longitudinal direction, say by λ < 1 in the e 3 direction, while its perpendicular directions e 1 and e 2 are kept traction free so that the deformation gradient and first Piola-Kirchhoff stress throughout the specimen take the spatially homogeneous forms Here, S and λ l are, respectively, the nominal axial stress and the lateral stretch induced in the specimen by the applied stretch λ. Now, according to the constitutive relation (19), the nominal axial stress S is given explicitly in terms of λ and λ l by while the lateral stretch λ l is defined implicitly in terms of λ as the smallest root larger than 1 of the nonlinear algebraic equation In these last two expressions, and c p is defined by the evolution equation (20), where now J = λλ 2 l in terms of the applied stretch and the induced lateral stretch.
In general, the nonlinear algebraic equation (25) does not admit an explicit solution for λ l and hence must be solved numerically. In the limit of small deformations as λ → 1, however, it is not difficult to show that the solution is given by where ν = 3κ − 2µ 2(3κ + µ) = 6 + c p 12 + 11c p denotes the Poisson's ratio of the elastomeric syntactic foam. In turn, it is not difficult to deduce that the stressstretch relation (24) reduces also explicitly to in the limit as λ → 1. The coefficient E is nothing more than the Young's modulus of the elastomeric syntactic foam.
The two material functions and two material parameters in (19)-(20) and their calibration from experiments. The proposed constitutive relation (19)-(20) contains two material functions, Ψ m (I 1 ) and g(J), and two material parameters, c and c 0 p . Physically, again, the function Ψ m (I 1 ) describes the elastic response of the elastomeric matrix in the foam. The function g(J) describes the evolution of the volume fraction c p of fractured/buckled microballoons along any given loading path. Moreover, c stands for the volume fraction of microballoons, while c 0 p denotes the initial volume fraction of fractured/buckled microballoons.
For the case of newly fabricated syntactic foams when the volume fraction c of microballoons is known from the fabrication process, one can presume that c 0 p = 0+, the function Ψ m (I 1 ) can be calibrated from a routine experiment such as uniaxial tension and/or compression on the matrix material, while the function g(J) can be calibrated from a uniaxial compression experiment on the foam. The last calibration process amounts to determining -for instance, by least squares -the value of c p in the stressstretch relation (24) that best matches the stress-stretch experimental data. Ideally, the uniaxial compression experiment should be carried out up to sufficiently large deformations that probe the entire evolution domain from c p = c 0 p to c p = c. For cases when no or only partial information is available about the fabrication process, one would have to make use of specific functional forms for Ψ m (I 1 ) and g(J), and then fit their material parameters together with c and c 0 p to the available experimental data on the foam. In this more elaborate calibration process, the formulas (22), (23), and (26) for the initial shear, bulk, and Young's moduli are likely to be useful in determining the volume fraction c of microballoons as well as the share c 0 p of these that are initially fractured/buckled alongside the initial shear modulus µ m of the underlying matrix material.

Comparisons with experiments
In this section, with the dual objective of illustrating its use and showcasing its descriptive and predictive capabilities, we confront the proposed model (19)-(20) with the experimental results of Brown et al. (2016; for silicone syntactic foams filled with glass microballoons. These pertain specifically to quasistatic uniaxial compression tests of syntactic foams with volume fractions of microballoons c = 0.46, 0.37, 0.30, 0.20, 0.10, wherein the matrix is the popular PDMS Sylgard 184 supplied by Dow Corning, with weight ratio 10:1 of base to curing agent, and the randomly distributed microballoons are the A16/500 glass bubbles supplied by 3M, which are made of sodalime-borosilicate glass and feature an average diameter of 60 microns and a thickness of 1 micron.  (3) for the PDMS elastomer used as the matrix material in the elastomeric syntactic foams. The associated initial shear modulus µm = µ 1 + µ 2 is also listed for convenience. Given that the volume fractions c = 0.46, 0.37, 0.30, 0.20, 0.10 of microballoons as well as the initial volume fraction c 0 p = 0+ of fractured microballoons in each of the fabricated foams are known from the outset, our first task is to calibrate the material function Ψ m (I 1 ). To do so, we choose Ψ m (I 1 ) to be given by the stored-energy function (3) and determine its four material parameters µ 1 , µ 2 , α 1 , α 2 by fitting the loading part of the stress-stretch data reported by Brown et al. (2016) for the uniaxial compression of an unfilled specimen; consistent with the presentation of the experimental results by Brown et al. (2016;, the stretch data is presented in the form of nominal compressive strain e = 1 − λ in all the results that follow. Table 1 lists the results obtained from that fit. Figure 6 compares the stress-stretch relation predicted by the stored-energy function (3), with the material parameters listed in Table 1, and the experimental data. Two observations are immediate. The first one is that the model is in good agreement with the experiment. The second one is that PDMS Sylgard 184 with 10:1 ratio of base to curing agent exhibits a small but noticeable hysteretic behavior. According to the recent investigations of Poulain et al. (2018) and Leonard et al. (2020) on the same type of elastomer, this behavior is not due to viscous dissipation, but instead to the Mullins effect. We do not model this effect here, but discuss it further in Section 5.
The initial linear elastic response of the foams. Having determined the elastic response of the PDMS matrix in the foams, we proceed by comparing the Young's modulus (26) predicted by the model with the experimental results at the initial time t = 0, when the initial volume fraction of the fractured microballoons is known from the fabrication process to be given by c p = c 0 p = 0+. Figure 7 shows that comparison. It is plain that the model predictions are in fairly good agreement with the experimental results for the entire range of volume fractions c of microballoons considered. As argued by Leonard et al. (2020) for the same type of PDMS elastomer filled with solid -as opposed to hollow -glass particles, the minor discrepancies between the rigorous homogenization result (26) and the experiments in Fig. 7 are likely due to the fact that the elasticity of PDMS Sylgard 184 is extremely sensitive to its cross-link density. Indeed, small differences in the base-to-curing-agent ratio and/or the curing time in the fabrication process of the PDMS Sylgard 184 specimens can lead to noticeable differences in their Young's modulus. The discrepancies may also be partly due to the clustering of microballoons.  Table 1.
The finite-deformation response of the foams with c = 0.37 and c = 0.20. We now turn to confront the model predictions with the experimental results at finite deformations for the foams with volume fractions c = 0.37 and c = 0.20 of microballoons, which unequivocally involve the evolution of the volume fraction c p of fractured microballoons.
To that end, we first need to calibrate the material function g(J) in the evolution equation (20) Table 2 for the elastomeric syntactic foams with microballoon volume fraction: (a) c = 0.37 and (b) c = 0.20. c = 0.37 w 1 = 0.6312 γ 1 = 15.00 β 1 = 0.4478 δ 1 = 5.00 w 2 = 0.3688 γ 2 = 80.55 β 2 = 1.00 δ 2 = 0.69 c = 0.20 w 1 = 0.4667 γ 1 = 14.99 β 1 = 0.4478 δ 1 = 5.00 w 2 = 0.5333 γ 2 = 24.91 β 2 = 0.6833 δ 2 = 5.00 the two foams. We do so by choosing g(J) to be given by the function (18) with R = 2 and by determining via least squares the values of its seven material parameters w 1 , γ 1 , β 1 , δ 1 , w 2 = 1 − w 1 , γ 2 , β 2 , δ 2 for which the stress-stretch relation (24) predicted by the model best matches the measured stress-stretch responses of the two foams in uniaxial compression. Table 2 lists the material parameters that result from such a fitting process. Precisely, the parameters for the foam with c = 0.37 were obtained by fitting the loading part of the single-cycle uniaxial compression test reported by Brown et al. (2016). The parameters for the foam with volume fraction c = 0.20 were obtained by fitting the loading parts of a ratcheted cyclic uniaxial compression test reported by Brown et al. (2018). Figure 8 presents plots of the resulting functions g(J) for the foam  Having fully calibrated both material functions Ψ m (I 1 ) and g(J), the model (19)-(20) for the elastomeric foams with c = 0.37 and c = 0.20 is now complete and can be readily deployed to determine the mechanical response of those foams under any loading path of choice. Figure 9(a) confronts the stress-stretch response (24) predicted by the model with the experimental data of Brown et al. (2016) for the foam with volume fraction c = 0.37 of microballoons that is loaded in uniaxial compression up to a compressive strain of 1 − λ = 0.5 and then unloaded. Figure 9(b) presents the corresponding evolution of the volume fraction c p of fractured microballoons predicted by the model as a function of the applied com-pressive strain 1 − λ.
It is plain that the stretch-stress model prediction is in fairly good agreement with the experimental data, especially for the loading part of the test. The larger discrepancies seen during the unloading may be due, at least in part, to the fact the model does not account for the Mullins effect in the underlying PDMS matrix, which, as noted in the discussion of Fig. 6 above, is noticeable. Another possible contributing factor is that fracturing of microballoons might also occur during the unloading part of the experiment, a process that the basic choice (18) for the material function g(J) does not account for. Indeed, Fig. 9(b) shows that the model predicts that c p remains essentially vanishingly small up to a compressive strain of about 1 − λ = 0.14, at which point c p starts to increase monotonically with increasing strain. At the maximum applied strain of 1 − λ = 0.5, the volume fraction of fractured microballoons reaches the value c p = 0.235. During the entire unloading part of the test, the model predicts that c p remains constant at c p = 0.235. Figure 10 shows analogous results to those presented in Fig. 9 for the foam with volume fraction c = 0.20 of microballoons. In this case, the results correspond to the tenth cycle in a ratcheted loading/unloading uniaxial compression test in which the maximum applied compressive strain 1 − λ in each successive loading cycle is increased by 0.05; see Fig. 7(a) in Brown et al. (2018). As was the case in Fig. 9, it is immediate from Fig. 10(a) that the stress-stretch response predicted by the model is in good agreement with the loading part of the test. It is equally clear that the model underpredicts the softening of the unloading part, significantly more so than for the foam with c = 0.37. Again, this is partly due to the Mullins effect of the matrix. As we elaborate in the next section, it may also be due to the fact that microballoons might continue to fracture during the unloading part of the test, which is a process, again, that the basic choice (18) for the material function g(J) does not account for.

Fracturing of microballoons during unloading
To gain definite insight into the respective extents to which the matrix Mullins effect and the fracturing of microballoons impact the unloading branches in the results presented in Figs. 9(a) and 10(a), one would have to first generalize the homogenization result (8) to account for the Mullins effect in the underlying matrix material and then generalize the free-energy function (14) for the elastomeric syntactic foams accordingly. This is an admittedly challenging endeavor that is beyond the scope of this work.
Nevertheless, one can readily gain partial insight by neglecting altogether the effect of the Mullins effect and considering that the only active dissipation mechanism during the unloading branches is the fracturing of microballoons. To do so, it suffices to determine the value of c p in the stress-stretch relation (24) predicted by the model that best matches the stress-stretch experimental data in Figs. 9(a) and 10(a) for both the loading and the unloading  Figure 11: (a) Comparison of the stress-stretch experimental data presented in Fig. 9(a) with the response (24) predicted by the model with function (3), the material parameters listed in Table 1, and the volume fraction cp of fractured microballoons presented in part (b). Note that cp continues to increase during the unloading.
branches. The results of such a fitting process are presented in Fig. 11 for the elastomeric foam with volume fraction c = 0.37 of microballoons and in Fig. 12 for the foam with c = 0.20. Both sets of results show that the model is in good agreement with the experimental data and therefore suggest that fracture of microballoonsrather remarkably -may indeed occur during the unloading of the foams.

Final comments
Summing up, the results presented in Figs. 7 through Fig. 12 point to the capabilities of the proposed framework/model to accurately explain, describe, and predict the mechanical behavior of elastomeric syntactic foams along arbitrary quasistatic loading paths from the bottom up, in particular, directly in terms of the elastic response of the underlying elastomeric matrix, the initial volume  (24) predicted by the model with function (3), the material parameters listed in Table 1, and the volume fraction cp of fractured microballoons presented in part (b). Note that cp continues to increase during the unloading.
fraction of microballoons, and their fracturing or buckling. The results also point to the critical need to carry out more experiments that go beyond uniaxial compression and, in particular, that explore the possibility of the fracturing of microballoons during "unloading" in order to guide the construction of accurate evolution equations for the volume fraction c p of fractured/buckled microballoons. Plans to carry out such experiments are underway. We close by noting that the mathematical richness and simplicity of the proposed free-energy function (14) for elastomeric syntactic foams may readily be transcribed to model other types of elastomeric foams.