Prestrain-dependent viscosity of a highly filled elastomer: experiments and modeling

Highly filled elastomers exhibit a complex microstructure made up of rigid fillers bounded by a thin layer polymeric matrix. The interactions between the fillers and the binder amplify locally the applied strains and induce a nonlinear viscoelastic behavior. The aim here is to analyze the influence of prestrain on the viscoelastic behavior. This paper proposes a prestrain-dependent viscoelastic constitutive model. The model is a superposition of three relaxation spectra, each corresponding to a family of polymer chains, and can be regarded in either its continuous or discrete expression. More specifically, one of these relaxation spectra is modified to assure the prestrain sensitivity. The parameters of the discrete model are identified from relaxation and DMA experiments performed on a solid propellant, and the obtained predictions match closely the experiments. The novelty of the analysis proposed in this paper is threefold. On the one hand, we report a new series of experimental measures, performed for a large range of frequencies for the DMA experiment and relaxation times for the relaxation experiment, and, on the other hand, we propose a constitutive law compatible with the principles of thermodynamics, which predicts closely the measurements. Finally, the analysis is performed comparing both relaxation and DMA experiments using the spectrum of relaxation times. A peculiarity of the present discussion is the novel identification method used, which identifies directly the relaxation times. This technique leads to models with smaller and optimum numbers of parameters than classical methods based on a logarithmic distribution of relaxation times.


