# On the identification of nonlinear behavior of highly filled elastomers

Dimitri Jalocha, Aurélie Azoug, Andrei Constantinescu, Robert Nevière

► To cite this version:

Dimitri Jalocha, Aurélie Azoug, Andrei Constantinescu, Robert Nevière. On the identification of nonlinear behavior of highly filled elastomers. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717834

HAL Id: hal-01717834

https://hal.science/hal-01717834

Submitted on 26 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# On the identification of nonlinear behavior of highly filled elastomers

D.Jalocha<sup>a,b</sup>, A.Azoug<sup>b</sup>, A. Constantinescu<sup>a</sup>, R.Nevière<sup>b</sup>

**Résumé** — Solid propellants are highly filled elastomers used as propulsion medium. Because of the high volume filler fraction of 80% – 90%, several sizes of fillers are used and as a consequence of the complex microstructure the material presents important nonlinearities in the macroscopic behavior. An abrupt increase in the viscoelastic storage and loss modulus has been experimentally observed in [2]. The objective of this paper is to explore numerically a highly filled elastic composite under the assumptions of "rigid" fillers, an elastic matrix and large strains. The present numerical analysis is based on several artificially created microstructures. The results show the influence of the material behavior of the matrix on the global stiffness. The computed mechanical response is in agreement with experiment observations. The mechanical response of the composite under a uniaxial loading presents a linear and a nonlinear part. Both are analyzed to understand the microstructural effects that lead to the macroscopic behavior. The apparent modulus of the linear part increases with respect to the volume fraction of fillers and the stiffening effect occurring in the nonlinear part is not dependent of the matrix behavior.

**Mots clés** — Propellant - Nonlinearity - Uniaxial tensile simulation - Microscopic/Macroscopic

# 1 Introduction

Composite materials are usual for various applications which exploit their weight or mechanical and energetic properties. Within the class of highly filled elastomers, the propellants are used for solid propulsion and they provide the combustion material for the oxidation-reduction reaction and the structure of the motor. Therefore, they should withstand both energetic and mechanical design criteria. Fillers represent 80% of the total volume and play the role of the reducer and the oxidizer. They are included in an HTPB viscoelastic matrix, along with plasticizers and bonding agents [2].

Filled and unfilled elastomers have been the object of several recent papers. Within the domain of hyperelastic and viscoelastic material behavior, models have been proposed, starting from both a phenomenological basis ([9] and [8]) and a homogenization theory ([11]). A special interest for the present study is the model of Lion and Retka which takes into account the prestress ([6]).

