Gérald Munoz

Alain Dequidt
email: alain.dequidt@uca.fr

Nicolas Martzel
email: nicolas.martzel@michelin.com

Ronald Blaak

Armel Mbiakop-Ngassa

Julien Devémy

Benoit Latour

Sébastien Garruchet

Florent Goujon

Etienne Munch

Patrice Malfreyt

Predicting mechanical constitutive laws of elastomers with mesoscale simulations

The prediction of the mechanical behaviour of randomly reticulated polymer networks remains a challenging task despite the number of functional and physically motivated models that have been proposed to address it. We introduce here a bottom-up approach where the dynamics of the network is computed at the scale of its topological constraints and the interactions between them are built in order to retain the most relevant microscopic features of the polymer. Using input found in any classical polymer handbook, this model can accurately reproduce stress-strain curves of vulcanized polyisoprene.

Introduction

In the designing process of manufactured objects with target mechanical properties, an extensive use of finite element simulations is made. These simulations need to be supplied with material-specific constitutive relations expressing how local stress depends on local strain. In the search for the best material for a given application, it is desirable to predict such constitutive relations of new materials and therefore it is necessary to model the mechani-cal behavior of the material at a smaller scale.

In the case of polymer networks like gels and vulcanized rubbers, many models have been proposed, but so far they are not able to predict the mechanical properties with sufficient accuracy.

Phenomenological models are able to reproduce the mechanical response of specific materials. [START_REF] Attard | Hyperelastic constitutive modeling under finite strain[END_REF][START_REF] Pucci | A note on the Gent model for rubber-like materials[END_REF][START_REF] Yeoh | A new attempt to reconcile the statistical and phenomenological theories of rubber elasticity[END_REF] However, these models usually require a large number of fitting parameters, that often lack a physical interpretation and are therefore not able to quantitatively predict properties for new materials. [START_REF] Attard | Hyperelastic constitutive modeling under finite strain[END_REF][START_REF] Ogden | Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids[END_REF][START_REF] Hart-Smith | Elasticity parameters for finite deformations of rubber-like materials[END_REF] In addition, these models are usually poorly transferable to different loading modes: uniaxial (UN), equibiaxial (EB) or pure shear (PS). [START_REF] Marckmann | Comparison of Hyperelastic Models for Rubber-Like Materials[END_REF][START_REF] Hossain | More hyperelastic models for rubber-like materials: consistent tangent operators and comparative study[END_REF] On the other hand, some micro-mechanical models are supported by theoretical approaches based on polymer physics and statistical methods (tube models, extended tube, [START_REF] Kaliske | An Extended Tube-Model for Rubber Elasticity: Statistical-Mechanical Theory and Finite Element Implementation[END_REF][START_REF] Rubinstein | Elasticity of Polymer Networks[END_REF][START_REF] Wu | On improved network models for rubber elasticity and their applications to orientation hardening in glassy polymers[END_REF] ). In general, their mathematical formulations are quite complicated and the fitted parameters are often far from their physically expected values. [START_REF] Darabi | A generalized tube model of rubber elasticity[END_REF] These physical models also lack generality and are not able to reproduce stress-strain curves quantitatively in the whole range of relevant strain. [START_REF] Hossain | More hyperelastic models for rubber-like materials: consistent tangent operators and comparative study[END_REF] Finally, at the microscopic scale, simulations of the molecular dynamics of polymer chains can be performed with a level of detail that is capable of differentiating chemical specificity. However, the time and length scales of such simulations are still very far from the requirements of finite element simulations. [START_REF] Uddin | Multiscale modeling of a natural rubber: Bridging a coarsegrained molecular model to the rubber network theory[END_REF][START_REF] Tang | Allatom molecular dynamics study of impact fracture of glassy polymers. II: Microscopic origins of stresses in elasticity, yielding, and strain hardening[END_REF][START_REF] Li | Molecular dynamics simulation of fracture mechanism in the double interpenetrated cross-linked polymer[END_REF] Even molecular simulations using coarse-grained descriptions hardly reach microseconds and big enough system sizes to investigate these problems. [START_REF] Uddin | Multiscale modeling of a natural rubber: Bridging a coarsegrained molecular model to the rubber network theory[END_REF] In this paper, we use an intermediate scale to describe the polymer network. The model we introduce is built upon a previous model, [START_REF] Munoz | Heterogeneity Effects in Highly Cross-Linked Polymer Networks[END_REF] where only the topological constraints were taken into account. In order to capture the strain softening regime at small deformation, [START_REF] Rivlin | Collected Papers of R.S. Rivlin[END_REF] it is necessary to introduce entanglements as special topological constraints. [START_REF] Masubuchi | Brownian simulations of a network of reptating primitive chains[END_REF][START_REF] Everaers | Rheology and Microscopic Topology of Entangled Polymeric Liquids[END_REF][START_REF] Svaneborg | Connectivity and Entanglement Stress Contributions in Strained Polymer Networks[END_REF] The extension to add the modeling of entanglements explicitly is a major improvement with respect to the previous version of the model. Indeed, it provides access to the quantitative reproductions of stress-strain curves under various loading modes using only a few parameters, which all have a physical interpretations and reasonable values.