Introduction
Composite materials are employed in various applications, in which they exploit their exceptional properties. Carbon fiber composites appreciated for the excellent stiffness-to-weight ratio are increasingly applied in the aerospace industry (Botelho et al. 2003). Reinforced rubbers are solicited in large elastic strain regime with damping and are recommended for the manufacturing of tires, sealing, etc. (Roeder and Stanton 1983). Solid propellants, another class of reinforced rubbers, are retained for their energetic properties in solid boosters (Ozupek 1989). The mechanical behavior of these different composites is generally viscoelastic and tailored with respect to the final application. Therefore, the modeling of the material behavior at different scales took an important place in the last decades (see, e.g., Huber and Tsakmakis 2000;Christensen 1980;Swanson and Christensent 1983).
The models describing the mechanical viscoelastic behavior of reinforced elastomers or rubbers are of different nature. Some studies propose hyperelastic constitutive models for filled elastomer (Arruda and Boyce 1993;Oscar 2010), and others reproduce the timedependent behavior of the material based on phenomenological consideration (Ozupek and Becker 1992;Simo 1987;Steinmann et al. 2012) or on homogenization theory (Xu et al. 2008). Specifically, nonlinear viscoelastic models have recently been explored by combining the physics at the scale of the microstructure with the mechanical observation at the macroscopic scale of the applications (Azoug et al. 2014a(Azoug et al. , 2014b(Azoug et al. , 2014c. Other studies (Lion and Kardelky 2004) proposed complex constitutive models in the large-strain regime taking into account the nonlinear effects due to the presence of high-volume fraction of fillers.
Reinforced polymers and rubbers will cover a large range of the volumetric filler fraction, lying between 20 % and 70 %. Small filler fractions are used in rubbers, whereas high filler fractions are found in propellants. We remark that generally the nonlinearity of the mechanical behavior increases with increasing filler fraction.
The aim of this paper is to explore the particular nonlinearity induced by a prestrain on the viscoelastic properties of a highly filled elastomer. This effect was observed during an experimental campaign on solid propellants, performed by Azoug (2010) and reported in several articles. In Thorin et al. (2012) the influence of prestrain on viscoelastic properties for different HTPB composite propellants was measured using dynamic mechanical analysis (DMA). The prestrain-dependent behavior was modeled using a modified generalized Maxwell model, where the stiffnesses of each Maxwell element was linearly dependent on prestrain. The paper (Azoug et al. 2013b) was devoted to a biaxial experiment having a prestrained DMA on one axis and a fixed prestrain on the other. The experiments showed that the storage and loss moduli increase with the prestrain under both loadings and that the nonlinear behavior is quantitatively modified by adding an orthogonal prestrain. Moreover, the modification of the behavior under a horizontal prestrain was canceled out by an increase of the vertical prestrain, a phenomenon that may be explained by the alignment of the fillers with the prestrain. The experimental results of a complete biaxial prestrained DMA experiment have been recently reported in Jalocha et al. (2015b). The experiments on propellant specimen confirmed the previous results. Moreover, it has been shown that the dependence of prestrain in biaxial experiments can be characterized using the second invariant of the prestrain. The limitation of the reported analysis is the unique frequencies where the viscoelastic nonlinearity has been measured. The novelty of the analysis proposed in this paper is threefold. On the one hand, we report a new series of experimental measures, performed for a large range of frequencies for the DMA experiment and large range of times for the relaxation experiment. On the other hand, we propose a constitutive law, compatible with the principles of thermodynamics, which predicts closely the measurements. Finally, the analysis is performed comparing both relaxation and DMA experiments.
The analysis of the data starts from the idea that relaxation spectrum and generalized Maxwell models represent continuous and approximate discrete viscoelastic models, respectively. As a consequence, the variations of the experimental observations as functions of the prestrain are considered as variations of the measured complex and relaxation moduli. In this paper, we propose a prestrain-dependent viscoelastic constitutive model. The model is a superposition of three relaxation spectra, each corresponding to a family of polymer chains, and can be regarded in either its continuous or discrete expression. More specifically, one of these relaxation spectra is modified to assure the prestrain sensitivity. The parameters of the discrete model are identified from relaxation and DMA experiments performed on a solid propellant, and the obtained predictions match closely the experimental. A particularity of the present discussion is the novel identification method used, which identifies directly the relaxation times. This technique leads to models with smaller and optimum numbers of parameters than classical methods based on a logarithmic distribution of relaxation times. A detailed mathematical description of this identification algorithm and a discussion of its robustness are reported in Jalocha et al. (2015a).
The paper is organized as follows. The first section presents the material and experiments. The second and third sections discuss the continuous and discrete models, respectively. The final sections present the results and conclusions.

Material
The material studied here is a solid propellant, used for the propulsion of rockets and classified as a highly filled elastomer. Fillers represent 70-80 % of the total volume and are bonded by a viscoelastic matrix (Azoug 2010), see Fig. 1(a). The performance of the rocket engine will depend on both energetic capacity of the filler and the mechanical properties of the propellant in order to guarantee the performance and integrity of the engine. The mechanism of the engine can be summarized as follows. The oxidation-reduction reaction decomposes the fillers and produces highly pressurized gases in the combustion chamber, and the gases are ejected through the nozzle of the engine and propel the rocket. Integrity of the engine imposes that the evolution of the propellant during the burning process stays stable and in accordance with the initial design in spite of the important mechanical loadings of the propulsion.
The high-volume fraction of fillers demands exceptional mechanical properties from the matrix that is a network of thin bridges surrounding the rigid fillers. The characteristic of the polymer matrix assured by its microscopic chemical and physical properties. The manufacturing process will generally create an over cross linked elastomer matrix. At the microscopic scale, one will encounter a multiphase binder as already discussed in Azoug (2010) and composed of: (i) a principal cross linked network polymer that bonds the solid particles together and (ii) nonlinked polymer chains floating freely through the principal network. An idealized structure is displayed in a schematic diagram in Fig. 1(b).

Experimental procedure
The viscoelastic behavior of the propellant at the macroscopic scale of the structure can be characterized using standard relaxation or DMA experiments. Let us next introduce some The viscoelastic behavior of a material is completely represented by the continuous spectrum of relaxation time H (τ ) as exposed, for example, in Findley et al. (1976). The spectrum H (τ ) is related to the relaxation modulus E(t) over the time domain and to the dynamical modulus E * (ω) = E (ω) + iE (ω) over the frequency domain by the following equations: According to Markovitz (1977), the relaxation and complex moduli permit to establish the classical viscoelastic constitutive equations, expressed in the time domain as and in the frequency domain as where σ and ε denote the stress and strain fields, respectively. The representations in (1) and (3) are equivalent, and the passage between the time and frequency domains is obtained through a direct or inverse Fourier transform. As a consequence, σ * is the direct Fourier transform of the stress σ , ε * is the direct Fourier transform of the strain ε, and the complex modulus E * is the direct Fourier transform of the relaxation modulus E. Let us further recall that the spectrum of relaxation times H (τ ) cannot be measured directly. It is only recovered from either the relaxation modulus E(t) or the complex modulus E * (ω), which are the direct outputs of the relaxation and DMA experiments; for a detailed explanation, see Knauss et al. (2006). As a consequence, we can conclude that the relaxation and the DMA experiments are dual to each other, one representing the time domain, and the other the frequency domain.
The time or frequency domain, covered by the relaxation or the complex modulus, respectively, tends to cover several orders of magnitude. Hence, it is difficult to cover experimentally the complete viscoelastic domain, and based on the time-temperature superposition principle, the so-called master curves of relaxation and complex modulus are constructed. The method is explained, for example, in Brinson (2008), Williams et al. (1955).
The objective here is to propose a viscoelastic constitutive model with prestrain sensitivity capable of predicting data from relaxation and DMA experiments. The main features of the experimental procedure are summarized in Fig. 2. In both cases, strains are imposed, stresses are measured, and the moduli are computed from the measured data. The data will be first collected in the format of a master curve encompassing several orders of magnitudes in the time or frequency range and used after for identification of the model. The prestrain is  Fig. 3. The mathematical duality of the measured data in the sense of the Fourier transform discussed implies that the effect of the applied prestrain should be represented by the same phenomena to cover both loading cases.
In the time domain, the relaxation modulus E(t) is measured in the following way: a tensile prestrain ε s is imposed, and then, after a short delay of 10 min considered large enough for the stress to reach a stabilized state, the relaxation test is started. The history of the imposed strain is displayed in Fig. 3(a). The relaxation test is performed imposing a tensile strain step superimposed with the prestrain: is the Heaviside function), and the decreasing stress σ (t) is recording in function of the time. The relaxation modulus is computed from the data as the secant modulus defined by E(t) = σ (t) ε 0 , where σ (t) is the difference between the relaxed stress and the residual stress as a function of the imposed prestrain.
In the frequency domain, the complex modulus E * (ω) is measured in the following steps: a tensile prestrain ε s is imposed, and, after a short delay of 10 min, the DMA test is started. The history of the imposed strain is displayed in Fig. 3(b). The DMA test is performed imposing sinusoidal strain around the given prestrain, ε(t) = ε s + ε d sin(ωt), and the stress σ (t) is recorded. The complex modulus is computed as the secant modulus using the formula where σ d is the amplitude of the sinusoidal stress, ε d is the amplitude of the sinusoidal strain, and δ is the time delay between the strain and stress signals measured during the stabilized period.
As mentioned before, both moduli are computed as secant moduli. This guarantees the consistency and correlation in the data processing between the two experiments.
In experiments involving viscoelastic filled materials under relaxation or DMA tests, one encounters the Payne effect (Payne 1962) or, close to the large-strain region, the Mullins effect (Mullins 1969). In the time domain, the Mullins effect characterizes an instantaneous softening of the stress-strain curve that occurs whenever the load increases beyond its prior all-time maximum value. For the propellant under scrutiny, we remark that a small softening of the Mullins type is generally observed. However, the effect is completely erased if a sufficiently long time lapse occurs between the two loadings. In the frequency domain, the Payne effect characterizes a softening of the norm of the complex modulus, which increases the applied dynamic strain amplitude. The effect is present in the material and was the reason to chose only small dynamic strain amplitudes. However, the chosen amplitude is smaller than that encountered in practical applications. Let us further remark that the two effects are generally considered to be physically closely related and could be at the origin of the small difference in viscoelastic properties measured from the relaxation and the DMA experiments. Except this fact, the process history of the material does not change the experimental results.
Mathematically, it is straightforward that the essential features of the model are expressed as the spectrum of relaxation times H (τ ), and the aim is to identify the spectrum from both relaxation and DMA experiments. Therefore, DMA results are transposed into the time domain using the inverse Fourier transform based on the empirical relation E(t) = E * ( t 2π ) .

Experimental results
The relaxation tests have been performed for the following prestrain values: ε s = 0; 1.2; 2.5; 3.7; 5; 6; 7.2; 8.4 % at temperatures between +20°C and +40°C. The choice of the maximum prestrain value of 8.4 % ensures that the experiments are performed in the incompressible domain of the material and that damage (i.e., filler debonding, matrix fracture, etc.) does not occur. The temperature range is coherent with the application constraints of the model. Moreover, results given at smaller temperature than the lower limit of +20°C can later be replaced with DMA experiments, which provide data of better accuracy. The upper limit +40°C covers, after the time-temperature superimposition, the application of the model. The master curves obtained at a reference temperature of T ref = 20°C for different prestrain values are presented in Fig. 4. Panel (a) presents the time evolution of the relaxation modulus E(t), and panel (b) displays the instantaneous values of the relaxation modulus at t = 0.03 s. Both representations show that the relaxation modulus increases nonlinearly as a function of the prestrain ε s .
The DMA tests are performed for the values of prestrain ε s = 0; 3; 4; 6 % at temperatures between −90°C to +90°C with an increment of 5°C. The imposed strain amplitude ε d for the DMA test was equal to 0.02 %, and, as a consequence, the ratio between the cyclic varying strain and the static strain load was in minimum 2 × 10 −3 .
The master curves obtained from the DMA data for the reference temperature of T ref = 20°C are presented in Fig. 5. Panel (a) represents the evolution of the norm of the amplitude of the complex modulus E * (ω) for increasing pulsation ω and different prestrains. One can remark that the influence of the prestrains occurs essentially at small pulsation. Moreover, a detailed analysis of the evolution of E * (ω) with increasing prestrains, at fixed pulsation ω = 31.4 rad/s, displayed in panel (b), exhibits a twofold increase of the amplitude of the complex moduli accompanied with a 30 % decrease of tan(δ) over the tested prestrain range. The discrepancy for low frequencies, smaller than 10 5 rad/s, is the direct effect of the prestrain. The different curves show that the amplitude of the complex modulus increases with respect to the prestrain value. As in the case of the relaxation modulus, the complex modulus amplitude increases as a function of prestrain ε s , but the phase angle δ decreases nonlinearly as a function of the same prestrain; see Fig. 5(b).
A complete master curve for the relaxation modulus relaxation in the absence of prestrain ε s = 0 % is exhibited in Fig. 6. The curve is an overlap of the relaxation data, which cover the large times, and the inverse Fourier transform of the DMA data, which cover small times. It will cover a total time range between 10 −15 s and 10 8 s for a reference temperature of 20°C. In panel (a) of Fig. 6, we displayed the relaxation curve, and one can observe a good match between the two data sets, in spite of a small slope change captured by the DMA relaxation data. Panel (b) of Fig. 6 exhibits a comparison between DMA and relaxation data in the evolution of the relaxation modulus with prestrain at t = 0.03 s. The data sets have a

Continuous model 3.1 Material without prestrain
Physically, for polymer material, the relaxation spectrum H (τ ) is related to the statistical distribution of the molecular mobility of polymer chains (Doi 1974). At the microscopic scale, the molecular mobility can now be measured by nuclear magnetic resonance (NMR) spectrometry to investigate the spin relaxation times of the chains. Moreover, the results from NMR published in Mowery et al. (2005), Litvinov and Steeman (1999) prove the existence of a multiphase matrix in reinforced polymers and explore the influence of the binder matrix interactions on the mechanical properties.
At the macroscopic scale, this spectrum can be computed directly from the complete master curve of the relaxation modulus using the inverse Laplace transform; see Eq. (1). Following Baumgaertel and Winter (1992), a linear relaxation modulus curve can be generally expressed by the following analytical expression: where E 0 , E g , t 0 , and β are adjustable parameters. Starting from the two remarks and applying the change of variable f (z) = H ( 1 z ) z in (1), the analytical expression of the spectrum H (τ ) is given by the inverse Laplace transform of (4): with Γ (β) the Gamma function of β. See Smith (1971) and Widder (1946) for a complete explanation of the operations performed.
An inspection of the curve of relaxation modulus of solid propellant without prestrain displayed in Fig. 6(a) shows that the relaxation modulus is not linear. Based on NMR measurements, Azoug et al. (2014cAzoug et al. ( , 2015 showed that propellant has a multiphased matrix due to his manufacturing process. Moreover, they were able to bring up three phases in the matrix of the propellant associated to families of polymer chains and mechanisms of mobility to explain the observed nonlinearity in the relaxation. Zimm and Lumpkin (1993) proposed to associate mechanism of molecular mobility to phases of the polymer matrix. As a consequence of the previous observation, we propose to associate to each a phase of the matrix of the propellant and correspondingly to each chain mobility mechanism a relaxation modulus of the form (4) and to sum up the contributions in the relaxation modulus: As a consequence, the modulus is approximated by the following analytic expression: From a physical point of view, as explained in Azoug et al. (2013a), this decomposition is coherent with the assumption of the existence of three families of different molecular mobilities in the binder: • The first family (subscript "1") represents the molecular mobility of the principal reticulated network polymer of the matrix. • The second family (subscript "2") represents the molecular mobility of the free polymer chains through the principal network of the matrix. • The third family (subscript "3") represents the long-time polymer chain mobility.
As the Laplace transform is a linear operator, from (7) it follows that the relaxation spectrum is the sum of the contributions for each phase H (τ ) = H 1 (τ ) + H 2 (τ ) + H 3 (τ ), or explicitly: Simhambhatla and Leonov (1995) assumed that the prestrain influences only the mobility of the free polymer chain, and the change induced into the fillers network will only affect the elastic part of the behavior, not the viscous part. Physically, one can visualize the effect of prestrain as a trapping mechanism of the fillers and the reticulated polymer network on the free chains, which have a growing restriction in the movement with increasing prestrain. Moreover, the NMR measurements on prestrained samples, reported in Azoug et al. (2015), show that the prestrain influences only the second family corresponding to free polymer chains, confirming the assumptions. This means that the prestrain affects the slope of the second contribution to the relaxation modulus, which implies that β 2 becomes a function of the prestrain ε s . The prestrain-modified spectrum H takes the following expression:  (7) master curve of the relaxation modulus and the evolution of the relaxation moduli of the three phases

Identification
The parameters of the three families, that is, (E gi , t i , β i ) 1≤i≤3 , are now identified from the experimental relaxation modulus curves in two steps: (i) identification of the parameters in the absence of prestress ε s = 0, that is, the parameter β 2 is considered constant, and (ii) identification of the prestress-dependent parameter β 2 (ε s ). Both steps are performed by minimizing the least square error between the model using the NonLinearModelFit operator of MATHEMATICA ® (Mathematica). The identified parameter values in the absence of prestress ε s = 0 are summarized in Table 1, and the comparison between prediction and measurements is plotted in Fig. 7. The model matches measurements over the total range of time, with the exception of the first decades, where a large error is observed. A closer inspection of plot and of the identified parameters shows that viscoelasticity is predominately driven by the molecular mobility of the first phase, representing the reticulated network. Nevertheless, in the time interval 10 −4 ≤ t ≤ 10 s, one can remark that the second phase of free polymer chains has a significant influence, even if its amplitude represents only 1 % of the amplitude of the principal reticulated network ( E g2 E g1 ≈ 0.01). The contribution of the third phase is very small when compared with the two others. However, it describes the viscoelastic behavior at long times.
The identified values of the model are now used to display in Fig. 8 the continuous spectrum of relaxation times and the contribution of the three phases of the polymer matrix. The phases representing molecular mobility will each contribute to a local maximum situation times τ of 1 × 10 −11 s, 1 × 10 −4 s, and 1.8 × 1 s, respectively. As expected, one can remark a similar trend with the spectra of spin-spin relaxation times of the chains obtained by NMR spectrometry.
The prestrain dependence of the complete behavior is embodied only in the parameter β 2 as mentioned earlier. The identification is performed using experimental data exhibited in  (8)) with a focus on the time between 10 −6 ≤ τ ≤ 10 6 in order to observe the second local maximum. Fig. 9 Evolution of β 2 in function of the prestrain: dots are identified from experimental results, a solid line is the proposed model (10) to catch this evolution Fig. 6(b), and a model is proposed, given by the following expression: The nonlinear regression leads to β 0 2 = β 2 (0) = 0.312 (see Table 1), ε r = 0.045, and q = 2.1. Figure 9 presents the comparison of the prediction of the model and experimental data, represented as a continuous line and dots, respectively.
The complete identification of the prestrain-dependent model permits to inspect the evolution of the spectrum of relaxation as a function of time and prestrain, that is, H (τ, ε s ). Figure 10 exhibits the time evolution of H (τ, ε s ) for the prestrain values ε s = 0, 2, 4, 6 %. The plot shows a finite support in time and a decreasing of the second protuberance with increasing prestrain. This protuberance represents the contribution of the second prestraindependent molecular mobility of free polymer chains. From the physical point of view, the decreasing contribution is the signature of the increasing trapping of the free chains with increasing prestrain.

Generalized Maxwell model with prestrain effect
Let us consider a viscoelastic body occupying the domain Ω. Using the standard notation in continuum mechanics, we shall denote the current point in the reference and the actual Fig. 10 Evolution of the spectrum of relaxation time between 10 −6 ≤ t ≤ 10 s for prestrain ε s = 0, 2, 4, 6 %. The second protuberance is decreasing with increasing prestrain configuration as X and x, respectively. The transformation gradient F is defined as F = ∂x ∂X . Under the assumption of small strain, the transformation gradient becomes F = I+Grad(u), where u = x − X stands for the field of displacements. The small strain tensor ε is computed as ε = 1 2 (Grad(u) + Grad T (u)). The viscoelastic constitutive law is presented in the framework proposed initially in Halphen and Nguyen (1975), Holzapfel (2006) and denoted as generalized standard materials. This formal framework is thermodynamical and assures the complete consistency with the principles of thermodynamics and conducts naturally to efficient and robust numerical integration algorithms. The constitutive law is defined in terms of the free energy ψ and the dissipation potential φ. For the viscoelastic bodies studied here, we consider a decoupled representation of the free energy ψ : where ψ ∞ and ψ i with 1 ≤ i ≤ n represent the elastic and viscoelastic parts of the free energy, and α i denote a series of internal variables; their physical significance will be clarified later in the presentation. In this framework, according to Coleman and Gurtin (1967), the total stress σ is defined as The formal construction assures that for any isothermal process, the Clausius-Duhem inequality with φ i the dissipation potentials. A detailed discussion of this point for the case of viscoelasticity is given, for example, in LeTallec and Rahier (1994).
In order to specify detailed expressions of the free energy and the dissipation potential, we make the following choice of a quadratic expression for the elastic part for the free energy: with 0 ≤ E ∞ and 0 ≤ E i . This choice leads to a classical model in linear viscoelasticity and is also presented, for example, in Simo (1987).
If we further assume that the internal variables α i represent the strain in the viscous dampers, then the fictitious stress q i acting on each damper will be associated to the internal variable α i . Similarly to the reasoning in Holzapfel (2006), by analogy with (12) we postulate that q i = − ∂ψ ∂α i . Therefore, applying (14) leads to The dissipation potentials are equally chosen as quadratic potentials: where the viscosity η i is a positive real number. The construction of the constitutive model does not restrict so far a strain-dependent viscosity parameter. A similar construction and extended discussion of the subject are presented in Amin et al. (2006). As a consequence, we define the constant relaxation time τ i as τ i = η i E i with η i the viscosity value at the origin of the process.
A series of simple algebraic transformations of Eqs. (12) and (14) lead to the following expression of stresses: Moreover, (13) and (16) lead to the following system of evolution equations for the internal variables α i : Complete equations of the constitutive model are obtained after a final inference, which starts from the assumption q ∞ = E ∞ ε, then uses the definition of the fictitious stress q i , and finally replaces the expression into (17) and (18). The final expression of the constitutive equations in terms of stresses and internal variables is The viscoelastic constitutive model described previously is denoted as a generalized Maxwell model. Moreover, it is consistent with the principles of thermodynamics by construction. It can also be described in terms of rheological elements by the so-called Prony series. The n viscoelastic elements are completely described by the characteristic parameters (τ i , η i ) 1≤i≤n , where τ i is the discrete relaxation time, and η i is the associated viscosity. It is generally schematically represented by simple rheological elements using springs and dampers; see, for example, (Findley et al. 1976;Simo and Hughes 1998) and the graphic representation in Fig. 11(a).
For a homogeneous one-dimensional specimen, the model can be easily integrated numerically. The backward Euler method leads to the following time stepping algorithm, where the stress σ m at the time instant t m is computed as a function of the stress σ m−1 at the preceding time instant t m−1 and a load increment expressed as the total strain increment in the following steps: This algorithm can then be generalized into full 3D evolution equations and implemented in standard finite element codes.

Results and discussion
The viscoelastic relaxation spectrum H (τ, ε s ) is a continuous function of the time τ already discussed. Moreover, the continuous spectrum accepts a discrete approximation in the form of the prestrain ε s dependent generalized Maxwell model, represented by the so-called Prony series. The identification and comparison of its predictions with experimental data will be discussed next.
Let us start with the continuous spectrum in the absence of prestrain H (τ, 0). The identification method of such a spectrum is performed in two steps: (i) the choice of the characteristic times of the viscoelastic elements and (ii) the determination of the viscosities.
The classical choice of characteristic times is an equidistant distribution on the logarithm scale as suggested, for example, in Baumgaertel and Winter (1992). Here we use an optimization method minimizing the distance between the continuous spectrum and its discrete approximation. The mathematical background of the method and its application has extensively been discussed in Jalocha et al. (2015a). In order to cover the features in this paper, we propose a short review in the Appendix.
If we consider the support of the continuous spectrum defined by a cut of 0.005 of its maximal value, 0.005 H max ≤ H (τ ) ≤ H max , then the corresponding time interval is defined for τ as 10 −14 s ≤ τ ≤ 10 20 s. The demanded accuracy in the optimization process for a relaxation spectrum H with three protuberances, one for each phase of the matrix, has led to n = 48 discrete characteristic relaxation times τ i . The Prony series is completed by the values of the viscosity η i obtained by the formula with r i = τ i+1 τ i−1 . The n elastic moduli E i are computed as the ratio of viscosities and characteristic times: To continue the identification method, we now introduce the prestrain dependence of the continuous spectrum H (τ, ε s ). The prestrain dependence is only carried by the viscosity η i of each element, and therefore, the prestrain-dependent viscosity naturally becomes Let us finally remark that the construction preserves prestrain-independent characteristic times, whereas the viscosity is prestrain dependent. The expression of the thermodynamic potentials preserve a quadratic free energy in terms of strains and is coherent with the principles.
A schematic view of the model is presented in panel (a) of the Fig. 11, and panel (b) illustrates the identified prestrain viscosity evolution of the 17th and 28th elements of the identified model. One can remark that the elements have different evolutions, which is also a signature of different nonlinearities in different viscoelastic elements.
The predictions of the identified prestrain-dependent generalized Maxwell model will be compared next with experimental data from the relaxation or the DMA tests.
The first comparison concerns the relaxation test starting from several values of prestrain, with an imposed strain history as displayed in Fig. 3(a). The value of the strain step was ε 0 = 1 %. The comparisons of model predictions and experimental data in terms of evolution of the relaxation modulus as a function of time are plotted in Fig. 12(a). The figure represents the time interval 0.03 ≤ t ≤ 1000 s, where the nonlinearity is strongest. The model matches both the evolution of the relaxation modulus as a function of the time and the order of magnitude imposed by the presence of different prestrains.
The second comparison concerns the DMA test. We only represent the evolution of the norm of the complex modulus and of the loss modulus as a function of prestrain at a given angular frequency for several prestrain values, with an imposed strain history as displayed in Fig. 3(b). The chosen pulsation and strain step values were ω = 31 rad/s and ε d = 0.02 %, respectively. The viscosity η i of the prestrain-dependent model is considered to follow the .  (23). We further assumed that the prestrain ε s can be reasonably approximated by the total strain ε.
A comparison, in terms of the norm of the complex modulus as a function of the prestrain, between numerical predictions and the measured experimental data is plotted in Fig. 12(b). One can easily observe that the model predictions match the experimental data. Moreover, the experimental increase of the amplitude of the complex modulus as a function of the prestrain is closely described by the model. The small panel in the same picture exhibits the capacity of the same model to correctly predict the decrease of the shift factor δ as a function of the prestrain. A close inspection of the values shows a small overestimation of the numerical prediction when compared with the experimental data.

Conclusion
We presented a series of relaxation and DMA experiments involving prestrain and proposed a thermodynamical coherent constitutive model to predict the measured data. The proposed model is a prestrain-dependent generalized Maxwell model. The prestrain dependence was carried by the viscosities of each viscoelastic element, whereas the characteristic relaxation times and elastic moduli stayed prestrain independent. The parameters of the model were then identified from experimental data with a robust identification method, and numerical predictions using a simple backward Euler integration algorithm allowed us to confront the model to the experiments.
The material under discussion was a high filled elastomer in the form of a solid propellant. It exhibited a high-volume fraction of the filler with respect to the matrix and, as a consequence, manifested a complex viscoelastic behavior.
The results show that the nonlinearity of the macroscopic material behavior exhibited during relaxation and DMA experiments could be matched by the predictions of the experiment. The nonlinearity discussed in this paper was the evolution of the viscoelastic properties as a function of the prestrain. The results presented here complete both in terms of experiments and modeling the initial results studied in Azoug (2010), Azoug et al. (2013b), Thorin et al. (2012).
The nonlinear evolution of viscoelastic behavior as a function of the prestrain is a direct physical signature of the mobility of the molecular chains in the polymeric binder. The experimental results discussed here confirm the NMR result obtained at a microscopic scale and reported in Azoug et al. (2015). Moreover, they also show that the prestrain influences the relaxation modulus up to time t = 10 −4 s.
The proposed prestrain-dependent generalized Maxwell model predicts closely the experimentally observed nonlinear viscoelastic behavior in terms of (i) complex and loss moduli or (ii) relaxation moduli. The thermodynamic consistent construction of the model is compatible with two or three spatial dimensions and can therefore easily be extended for multi-axial loading configurations. An interesting perspective is the exploration of the prestrain-dependent biaxial experimental. A further challenge is the programming of a complete numerical integration algorithm.
The influence of the prestrain on the second phase, characterizing the molecular mobility of the free polymer chains, seems to be a characteristic of highly filled elastomer. Preliminary finite element computations on 3D microstructures have shown a tenfold amplification of the strain inside the matrix when compared with the applied macroscopic strain. This strain amplification diminishes the space between the fillers and therefore restrains the movements of the free polymer chains. A work in progress analyzes the evolution of this effect as a function of the filler volume fraction.
This papers opens new perspectives both in terms of experiments and theory for the nonlinear viscoelastic behavior of solid propellant and, more generally, for highly filled elastomer under multiaxial prestrain.
The optimal distribution of the relaxation times τ i is defined in terms of the endpoints of each unit box function τ i = t i−1 +t i 2 , i = 1, . . . , n, for given n. Moreover, it will minimize the residual function R. The end points of the time series t 0 and t n are supposed to be a priori chosen. See Fig. 13 for a schematic view.