In the experimental domain, Azoug has explored the behavior of different propellants with and without prestrain (see [10]). When results are analyzed in terms of storage modulus ($E'$) and loss modulus ($E''$) a nonlinear evolution is observed (see Fig.1). The evolution of $E'$ and $E''$ has been modeled with a phenomenological model but no relation has so far been given between the mechanical behavior of the microstructure and the different loading parameters. The objective of this paper is to link the microscopic and the macroscopic effects. The method proposed is based on numerical simulations of artificial microstructures and limited to the elastic behavior under a large strain assumption. A DMA is a rate dependent experimentation. To simplify the numerical resolution, assumptions are made. The storage modulus (in dynamic) is linked with the apparent modulus (in quasi-static), so an increase of the storage modulus implies an increase of the apparent modulus (stiffening of the material). So hyperelastic simulation are performed and only the apparent modulus is analyzed now, the loss modulus is forgotten. The term "modulus" denotes now the apparent modulus (slope at a point of the Cauchy stress strain curve). The objective is to understand the relations between the microscopic and macroscopic strains and the role of the constitutive law of the matrix in the apparition of the nonlinear macroscopic elastic moduli. We

<footnote><sup>a</sup> LMS/CNRS UMR 7649, Ecole Polytechnique, France
constant@lms.polytechnique.fr; jalocha@lms.polytechnique.fr</footnote>

<footnote><sup>b</sup> Herakles, Centre de Recherches du Bouchet, 91710 Vert-le-Petit, France</footnote>
shall analyse a series of tensile and shear experiments on different microstructures and on different filler volume fractions (denoting by *f* in this paper). Two particular polynomial constitutive laws are explored for the matrix : Neo Hookean (NH) and Arruda-Boyce (AB) model.

<figure><img src="image_1.png" /><figcaption>Fig. 1 – The evolution of the storage and loss moduli $E'$ and $E''$, with respect to [10]</figcaption></figure>

## 2 Approach

The propellant is a HTPB polymeric matrix filled with Aluminum and Ammonium particles. The numerical microstructures are random polydisperse distributions of spheres into a cube, as described in [7]. A volume fraction of 80% can't be obtained numerically but a maximum of 70% is possible. In order to keep the element size homogeneous in the numerical microstrucuture, the ratio between the largest and the lowest diameter of the fillers is fixed to 10. Fig.2 shows the real and a simulated microstructure, three types of microstructures are created and their parameters are defined in the Table 1. Only the radius of the fillers are chosen, the spacial distribution is random.

<figure><img src="image_2.png" /><figcaption>Fig. 2 – Real spacial distribution (a) and simulated spacial distribution with periodic boundary conditions (b), corresponding to an enlarged zone</figcaption></figure>

<table><thead><tr><td>$f$ (%)</td><td>Radius (mm)</td><td>Cube face length (mm)</td><td>Nodes</td><td>Limit strain</td></tr></thead><tbody><tr><td>30</td><td>0.1</td><td>0.5</td><td>7097</td><td>40%</td></tr><tr><td>50</td><td>0.1 - 0.05 - 0.02</td><td>0.3</td><td>7005</td><td>26%</td></tr><tr><td>70</td><td>0.1 - 0.05 - 0.02 - 0.015 - 0.01</td><td>0.3</td><td>33638</td><td>13%</td></tr></tbody></table>

Tableau 1 – Size parameters of the microstructures, number of nodes, and possible computation limit strain

The microstructure presents periodic boundary conditions given by Equation (1) is the mathematical translation of periodic boundary conditions (where $U^r$ denotes the displacement of the right face of the cube, $U^l$ the displacement of the left face, $\underline{F}$ the deformation gradient, $\underline{I}$ the identity matrix and $\delta X$ the difference of coordinates in the reference configuration between two points face to face). Therefore, the macroscopic strain is imposed by defining the displacement of only four corners.

$$U^r - U^l = (\underline{F} - \underline{I}) \delta X \qquad (1)$$
Compared to the matrix, the filler behavior is assumed to be rigid. A linear elastic behavior is chosen for this component, with a high value of Young modulus ($E_f = 6.6 \times 10^{10}$ Pa) and a Poisson ratio of $\nu = 0.159$ according to [5]. The mechanical behavior of the matrix is complex due to the presence of under-reticulation and plasticizer agents in the matrix, well developed in [2]. As a consequence the particular behavior is obtained only in the presence of the fillers. As the extraction of a microscopic specimen is not possible, we have experimentally tested a HTPB sample which presents some similitude with the binder, to have an order of magnitude of what the inside propellant matrix hyperelastic behavior could be, shown in Fig.3. For the mechanical behavior of the matrix we propose a polynomial model, described by the energy function given in Equation (2).

$$W = \sum C_{\alpha}(\bar{I}_{1}^{\alpha} - 3^{\alpha}) \qquad (2)$$

$$\bar{I}_{1} = \bar{\lambda}_{1}^{2} + \bar{\lambda}_{2}^{2} + \bar{\lambda}_{3}^{2} \qquad (3)$$

$$\bar{\lambda}_{i} = J^{-\frac{1}{3}} \lambda_{i} \qquad (4)$$

<figure><img src="image_3.png" /><figcaption>Fig. 3 – Experimental measures during an uniaxial tensile test of HTPB polymer at natural state</figcaption></figure>

$\bar{I}_1$ is the first invariant of the deviatoric part of the deformation gradient, defined by Equation (3), $\lambda_i$ the principal strain in the $i$ direction (see Equation (4)) and $J = \det(F)$. Next, we shall to focus this study on two particular polynomial models: the Arruda-Boyce (AB) model, defined in [1]

$$W = C_1 \left( \frac{1}{2}(\bar{I}_1 - 3) + \frac{1}{20\lambda_m^2}(\bar{I}_1^2 - 9) + \frac{11}{1050\lambda_m^4}(\bar{I}_1^3 - 27) + \frac{19}{7000\lambda_m^6}(\bar{I}_1^4 - 81) + \frac{519}{673750\lambda_m^8}(\bar{I}_1^5 - 243) \right) \qquad (5)$$

and the Neo-Hookean (NH) model

$$W = C_1(\bar{I}_1 - 3) \qquad (6)$$

It is interesting to compare the composite response for two matrix behaviors: in one hand a stiffening behavior (AB model) and in the other hand a non-stiffening behavior (NH model). We recall the propellant is assumed incompressible ($\det(F) = 1$). Thus during a uniaxial tensile test (supposed in the direction 1), a unique strain parameter $\lambda$ drives the behavior:

$$\lambda = \lambda_1 ; \lambda_2 = \lambda_3 = \frac{1}{\sqrt{\lambda}} \qquad (7)$$

the invariants are:

$$\bar{I}_1 = \lambda^2 + \frac{2}{\lambda} \qquad (8)$$

$$\bar{I}_2 = \bar{I}_3 = 1 \qquad (9)$$

the relation between the stretch $\lambda$ and the Cauchy stress $\sigma$ is given by

$$\sigma = \lambda \frac{\partial W(\lambda)}{\partial \lambda} \qquad (10)$$

The computations have been performed using the ABAQUS® finite element code and the standard implementation of the constitutive laws, see Fig.4. The numbers of nodes are recapped in the Table 1.
Large strains are imposed, so nonlinear geometry is used. For each microstructure, a limit strain value is reached. At higher strain levels, the computation does not converge anymore, elements are too warp (see limit values in Table 1).

The finite strain assumption permits a rearrangement of the microstructure during loading. With increasing load, one can observe relative movements and rotations of the particles in the polymeric matrix, which conduct to the creation of a high strain chains and high strain bands. A similar structure was observed during a two dimensional analysis as reported in [2]. The ratio between the macroscopic strain and the average of the microscopic strains in the chain is around 4.5 (see next part).

<figure><img src="image_4.png" /><figcaption>Fig. 4 – Strain in the load direction for different steps during a uniaxial tensile simulation on 70% filled composite, macroscopic strain of (a) 2%, (b) 7% and (c) 13%</figcaption></figure>

## 3 Results and discussion

In a preliminary stage we verify the isotropy of the created microstructures, in order to use general results on isotropic materials. One microstructure with 70% of filler is submitted to uniaxial tensile tests in the three directions and three shear tests. Results (in Fig.5) exhibit a small scattering of the mechanical response according to the direction. For an imposed stretch of 1.13 the average of the Piola stress is $3.95^{\pm 0.26} \times 10^8$ Pa. For an imposed shear of 0.03, the average of the Piola stress is $5.4^{\pm 0.22} \times 10^6$ Pa. As a consequence microstructures will be considered as isotropic and $\pm 5\%$ of error in the tangent stiffness can be considered as normal.

<figure><img src="image_5.png" /><figcaption>Fig. 5 – (a) Uniaxial tensile simulations in the direction Ox, Oy and Oz and (b) shear simulations in the direction OyOz, OxOz and OxOy for the 70% filled microstructure to evaluate the isotropy</figcaption></figure>

For all uniaxial tensile simulations of each microstrucutre, the strain ($\epsilon = \frac{\delta l}{l}$) and the Cauchy stress ($\sigma = \frac{F}{a}$, where a is the real surface) are computed and the results are shown in Fig.6(a). The Arruda-Boyce model is used for the first part here. The errors due to the random distribution are displayed in the Table 3. By the centered difference method, the slope at each point of the stress/strain curve is computed, to obtained the modulus, plotted in Fig.6(b).

The stiffening effect of the composite when fillers are added in the matrix is a well known process, as described in [3]. More fillers are added, more the composite is stiffer, and this effect evolve nonlinearly with respect to the volume fraction of fillers. Results of this simulation demonstrate it one more time.
<figure><img src="image_6.png" /><figcaption>Fig. 6 – (a) Results of the uniaxial tensile simulation for the 30%, 50% and 70% filled composite (b) Apparent modulus for each microstructure, computed by the centered difference method from the results of uniaxial tensile simulation. The Arruda-Boyce model for the matrix is used.</figcaption></figure>

<table><caption>Tableau 2 – Errors due to the random distribution for the uniaxial tensile simulation of each microstructure</caption><thead><tr><td>f (%)</td><td>Strain</td><td>Cauchy stress (Pa)</td></tr></thead><tbody><tr><td>30</td><td>30%</td><td>6,7<sup>±1</sup> × 10<sup>7</sup></td></tr><tr><td>50</td><td>23%</td><td>1.48<sup>±0.07</sup> × 10<sup>8</sup></td></tr><tr><td>70</td><td>11%</td><td>3.33<sup>±0.27</sup> × 10<sup>8</sup></td></tr></tbody></table>

The comparison of the nonlinearities of Fig.6(b) with experimental results of Fig.1 presents several similarities of shape. A linear part, occurring for small strains (ε < 0.01, see Fig6(b)), and a nonlinear part where an elbow is visible.

The first analysis corresponds to the apparent elastic moduli at small strain. Denoting by $\bar{E}$ and $E$ the modulus of the composite and the matrix respectively, it has been shown in [7] that, if the matrix has an energy function depending only on the first invariant $\bar{I}_1$:

$$ \frac{\bar{E}}{E} = \frac{1}{(1-f)^{5/2}} \qquad (11) $$

This analytical result has been verified for FEM simulations with volume fractions of fillers up to $f = 30\%$ in [7]. Fig.7(a) shows the evolution of $\bar{E}$ up to a volume fraction of $f = 70\%$. The modulus increases nonlinearly with the filler volume fraction. The analytical function tends to be still valid up to 70%. Microstructures are considered isotropic, so a comparison with the Voigt, Reuss and Hashin and Shtrikman (HS+ and HS-) bounds can be made (see [4]). Results presented in Fig.7(b) tend to confirm the validity of the obtained modulus in the linear elastic phase of the composite behaviors.

<figure><img src="image_7.png" /><figcaption>Fig. 7 – (a) Comparison between the analytical model given by [7] and highly-filled propellant FEM data (b) Comparison between analytical model and Voigt, Reuss, and HS- and HS+ bounds. AB model is used for the matrix behavior</figcaption></figure>
The next question concerns the apparition of the nonlinearity characterized by the "elbow" in the representation of the apparent modulus versus the strain (see Fig.6(b)). Under a certain strain value, the apparent modulus is constant. But over this limit, the material becomes stiffer, and the modulus increases strongly. We shall explore the influence of the matrix behavior on the onset of the nonlinearity. Two mechanical behaviors are focused on, a Neo-Hookean model and an Arruda-Boyce model.

First, let us consider that the matrix follows an Arruda-Boyce model (see Equation 5). The supposed inside propellant matrix behavior is linear up to a very large strain (for $\epsilon = 75\%$). After this phase an elbow occurs. That implies a constant modulus until the end of the linear phase, and an increase of the modulus value at higher strains. In unfilled elastomers, this limit appears for large strain (for $\epsilon = 75\%$). In the case of filled composites, macroscopic behavior follows the same outline and a strong stiffening takes place. The strain at which the nonlinearity appears decreases when the filler fraction increases (see Fig.6(a)). The shape of the curve for filled composites is the same as the one for the unfilled matrix, but it seems to undergo an homothety along the strain axis.

Now, suppose the inside propellant matrix can be modeled by a Neo-Hookean model. After a certain limit value, the "elbow" appears again. The stiffening effect is independent of the choice of the matrix behavior. This is a microstructural effect. To compare between the two models, computation for the 30% filled microstructure is shown in. Fig.8. When the matrix behavior is naturally not stiffening (typically a Neo-Hookean model), the composite behavior is stiffening. When the matrix behavior is stiffening (for example an Arruda-Boyce model) the composite behavior is also stiffening. The difference appears only to the strain value of the stiffening effect. For the Neo-Hookean model, the stiffening occurs at a larger strain than for the Arruda-Boyce model. Whatever the matrix behavior, the composite behavior will present an increase of the modulus value with respect to the strain. Calculations are made with a stiffening polynomial energy function (Arruda-Boyce) for the matrix in order to observe this effect for reasonable strain value. A different choice of the matrix behavior affects the composite mechanical behavior only by a shift of the stress/strain curve along the strain axis.

<figure><img src="image_8.png" /><figcaption>Fig. 8 – Results of the computation for the 30% filled composite. (a) represents the response of the matrix only, for an uniaxial tensile simulation, modeled with a NH model and an AB model. (b) shows the behavior of the 30% filled composite with a NH model and an AB model for the matrix behavior</figcaption></figure>

Another question concerns the strain distribution in the matrix for a given macroscopically applied strain. It is well known that, due to the high contrast in elastic moduli between the filler and the binder, only the binder will support the applied strain. The imposed macroscopic strain or stretch ($\epsilon_M$ or $\lambda_M$) is not the microscopic strain or stretch ($\epsilon_m$ or $\lambda_m$) of each component of the composite. Fig.9 presents the distribution of the strain inside the matrix. The volume of each discrete finite element of the matrix and the associated stretch in the Gauss point of this element are recorded. A volume distribution of the microscopic strain can be built. The same recording is made for all the microstructures, with an imposed macroscopic stretch of 1.13. The distributions seem to follow a Gaussian Law and the average of the microscopic stretch increases with the filler volume fraction.

Defining the ratio $\frac{\epsilon_m}{\lambda_M}$, a comparison between macroscopic (imposed) strain ($\epsilon_M$) and the average of the microscopic strain ($\epsilon_m$) can be made with respect to the percentage of filler. The evolution of this ratio is plotted in Fig.10. A fit is made to find an empirical model given the evolution of this ratio in function of the volume fraction of filler $f$, see Equation (12). The model is empiric, but presents the same form as the previous analytical model Equation (11).
<figure><img src="image_9.png" /></figure>

(a)

<figure><img src="image_10.png" /></figure>

(b)

<figure><img src="image_11.png" /></figure>

(c)

<figure><img src="image_12.png" /><figcaption>Fig. 10 – $\langle \epsilon_m \rangle / \epsilon_M$ with respect to the volume fraction of filler and its model</figcaption></figure>

$$ \frac{\langle \epsilon_m \rangle}{\epsilon_M} = \frac{1}{\left(1 - \frac{3f}{4}\right)^2} \qquad (12) $$
# 4 Conclusion

The principal performances of propellants are energetic, but the mechanical behavior must be known to design an efficient motor. Due to the high volume ratio of filler, some nonlinearities appear, more especially dynamical nonlinearities in function of the prestrain of the propellant sample. Although it is possible to measure those nonlinearities with DMA tests, explaining why they occur is not trivial. One approach consists in using numerical simulations. Microstructures are generated, and following some simplifying assumptions (considering an hyperelastic problem, choosing the propellant matrix behavior), uniaxial tensile simulations are performed. Modeling the propellant matrix behavior with either a Neo-Hookean or an Arruda-Boyce model, the resulting strain/stress curve shows two parts. In the linear part, the apparent modulus is constant with respect to strain and increases according to the filler volume fraction. In the nonlinear part, the apparent modulus exhibits an "elbow" and increases with respect to strain. Since the stiffening effect does not qualitatively depend on the chosen matrix constitutive law, the results prove that the nonlinear macroscopic behavior originates from the heterogeneity of the microstructure. Finally, the microscopic matrix strain is not equal to the imposed macroscopic strain. Since the fillers are rigid, the binder undergoes the entire displacement. Numerically, the average of the microscopic stretch can be 1.6 times the macroscopic stretch of 1.13.

## Acknowledgments :

This work has been partially financed by the DGA. The authors would like to thank $M^{me}$ Amiet (DGA) for supporting the project and K. Danas for fruitful discussions.

## Références

[1] E Arruda and M Boyce. A three-dimensional model for the large stretch behavior of rubber elastic materials. *Journal of Mechanics of Physics and Solids*, 41(2):389-412, 1993.

[2] Aurelie Azoug. *Micromecanismes et comportement macroscopique d'un elastomere fortement charge*. PhD thesis, Ecole Polytechnique, 2010.

[3] Yoshihide Fukahori. Generalized concept of the reinforcement of elastomers. *Rubber Chemistry and Technology*, 80:701-725, 2007.

[4] Z Hashin and S Shtrikman. A variational approach to the elastic behavior of multiphase minerals. *Journal of Mechanics of Physics and Solids*, 11:127-140, 1963.

[5] Jinn-Shing Lee and Chung-King Hsu. Thermal properties and shelf life of hmx htpb based plastic bonded explosives. *Thermochimica Acta*, 393:153-156, 2002.

[6] A Lion, J Retka, and M Rendek. On the calculation of predeformation-dependent dynamic modulus tensors in finite nonlinear viscoelasticity. *Mechanics Research Communications*, 36:500-520, 2009.

[7] Oscar Lopez-Pamies, Taha Goudarzi, and Kostas Danas. The nonlinear elastic response of suspensions of rigid inclusions in rubber :ii a simple explicit approximation for finite-concentration suspensions. *Journal of the Mechanics and Physics of Solids*, 61:19-37, 2013.

[8] S Ozupek and E Becker. Constitutive modeling of high-elongation solid propellants. *Journal of Engineering Materials and Technology*, 114:111-115, 1992.

[9] Paul Steinmann, Mokarram Hossain, and Gunnar Possart. Hyperelastic models for rubber-like materials : consistent tangent operators and suitability for treloar's data. *Archive of Applied Mechanics*, 82:1183-1217, 2012.

[10] Anders Thorin, Aurelie Azoug, and Andrei Constantinescu. Influence of prestrain on mechanical properties of highly-filled elastomers : Measurements and modeling. *Polymer Testing*, 31:978-986, 2012.

[11] F Xu, N Aravas, and P Sofronis. Constitutive modeling of solid propellant materials with evolving micro-structural damage. *Journal of the Mechanics and Physics of Solids*, 56:2050-2073, 2008.