Our new model can be compared to other models at the same scale. The "Twentanglement" model, [START_REF] Padding | Uncrossability constraints in mesoscopic polymer melt simulations: Non-Rouse behavior of C120H242[END_REF] for instance, tracks the motion of every entanglement from birth to death, which is technically complicated and time consuming. Our model is similar, except that we include reticulation nodes, and we model entanglements in a more limited fashion by neglecting their creation and annihilation. This makes the model much simpler to implement and faster to run. The very popular slip-link model [START_REF] Edwards | The effect of entanglements in rubber elasticity[END_REF][START_REF] Tzoumanekas | From atomistic simulations to slip-link models of entangled polymer melts: Hierarchical strategies for the prediction of rheological properties[END_REF] is another model at a similar scale that takes entanglements into account. The drawback of this type of model is that it lacks a consistent force balance between pairs of entanglements (sliplinks), so that the topology of the network is not fully taken into account. In addition, it can be difficult to choose the appropriate slip-link properties (friction, spring stiffness,...) for new materials. This may be the reason why it is not used more often to form constitutive relations in continuum models. In the following, we will introduce the new polymer network model, analyze the impact of entanglements and compare the mechanical response with standard experiments from the literature.

Materials and Methods

The reticulated polymers are represented by a classical network with vertices and edges. Vertices represent topological constraints of the system, while edges correspond to the phantom polymer chains that connect pairs of topological constraints. We consider two types of topological constraints, chemical crosslinks and chain entanglements, that are modeled by two types of particles. The result is a network of crosslinks connected by phantom polymer strands that are subdivided into segments when passing through entanglement particles, similar to what is obtained by performing a primitive path analysis. [START_REF] Everaers | Rheology and Microscopic Topology of Entangled Polymeric Liquids[END_REF][START_REF] Svaneborg | Connectivity and Entanglement Stress Contributions in Strained Polymer Networks[END_REF] It is important to consider crosslinks and entanglements as different particles in order to allow a specific dynamics for both types of constraints. In the following, we call "strand" the portion of a polymer chain between successive crosslinks. We call "segment" the portion of a strand between successive entanglements or crosslinks.

From a numerical point of view, the system is a cubic simulation box which contains a network of particles connected with phantom segments. The parameters needed to build the system are the box size a, the number of cross-links and entanglements, their connectivity, the total number of Kuhn segments and their Kuhn length b. In practice, the box size is chosen large enough to have good statistics and the number of cross-links and entanglements particles are estimated from the experimental cross-link density and entanglement mass. The number of Kuhn segments is fixed by the polymer density and the value of the Kuhn length.

Results in this paper are produced with a box size a = 60 nm, 7776 crosslink particles and 20304 entanglement particles with a connectivity equal to 4, which results in a system that contains 15552 phantom strands. The Kuhn segments or monomers are not considered explicitly, but for every segment we use only their number as a variable quantity that determines the strength of the corresponding force in the simulations.

After generating the system and a dynamical equilibration procedure at a constant temperature T = 298.15 K, the box is stretched with a constant strain rate slow enough to ensure the validity of a quasi-static approximation. In order to check this point, we verified that the results are the same with a strain rate twice as slow. Depending on the mechanical solicitation mode, the contraction in one or two directions is made under the constraint of total fixed system volume.

Experimentally, the nominal stress σ 0 is obtained by the measurement of the macroscopic force F needed to stretch the elastomer and is normalised by the initial cross section A 0 of the sample, i.e., σ 0 = F/A 0 . The local and global nominal stress σ 0 can be computed at every step knowing the force f of each phantom segment by using the Irving-Kirkwood formula. [START_REF] Irving | The Statistical Mechanical Theory of Transport Processes. IV. The Equations of Hydrodynamics[END_REF] In practice, the segment contribution to the stress tensor σ is computed by:

σ seg = -1 V nodes i,j r i→j ⊗ f i→j ( 1
)
where V is the total volume of the simulation box, r i→j the vector from node i to node j, f i→j the force exerted by i on j and the sum runs over all node pairs connected by a segment. As we show in the appendix, σ 0 can be expressed as:

σ 0 = F A 0 = σ seg -σ seg ⊥ A A 0 , ( 2
)
where A is the instantaneous cross section of the sample. σ seg is the segment stress in the direction of traction while σ seg ⊥ is the segment stress in the direction of the free faces.

In order to characterize the macroscopic mechanical behavior we can examine the stress/strain relations. These curves provide two major quantities that characterize a material, the modulus at small deformation and the finite extension value λ max at which σ 0 diverges.

Network Generation

As demonstrated recently, the procedure by which the initial network is generated plays a prominent role with respect to a subsequent mechanical response. Here we have extended the process described before 15 in order to include entanglements.

This process is illustrated in Figure 1. The initial topology is build by placing nodes randomly, connecting the reticulation nodes, and entangling the segments. Thereafter, the topology is statistically equilibrated by swapping segment ends and transferring monomers between random pairs of segments according to the target probability:

P = P (r i |n i ) Q(N I ) P (r j |n j ) Q(N J ) , (3)
where r i is the end-to-end distance of segment i, n i is the number of Kuhn monomers of segment i and N I is the total number of Kuhn monomers of the strand to which segment i belongs. The conditional probability P (r i |n i ) is assumed to correspond to a Gaussian chain with the probability of over-extended chains set to 0 (i.e. if r i > n i b). Q(N I ) is the geometric probability for strand I to have N I monomers in total

P (r i |n i ) = 3 2π b 2 n i 3/2 exp - 3 r 2 i 2 b 2 n i , Q(N I ) = (1 -p) p N I -1 , p = 1 - 1 N I , ( 4
)
where b is the Kuhn length of the elastomer and N I is the mean number of Kuhn monomers between successive crosslinks. Some particular cases need to be excluded, namely those resulting in self-connected segments (loops) or in cycles of segments without reticulations (like a snake biting its own tail). The total end-toend distance of all segments is monitored and its convergence is used to decide when this step in the procedure can be terminated. Finally, a mechanical equilibration is performed by dis-placing the nodes according to the forces they experience and which are detailed in next paragraph.

Network dynamics

We use Brownian dynamics to model the displacements of the nodes. The velocity of a node i is given by:

γ v i = nodes j f j→i + 2 γ k B T η i , η i = 0 , η i η j = δ(t -t ) δ ij I , (5)
where the sum runs over all the nodes j connected to node i, η is a Langevin-like normalized random force, k B is the Boltzmann constant and T is the absolute temperature. The dynamics ( 5) is invariant under any transformation that leaves the product γ t invariant, such that results from properly done quasi static simulations are independent of γ. f is the finite extension version of the entropic force of a segment. Using the Cohen approximation of the inverse Langevin function, which is sufficient for the current purpose, [START_REF] Kröger | Simple, admissible, and accurate approximants of the inverse Langevin and Brillouin functions, relevant for strong polymer deformations and flows[END_REF] the force is given by:

f j→i = - k B T n ij b 2 3 -x 2 1 -x 2 r j→i , (6)
where x = r j→i /(n ij b) and n ij is the number of Kuhn monomers constituting a segment. Entanglements add some extra relaxation processes to the network. In the quasi-static limit all the segments of a strand between two crosslinks should experience the same tension. To this end a flow of Kuhn monomers is allowed from the more stretched segments of the strand to the less stretched ones. This flow corresponds to the sliding of the entanglements along the chain. Technically, we impose that the tension in every segment of a strand is the same and is determined by the total path length and the total number of Kuhn monomers of the strand.

In practice we will often have more than one entanglement node between two crosslinks.

During stretching, two such entanglements can get closer and stick together, effectively blocking each other from sliding further (Figure 2b). To avoid the formation of such entangled clusters, we have introduced the possibility for two entanglements to cross each other and keep on sliding along the chain, as illustrated in Fig- ure 2c. Two entanglement particles are assumed to cross each other when the scalar product of their next hypothetical distance vector r i→j and the current distance vector r i→j is negative

r i→j • r i→j < 0 . (7)
Note that no entanglements are created or destroyed when phantom segments move. This is consistent with the fact that the linking number is a knot invariant, if all the chains are connected at both ends. [START_REF] Horner | Knot theory in modern chemistry[END_REF][START_REF] Panagiotou | A Study of the Entanglement in Systems with Periodic Boundary Conditions[END_REF] It is important to realise that we propose an entanglement model while knowing that it is difficult to define and localise entanglements from an experimental point of view. The choice to model entanglements as particles makes it possible to connect results with those of a graph of an entanglement network as it can be obtained by for example a primitive path analysis. Whereas the Twentanglement model locates entanglements when they are created or annihilated, in our case they persist, but this is in agreement with the notion that in equilibrium the average number of trapped entanglements in a cross-linked network stays constant. From a topological point of view, entanglement crossing is not necessarily the most probable event, but it is not impossible either. Therefore when we allow all entanglement crossings along a strand between two crosslinks, we overestimate this relaxation mechanism for the phantom strands in the network. In real systems, however, it is also sometimes possible for entanglements to slide over a crosslink, which is not allowed in our model, and results in our model underestimating such a relaxation mechanism of the network structure. In any case this approximation will facilitate the relaxation of the network and results in a good agreement with the experimental data available. Pick one segment and one entanglement at random. e) Make the selected segment pass through the chosen entanglement (split in two segments). f) Repeat until all entanglements are saturated and distribute Kuhn monomers evenly between the segments. g) Pick two segments at random. h) Consider every reconnection and Kuhn monomer exchange and select one according to their relative probability. i) Repeat until convergence (see text). j) Displace both kinds of nodes to equilibrate forces.

Results and discussion

In order to study the network deformation, the simulation box is stretched in one or two dimensions by small constant increments while keeping the volume constant. The deformation increment is chosen small enough as to remain in the quasi-static limit. The stress tensor is computed as explained above. To investigate the effects of the different components of our model, we complement the simulation results (Sim) for the cross-linked network with entanglements and the possibility for entanglements to pass each other, with three variations. In system A we disable the possibility of entangle-ments to pass each other. In system B all entanglements have been changed into permanent cross-links and hence this system has a higher cross-link density, and in system C the entanglements are removed and the cross-link density remains the same as the reference simulation.

The Mooney-Rivlin representation 16,27,28 of the stress-strain curve is very useful to characterize the softening caused by the presence of entanglements (See Figure 3). A simple example to show the origin of the mathematical expression of the reduced stress is provided in the supplementary information. While entanglement-free polymer networks exhibit a monotonic Mooney-Rivlin stress-strain curve (systems B and C), networks with entanglements (curves Sim and A) display a nonmonotonic behavior and a negative slope (from right to left) in the low strain regime, in agreement with experiments. This demonstrates that this simple model of entanglements is enough to account for the experimental softening. We also note that for systems with the same crosslink density (curves Sim,A,C), the strain value at which the stress diverges is not really affected by the addition of entanglements. These two results suggest that entanglement nodes contribute less and less to the stress when the strain increases, which we believe is due to chain alignment.

In order to compare the results of our model with the reference experimental data from Treloar, 29 we used the parameters in Table 1. We show in Figure 3 how the model recovers this experimental trend with these parameters, along with the even more enhanced softening caused by the entanglements swapping mechanism in the curve labeled by Sim.

The crosslink density and entanglement mass have been chosen to give the best fit to the experimental data. The value of these parameters are in the range of what is found in the literature (See Table 1).

Figures 3 and4 compare our uniaxial tensile simulation results to experimental data. With only 4 parameters, our model matches well the experimental data, over the whole strain domain. An estimated error can be quantified by:

Error = 1 M M i=1 σ model 0 (λ i ) -σ expt 0 (λ i ) 2 (8
) with λ i the i th strain value and M the number of experimental data points available.

Figure 4 shows three types of mechanical deformation. Generally speaking, the model produces qualitatively a good estimate of the multi-axial behavior. Note that the same parameters are used for every deformation mode. The UN and PS simulations give a quantitative agreement in the whole strain domain. In the case of the EB simulations, the agreement with experiments is also very good at low strain, but the stress is slightly underestimated for strain values higher than 300 %.

To quantify the significance of the model, we compare our model with state of the art results produced by the Extended tube model [START_REF] Kaliske | An Extended Tube-Model for Rubber Elasticity: Statistical-Mechanical Theory and Finite Element Implementation[END_REF][START_REF] Darabi | A generalized tube model of rubber elasticity[END_REF] and the modified Flory-Erman model. [START_REF] Flory | Theory of elasticity of polymer networks. 3[END_REF][START_REF] Erman | Relationships between stress, strain, and molecular constitution of polymer networks. Comparison of theory with experiments[END_REF] In Table 2, we give a general indicator to compare our model to the two others. Note that the simulation parameters for every model are the set that gives the best fit of the uniaxial exten- A more detailed graphical comparison is made in Figures 5, 6 and7. In fact, Figure 5 highlights the conclusion that our model gives equivalent results compared to the two others models.

For the equibiaxial mode, Figure 6 complements the information of Table 2. It reveals the difference in behavior with respect to the extended tube model, which performs better than our model in the high strain limit, but overestimates the stress at λ < 3.

Finally, Figure 7 confirms that our model reproduces quantitatively the experimental shear data, whereas the modified Flory-Erman model produces a good qualitative approximation and the extended tube model overestimates the mechanical behavior for the full strain range.

Not only is the performance in predicting mechanical constitutive laws of elastomers by our proposed model favourable, it also sheds more light on the link between structure and mechanical properties of the cross-linked network. The model enables us to investigate the microstructure at the level of topological constraints, that can for instance be homogeneously distributed or more localised in clusters. In addition, it allows access to the topology of the network and in particular the identification of the most relevant nodes in the network where the stress is concentrated and which ultimately are the locations that weaken the material by ruptures.

An example of additional information that can be obtained from our model is the correlation between the force distribution and the number of Kuhn monomers per segment. Figure 8 presents this 2D correlation for a not entangled system and for an entangled network. First, it is interesting to observe that the Q segments (n) distributions are in both cases different from the theoretically expected distribution Q geom (n) (Red dotted lines on histograms in Figure 8). In the not entangled case, the decrease in the density of segments composed of less than ten monomers is explained by the restrictive condition of equation 3, which prohibits connecting two crosslinks at a distance r greater than the maximum length of a strand consisting of n monomers of length b (r > n b). Since the particles are randomly distributed in the box, it is rare for two crosslink particles to be close enough to be connected by a small number of monomers.

When we build our entangled systems following the procedure outlined in the section Network Generation, we obtain the same shape of the Q strand (N ) distribution as found for the fully crosslinked case (See Figure S-1). However, in the entangled system the Q segment (n) distribution after the dynamical equilibration, has an even smaller number of segments with fewer than ten monomers. This is due to our procedure of dynamic balancing of the entanglements, which allows a flow of monomers between segments belonging to the same strand. Indeed, the balance of tensions tends to homogenize the distribution of forces. Since the tensest segments are the shortest ones, their number decreases when monomers flow (see figure 8).

If we consider the correlation based on the two-dimensional histogram isodensity lines, the "shortest" segments in number of monomers experience higher forces in both cases. When we compare the not entangled and entangled case, we observe that the force distributions have the same mean value. However, for the entangled system this distribution is much less dispersed, which is to be expected since in our model the tension of the segments is the same for all the segments of the same strand, which results in more uniform forces. The two vertical 1D histograms of the force distribution seem to follow approximately a log-normal pattern. The force distribution in the entangled case appears more symmetrical as it is narrower.

Conclusion

In summary, we have introduced a model for a polymer network at the scale of topological constraints, including crosslinks and a simple and efficient implementation of entanglements. histogram is the same system, but with all entanglement particles replaced by cross-link particles. The bottom 2D histogram corresponds to the cis-polyisoprene system with parameters of Table 1.

As we have demonstrated, our model compares favorably to state-of-the-art mechanical models under various loading modes. It has the major advantage, that the four model parameters that are being used have a physical interpretation and for two of them values can be extracted from literature or experiments without fitting. This level of coarse-graining makes it possible to reach deformation speeds close to those of experiments and substantial sample sizes. Taking a few hours or days to obtain a stress-strain curve on 16 cpu, the current model is computationally slower than micro-mechanical or semiempirical models. However, it enables us to correlate the spatial structure of the topological constraints and the topology of the network to the macroscopic mechanical properties of the material. This makes it suitable for producing constitutive laws, that can be used in finite element simulations, with a predictive power. Even more realistic simulation results at large deformations are expected to be found by tuning the strand-length distribution and introducing a breaking mechanism for highly stretched strands, both of which form a natural extension of the current work.

Figure 1 :

 1 Figure 1: Sketch of the network preparation. Black dots are reticulation nodes, white dots are entanglements. For the sake of clarity of the pictures, the reticulation nodes have a connectivity of only 3 instead of 4 that is used in the simulations. a) Place nodes randomly. b) Choose a pair of reticulation nodes and connect them. c) Repeat until all reticulation nodes are saturated. d)Pick one segment and one entanglement at random. e) Make the selected segment pass through the chosen entanglement (split in two segments). f) Repeat until all entanglements are saturated and distribute Kuhn monomers evenly between the segments. g) Pick two segments at random. h) Consider every reconnection and Kuhn monomer exchange and select one according to their relative probability. i) Repeat until convergence (see text). j) Displace both kinds of nodes to equilibrate forces.

Figure 2 :

 2 Figure 2: a) Entanglements sliding towards each other. They may get blocked (b) or pass over (c).

Figure 3 :

 3 Figure 3: Uniaxial Mooney-Rivlin plot with σ 0 the nominal stress and λ the strain value. Exp: Experimental data; 29 Sim: Simulations; A: Same as Sim without entanglement swapping; B: Same as Sim with entanglement nodes replaced by crosslinks; C: Same as Sim with entanglement nodes removed.

Figure 4 :

 4 Figure 4: Simulation results with parameters adjusted to reproduce the experimental data by Treloar. 29 UN: uniaxial, EB: equibiaxial, PS: pure shear. All the simulations are performed with a single set of parameters obtained by fitting to uniaxial data.

Figure 5 :

 5 Figure 5: Comparison of the performance of the different models with respect to the uniaxial experimental data by Treloar.[START_REF] Treloar | Stress-Strain Data for Vulcanized Rubber under Various Types of Deformation[END_REF] Simulation results with parameters adjusted to reproduce the experimental data.[START_REF] Hossain | More hyperelastic models for rubber-like materials: consistent tangent operators and comparative study[END_REF]

Figure 6 :

 6 Figure 6: Comparison of the performance of the different models with respect to the equibiaxial experimental data by Treloar.[START_REF] Treloar | Stress-Strain Data for Vulcanized Rubber under Various Types of Deformation[END_REF] All the simulations are performed with the set of parameters obtained by fitting to uniaxial data.[START_REF] Hossain | More hyperelastic models for rubber-like materials: consistent tangent operators and comparative study[END_REF]

Figure 7 :

 7 Figure 7: Comparison of the performance of the different models with respect to the pure shear experimental data by Treloar.[START_REF] Treloar | Stress-Strain Data for Vulcanized Rubber under Various Types of Deformation[END_REF] All the simulations are performed with the set of parameters obtained by fitting to uniaxial data.[START_REF] Hossain | More hyperelastic models for rubber-like materials: consistent tangent operators and comparative study[END_REF]

Figure 8 :

 8 Figure 8: 2D correlation histograms between segment force distribution and the number of Kuhn monomers per segments in equilibrium before stretching. 1D histograms are the projections of the 2D histograms. The 1D horizontal histogram is the distribution of Kuhn monomers per segment. The red dashed curve is the theoretical Q(n) distribution. The 1D vertical histogram is the distribution of the force per segment. The black dashed curve is a fit by a log-normal distribution. The top 2Dhistogram is the same system, but with all entanglement particles replaced by cross-link particles. The bottom 2D histogram corresponds to the cis-polyisoprene system with parameters of Table1.

Table 1 :

 1 Parameters for cis-polyisoprene vulcanized specimens containing 8 % of sulfur by mass of elastomer. × 10 19 cm -3 2 to 8 × 10 19 cm -3 entanglement mass 2919 g mol -1 3680 g mol -1

		simulation	expected
	density	0.91 g cm -3	0.91 g cm -3
	Kuhn length	0.934 nm	0.934 nm
	crosslink density	3.6

Table 2 :

 2 Deviation from experimental data computed using(8) for our simulations and state-of-the art models[START_REF] Hossain | More hyperelastic models for rubber-like materials: consistent tangent operators and comparative study[END_REF][START_REF] Kaliske | An Extended Tube-Model for Rubber Elasticity: Statistical-Mechanical Theory and Finite Element Implementation[END_REF][START_REF] Flory | Theory of elasticity of polymer networks. 3[END_REF][START_REF] Erman | Relationships between stress, strain, and molecular constitution of polymer networks. Comparison of theory with experiments[END_REF]

	λ < 6 Our results Modified Flory-Erman Extended tube model
	UN	0.0370	0.0356	0.0580
	EB	0.1707	0.2185	0.0858
	PS	0.0128	0.0693	0.0944
	λ < 3 Our results Modified Flory-Erman Extended tube model
	UN	0.0299	0.0361	0.0275
	EB	0.0428	0.0711	0.0953
	PS	0.0123	0.0433	0.0652
	sion data, and which are used for every other
	deformation mode. Generally speaking, for the
	whole strain range (λ < 6), our model is equiva-
	lent to the two references. But if we restrict the
	comparison to values of λ < 3, the Error values
	as computed in Table 2 indicate that the re-
	sults of our model are significantly better than
	the best existing models from the literature, i.e.
	the model by Flory and Erman 30,31 (4 param-
	eters and modified by the 8-chain energy func-
	tion) and the Extended tube model 8,11 (3 to 6
	parameters) in the regime of low and interme-
	diate strains.

Acknowledgement This work was performed in SimatLab, a joint public-private laboratory dedicated to the modeling of polymer materials. This laboratory is supported by Michelin, Clermont Auvergne University (UCA), SIGMA Clermont and CNRS. We are grateful to the Mésocentre Clermont Auvergne University for providing computing and storage resources.

Appendix Computation of the total stress

The total stress in the simulation box has two origins:

• the first one, denoted by σ seg , is the contribution of the entropic elastic forces stretching the segments. These forces are purely attractive resulting in a negative contribution to the total pressure.

• the second one, denoted by σ vol , comes from the incompressibility of the system. This term is isotropic and corresponds to a Lagrange multiplier associated with the constraint of constant volume.

The value of σ vol is determined by the condition that the total stress on the free surfaces vanishes, which imposes that

Therefore the total stress in the direction of traction is

Note that the total isotropic pressure 1 3 Tr(σ) is different from σ vol because Tr(σ seg ) = 0 (in fact, < 0). Note also, that there are different ways to achieve a given deformation, for example by stretching along x only, or by compressing along y and z only. The final state has the same σ seg but different σ vol , resulting in different total stress σ and different pressure.

TOC Graphic