Antoine Gloria
email: agloria@ulb.ac.be

Patrick Le Tallec
email: patrick.letallec@polytechnique.edu

Marina Vidrascu
email: marina.vidrascu@inria.fr

FOUNDATION, ANALYSIS, AND NUMERICAL INVESTIGATION OF A VARIATIONAL NETWORK-BASED MODEL FOR RUBBER

may come

Foundation, analysis, and numerical investigation of a variational network-based model for rubber

Introduction

Two theories at different scales allow us to describe rubber-like materials. The first and more classical one is the continuum mechanics framework, ideal rubbers being prototypical of hyperelastic materials in finite deformation. At the other end of the spectrum, we find statistical mechanics approaches for polymer chains. In the continuum theory of rubber elasticity, constitutive laws are phenomenological in essence, see for instance the survey article [START_REF] Arruda | Constitutive models of rubber elasticity : a review[END_REF] (even if somewhat guided by the mathematical analysis, see [START_REF] Ciarlet | Mathematical elasticity. Volume I: three-dimensional elasticity[END_REF]). In statistical mechanics, the free energy of a polymer chain can be derived from first principles, as in [START_REF] Flory | Statistical mechanics of chain molecules[END_REF]. It is therefore tempting to relate the adhoc energy densities encountered at the continuum level with the physically motivated (free) energy of polymer chains. The difficulty of such a program is to pass from the free energy of one polymer chain to the free energy of a network of cross-linked polymer chains.

A first very strong assumption one usually makes is that the free energy of a network of polymer chains is the sum of the free energies of the chains as if they were isolated. This is of course not true and this assumption should be justified. Besides this assumption, one needs to specify how these chains interact together. To this aim, some constitutive assumption is made to pass from the polymer chain to the polymer network. The easiest and first assumption is due to Treloar in [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF] and consists in imposing that the network deforms locally in an affine manner. Hence all the polymer chains are deformed according to the macroscopic strain gradient. This is the so-called affine assumption (or in another context the Cauchy-Born rule). The associated energy density is in good agreement with mechanical experiments in small deformation. In large strain, however, this model overstimates the stress. This model overestimates the energy as well, and a natural way to make corrections to this model is to relax the affine assumption. A first possibility consists in postulating a geometric response of a representative volume element of polymer chains (yet somewhat unrelated to the actual polymer network). The eight-chain model by Arruda and Boyce in [START_REF] Arruda | A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF] is of this type. This model can be seen as a "cheap" relaxation of the Treloar model, and is indeed in much better agreement with real experiments (with the same very limited set of free parameters). Another way to relax Treloar's model is to replace the affine assumption by a "minimization principle": Instead of following the macroscopic deformation, we assume that the cross-links minimize the overall free energy Date: October 26, 2012.

1 of the network. This is what we call the variational model. It shares features with a model introduced by Böl and Reese [START_REF] Böl | Finite element modelling of rubber-like polymers based on chain statistics[END_REF].

In the present work we study in detail the variational model, and compare it both to the Treloar model (which it relaxes) and to the Arruda-Boyce model (which is another relaxation of the Treloar model). The aim of this work is to investigate how close a model relying only the very few physical ingredients of the Treloar model can bring us to rubber elasticity. Of course such a point of view excludes the comparison to other very interesting (and more advanced) models. This study also allows us to identify what is missing in the variational model. This can serve as a starting point to enrich the variational model using the insight of more advanced models such as the non-affine microsphere model [START_REF] Miehe | A micro-macro approach to rubber-like materials -Part I: the non-affine micro-sphere model of rubber elasticity[END_REF], the maximal advance path constraint model [START_REF] Tkachuk | The maximal advance path constraint for the homogenization of materials with random network microstructure[END_REF], the non-affine tube model [START_REF] Heinrich | Theoretical and numerical formulation of a molecular based constitutive tube-model of rubber elasticity[END_REF], or the constraint junction models [START_REF] Flory | Network topology and the theory of rubber elasticity[END_REF] -which all take into account topological effects.

The paper is organized as follows. In Section 2 we describe the Treloar, the Arruda-Boyce, and the Böl-Reese models. The variational model is precisely defined in Section 3. We then address the main mechanical and mathematical properties of the models in Section 4. For the Treloar and the Arruda-Boyce models, analytical formulas can be written in terms of the principal invariants of the Cauchy-Green strain tensor. For the variational model, one has to appeal to numerical simulations. The numerical solution method and some convergence issues are analyzed in Section 5. The dependence of the variational model upon some physical parameters is also investigated in this section ; this is the core of this article. The last section is dedicated to a quantitative comparison of the different models on Treloar's data, starting with the same characteristic parameters at the polymer-chain level.

Throughout the text we use the following notation:

• {e 1 , e 2 , e 3 } denotes the canonical basis of R 3 ; • M 3 is the set of 3 × 3 real matrices (strain gradients in R 3 ), and M 3 + is the subset of those matrices with positive determinants; • For all ξ ∈ M 3 , I 1 (ξ), I 2 (ξ), and I 3 (ξ) denote the first, second, and third principal invariants of ξ T ξ, respectively, that is:

I 1 (ξ) := tr(ξ T ξ), I 2 (ξ) := tr(cof(ξ T ξ)), I 3 (ξ) = (det ξ) 2 ;

• SO(3) is the set of rotations of R 3 ;

• S 2 is the unit sphere in R 3 ;

• D is an open bounded domain with a Lipschitz boundary;

• For all 1 < p < +∞, L p (D), W 1,p (D), and W 1,p 0 (D) denote the Lebesgue space of p-integrable functions on D, the Sobolev space of p-integrable functions on D whose distributional gradients are also p-integrable in D, and the Sobolev space of such functions which additionally vanish on the boundary ∂D of D, respectively.

Review of some discrete models for rubber

Before introducing the discrete models, we first quickly review the polymer chain free energies usually used as a fundamental brick in the polymer chain models.

2.1. Polymer chain free energy. Given a polymer chain made of N rigid segments of length l at absolute temperature β = 1 k B T , with a chain density n, the free energy (of entropic origin) for a chain of length r c can be modeled by

Wc (r c , N ) = n β N r c N l θ r c N l + log θ rc N l sinh θ rc N l - c β , (1)
where c is a constant and θ the inverse of the Langevin function L : α → coth α -1 α . The energy [START_REF] Alicandro | Mathematical derivation of a rubber-like stored energy functional[END_REF] has been derived by Kuhn and Grün in [START_REF] Kuhn | Beziehung zwischen elastische Konstanten und Dehnungsdoppelberechnung Eigenschaften hochpolymerer Stoffe[END_REF] using a non-Gaussian assumption: each segment of the chain obeys a non-Gaussian random walk. We refer to [START_REF] Flory | Statistical mechanics of chain molecules[END_REF] for details. In particular, the energy is infinite as soon as r c > N l, the total length of the chain. For discrete to continuum derivations, θ is usually replaced by the first terms of its series expansion:

θ(r) = 3r + 9 5 r 3 + 297 175 r 5 + 1539 875 r 7 + 672 359 r 9 + O(r 11 ), (2)
although this simplification is not essential for our discussion (e.g. Padé approximations behave better close to the finite extensibility limit). A series expansion of Wc then reads:

Wc (r c , N ) = n β N 3 2 r c N l 2 + 9 20
r c N l 4 + 9 350
r c N l 6 + 81 7000
r c N l 8 + 243 673750 r c N l 10 + O r c N l 12 . (3)
The behavior of the polynomial approximation at infinity satisfies the classical coercivity assumption on hyperelastic materials at infinity. Replacing the inverse of the Langevin function by the first terms of a series expansion is a rather good modeling at high temperature (see [START_REF] Tallec | Mécanique des Milieux Continus[END_REF]). A remarkable property of such an energy is Wc (0) = 0 and Wc (1) > 0.

In particular the prefered configuration of a polymer chain satisfies r c = 0.

When N is fixed, we simply write Wc (r c ) instead of Wc (r c , N ).

2.2.

Treloar type models. We consider a variant of the orginal model by Treloar. As in [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF], we assume that the polymer network deforms in an affine manner with deformation gradient ξ. We also assume that the network is isotropic and we give ourselves a distribution function ν(N ) of the polymer chains made of N segments in the network. We then associate general energies with the polymer chains, such as the one recalled in the previous subsection. To make precise the reference configuration of the chains, we use a distribution function ρ(L, N ) which gives the probability that a chain made of N segments has length L in the reference configuration. This yields for all ξ ∈ M 3 ,

WT (ξ) = R + R + S 2 Wc (L|ξ • e(X)|, N )dσ(X)dρ(L, N )dν(N ),
where S 2 is the unit sphere centered at the origin O, e(X) is the vector OX for X ∈ S 2 , and σ is the uniform measure on S 2 .

In what follows, for simplicity, we assume there is only one type of chains in the network (that is, all the chains have a fixed number N of segments), and we assume that the length of the polymer chains in the reference configuration is √ N l, which corresponds to the average distance of a random walker from the origin after N steps of length l, see [START_REF] Flory | Statistical mechanics of chain molecules[END_REF]. The Treloar energy density then turns into

WT (ξ) = S 2 Wc ( √ N l|ξ • e(X)|)dσ(X). ( 4
)
When Wc is replaced by its Taylor expansion (3) in (4), we denote by W 5 T the associated Treloar model (that is a Taylor expansion of WT ). As one can check on (4), W 5

T is minimal for ξ = 0. Actually one has not taken repulsion of chains into account in the model. To include volumetric effects in the description, one may add the so-called Helmholtz volumetric energy

W Helm (ξ) = K(det(ξ) 2 -1 -2 log(det(ξ))), (5)
which yields the final form we will consider:

W 5 T (ξ) = W 5 T (ξ) + W Helm (ξ)
. This can be generalized to any order p of the Taylor expansion.

2.3. Arruda-Boyce type models. In their original paper [START_REF] Arruda | A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF], Arruda and Boyce have proposed a model which relates the energy of the material in terms of its deformation gradient through the free energy of polymer chains. This model is called the eight-chain model. It basically amounts to taking as a representative volume element for the polymer network the set of eight chains linking the center of a cube to its eight corners. For a deformation gradient ξ ∈ M 3 + , the cube is supposed to align itself according to the principal directions of the deformation (eigenvectors of ξ T ξ), and to deform itself according to the principal stretches (square-roots of the eigenvalues of ξ T ξ). The energy density of the material is then proportional to the energy of the chains in the deformed cube. Note that 3 , where I 1 (ξ) = tr(ξ T ξ) is the first principal invariant. For isochoric deformations, the energy density is thus given by

WAB (ξ) = n β N   I 1 (ξ) 3 √ N l N l θ( I 1 (ξ) 3 √ N l N l ) + log   θ( I 1 (ξ) 3 √ N l N l ) sinh θ( I 1 (ξ) 3 √ N l N l )     = n β N   I 1 (ξ) 3N θ( I 1 (ξ) 3N ) + log   θ( I 1 (ξ) 3N ) sinh θ( I 1 (ξ) 3N )     . (6)
where typically n β = 0.27 MPa and N = 26.5 to fit Treloar's experiments. As for the Treloar model, when Wc is replaced by its Taylor expansion up to order p (see [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF] for p = 5), this gives rise to a Taylor expansion of (6) that we denote by W p AB . In order to consider non isochoric deformations, one should add the Helmholtz volumetric energy as for the Treloar type energies. We thus set

W p AB (ξ) = W p AB (ξ) + W Helm (ξ). ( 7
)
This variant of the Arruda-Boyce model is not the most commonly encountered. One usually uses reduced invariants:

W p AB (ξ) = W p AB ((det ξ) -1/3 ξ) + W Helm (ξ), (8)
which amounts to assuming that the energy density can be additively split into an isochoric part and a volumetric part. The advantage of such a decomposition is that the identity is the unique natural state of the energy density. Recalling however how the model has been built, it may seem more natural to take into account changes of volume in the representative volume element itself, as it is done in [START_REF] Böl | Finite element modelling of rubber-like materials -a comparison between simulation and experiment[END_REF]. Yet, due to the value of K in [START_REF] Arruda | Constitutive models of rubber elasticity : a review[END_REF] (which measures the strength of the volumetric energy), which may be typically 10 4 times larger than n β ∼ 0.1 MPa, the difference is not substantial in practice.

As we shall see in the last section, the Arruda-Boyce model is an elementary relaxation of the Treloar model, obtained using Jensen's inequality. This elementary observation allows to relate the Arruda-Boyce model to the Treloar model, and may give a better understanding of the nature of the modelling assumption: the network is allowed to relax by "convexification". The Arruda-Boyce model therefore belongs to the same family of models as the variational model we are interested in.

2.4. Böl-Reese model. The model developed by Böl and Reese in [START_REF] Böl | Finite element modelling of rubber-like polymers based on chain statistics[END_REF] is a hybrid model at the macroscopic scale. In particular, the energy of a macroscopic sample D is associated with a tetrahedral mesh of D: it is the combination of a volumetric energy on each element T of the mesh (typically the Helmholtz energy) and an energy associated with the edges of the mesh, which model "bundles" of polymer chains. The energy of one edge is given by a constant times the energy Wc (r) in [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF], where r = L L 0 , L 0 is the length of the edge in the reference configuration, and L, the length of the edge in the deformed configuration. The other parameters in Wc are then to be fitted. This model is hybrid since it does not enter the classical framework of continuum mechanics: the macroscopic energy is split into a continuum volumetric energy, and a discrete energy associated with edges of a mesh. As shown in [START_REF] Böl | Finite element modelling of rubber-like materials -a comparison between simulation and experiment[END_REF], numerical experiments are in good agreement with other classical models and experimental data. Yet this does not provide a pure continuum description of rubber elasticity. As proved in [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF], the model "converges" to a continuum model as the meshsize goes to zero. However, the limiting model highly depends on the geometry of the mesh in terms of isotropy properties for instance. This is an important issue which was not addressed in [START_REF] Böl | Finite element modelling of rubber-like materials -a comparison between simulation and experiment[END_REF][START_REF] Böl | Finite element modelling of rubber-like polymers based on chain statistics[END_REF].

The model we propose in the following section coincides with the limit of the Böl-Reese model as the meshsize goes to zero for a suitable mesh. It will therefore not be surprising that the two models yield similar results. Nevertheless, our modeling assumptions are different: we do not introduce bundles of polymer chains at a larger scale, we only work at the scale of the polymer chains.

Variational model: Motivation and rigorous results

3.1. Description and motivation. We consider a macroscopic sample of natural rubber D, whose boundary is linearly deformed through the map x → ξ • x, ξ ∈ M 3 + . The sample is made of a network of cross-linked polymer chains. The cross-links are assumed to be permanent. In this first (rough) model, we neglect entanglements, that is, we neglect topological constraints (this will be made clear in the definition of the network). Each polymer chain is itself made of a given number of monomers: for a chain the energy of a configuration is obtained through the probability density of a random walk (see for instance [START_REF] Kuhn | Beziehung zwischen elastische Konstanten und Dehnungsdoppelberechnung Eigenschaften hochpolymerer Stoffe[END_REF], [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF]). We assume that each monomer is surrounded by a fixed volume (from which other monomers are excluded), and that the network of chains is packed and almost incompressible. This assumption adds a volumetric term to the energy which depends on the configuration of the network. This volumetric term accounts for the interaction between the chains (which does not appear in the energy of one single chain). Note that the relevant scale associated with this contribution is much smaller than the one corresponding to the contribution associated with the random walk variable.

In the description of the rubber model we consider, we denote by u the positions of the cross-links, and by s = {s i } the positions of the monomers of the chain i. The Hamiltonian of the system can be split into two parts:

H(u, s) = H vol (u, s) + i H i (u, s i ).
The first part H vol (u, s) is the volumetric energy of the network, which models the interactions between the chains, whereas the second part H i (u, s i ) is the energy of each chain as if it were isolated (and for which u prescribes the end-to-end vector, and s i describes the positions of the monomers constituting the chain).

At finite temperature β = 1 k B T , the Boltzmann distribution yields the following formula for the free energy of a given deformed network:

F (ξ, D) = - 1 β ln Z = - 1 β ln U S i (u) exp -βH vol (u, s) - i βH i (u, s i ) du i ds i ,
where Z is the partition function, U is the set of admissible positions of the cross-links (satisfying the constraint on the boundary), and S i (u) denotes the set of admissible positions of the monomers composing the chain i whose head and tail are prescribed by u.

This free energy is far from being explicit. However, it is possible to further simplify the problem and still capture some interesting features. We present a heuristic argument which leads to the decoupling of the s i variables. We first assume that H vol (u, s) = H vol (u) only depends on u and not on s, which amounts to replacing the excluded volume constraint around monomers by an excluded volume constraint between cross-links. Note that this is a rather strong assumption whose effect is to make chains interact via their cross-links only: this decouples the variables s i from one another. We may then rewrite the free energy as follows:

F (ξ, D) = - 1 β ln U exp -βH vol (u) +β i 1 β ln S i (u) exp -βH i (u, s i ) ds i du . (9)
This amounts to taking the free energy of the following coarse-grained Hamiltonian:

H ξ (u, β) := H vol (u) - i 1 β ln S i (u)
exp -βH i (u, s i ) ds i .

We are now in position to interprete the Treloar assumption, and define the variational model. The affine assumption consists in replacing the integral over the admissible configuration set U in (9) by the evaluation of the integrand at u(x) = ξ • x. The variational model is defined by taking the minimum of the coarse-grained Hamiltonian, instead of the free energy:

F (ξ, D) |D| ≃ inf u H ξ (u, β) |D| . (10)
Formally, this amounts to considering two temperatures: the physical temperature β for the polymer chains and an "artificial" temperature β ′ for the cross-links (which is considered small, β ′ → ∞):

lim β ′ →∞ - 1 β ′ ln U exp -β ′ H ξ (u, β)(u) du = inf u H ξ (u, β).
It is not yet clear whether [START_REF] Delone | A local test for the regularity of a system of points[END_REF] is a valid approximation. This problem is currentlty under investigation on a simplified model.

In terms of orders of magnitude, recall that polymer chains are typically 100nm long whereas the macroscopic sample is of the order of the cm, which yields a factor 10 5 . Hence, provided D is a macroscopic sample, (10) will be close to the "thermodynamic limit"

W V (ξ) := lim |D|→∞ F (ξ, D) |D| , (11)
where D properly invades R 3 (see for instance [START_REF] Ruelle | Statistical mechanics. Rigorous results[END_REF]). For such a limit to exist, the network of polymer chains should have some ergodic property: either the network has some periodic structure (yet we are not dealing with crystals), or the network should yield spatial decorrelations (in a statistical or stochastic framework) -although other less physically relevant properties could also be considered stricto sensu.

Such a limiting process has been studied in [START_REF] Alicandro | Mathematical derivation of a rubber-like stored energy functional[END_REF][START_REF] Alicandro | Convergence analysis of the Böl-Reese discrete model for rubber[END_REF]. There, the limit ( 11) is proved to exist provided the stochastic network satisfies some structural and statistical properties.

In addition, the link between boundary value problems associated with the network and boundary value problems with the energy density W V is also made rigorous, thus generalizing the validity of the thermodynamic limit when boundary conditions and forcing terms (loads) are taken into account. We recall the relevant results in the following subsection.

3.2. Homogenization of stochastic discrete systems. In order to present the results of [START_REF] Alicandro | Convergence analysis of the Böl-Reese discrete model for rubber[END_REF] on the homogenization of stochastic discrete systems, we need to define the notion of stochastic network, and make precise the associated energy functional.

Definition 1. We say that a stochastic point process L in R 3 (that is a sequence of random points in R 3 ) is admissible if:

• (regularity) There exist r ≥ r > 0 such that almost surely:

any two points of L cannot be closer than r, -any ball of radius r contains at least one point of L; • (stationarity) L and x + L have the same statistics for all x ∈ R 3 ; • (ergodicity) L is ergodic.

We further assume that the Delaunay tessellation T of R 3 into tetrahedra associated with L (that is, the vertices of T are given by L) is almost surely unique (see [START_REF] Delone | A local test for the regularity of a system of points[END_REF]).

For rigorous definitions of admissible stochastic lattices and their Delaunay tessellations, we refer the reader to [START_REF] Alicandro | Convergence analysis of the Böl-Reese discrete model for rubber[END_REF], and to the references therein. Let us also introduce a rescaling of L and T . For all ε > 0, we set L ε := εL, which satisfies Definition 1 with εr and εr in place of r and r, and with T ε := εT in place of T .

Given a tessellation T ε of R 3 , one may define the space S(T ε ) of continuous and piecewise affine deformations u ε on T ε . Such deformations u ε , which are characterized by their values at vertices, are such that their gradients ∇u ε are piecewise contant on T ε . In particular, for every element T ε (tetrahedron) of the tessellation T ε , det ∇u ε|T ε measures the ratio of volume between u ε (T ε ) and T ε .

We are now in position to associate an energy with any deformation field u ε ∈ S(T ε ), on an open bounded domain D of R 3 . We consider two contributions: an energy associated with the changes of length of the edges of the tessellation, and an energy associated with the changes of volume of the elements of the tessellation. More precisely, denoting by E the set of edges of T , we define the energy of u ε ∈ S(T ε ) on D by

F ε (u ε , D) = ε 3 e∈E,e⊂D/ε W nn |e 1 -e 2 |, |u ε (εe 1 ) -u ε (εe 2 )| ε|e 1 -e 2 | + T ∈T ,T ⊂D/ε |εT |W vol (det ∇u ε|εT ), (12)
where e = (e 1 , e 2 ) (e 1 and e 2 are the two vertices of the edge e), W nn : R + × R + → R + is the energy of the deformed edges, and W vol : R → R + is the volumetric energy. We make the following assumptions on W nn and W vol :

Hypothesis 1. There exist p > 1 and positive constants c, C such that for all r ≤ γ ≤ r, λ ≥ 0, and ξ ∈ M 3 ,

cλ p -1/c ≤ W nn (γ, λ) ≤ C(λ p + 1), (13) W vol (det ξ) ≤ C(|ξ| p + 1). ( 14
)
We then have the following convergence result (see [START_REF] Alicandro | Convergence analysis of the Böl-Reese discrete model for rubber[END_REF]Theorem 5]).

Theorem 1. For all ε > 0, let L ε and T ε be the rescaled stochastic point process and associated Delaunay tessellation of Definition 1. For every open bounded subset D of R 3 , we consider the energy F ε (u ε , D) defined by [START_REF] Flory | Network topology and the theory of rubber elasticity[END_REF] for u ε ∈ S(T ε ), and extended by +∞ on W 1,p (D) \ S(T ε ), for p > 1, W nn , and W vol as in Hypothesis 1.

Then the functional

F ε (•, D) Γ(L p (D))-converges on W 1,p (D) as ε → 0 to the functional F V (•, D) defined by F V (u, D) = D W V (∇u(x)) dx, (15)
where W V : M d → R + is quasiconvex, satisfies a standard growth condition of order p:

there exists C > 0 such that for all ξ ∈ M d , 1 C |ξ| p -1 ≤ W V (ξ) ≤ C(1 + |ξ| p ),
and it is given by the asymptotic homogenization formula

W V (ξ) = lim R→∞ 1 R 3 inf{F 1 (u, Q R ), u ∈ S(T ), u(x) = ξ • x if dist(x, ∂Q R ) ≤ 2r}, ( 16
) with Q R = (-R/2, R/2) 3 , almost surely.
In addition, if L is isotropic in the sense that for all rotations R ∈ SO(3), L and R(L) have the same statistics, then W V is an isotropic energy density.

Remark 1. Note that we can take spheres in place of cubes in the asymptotic formula [START_REF] Gloria | Numerical approximation of effective coefficients in stochastic homogenization of discrete elliptic equations[END_REF], which we actually do in numerical approximations.

Let us comment on Theorem 1. First, ( 16) is a mathematical version of [START_REF] Flory | Statistical mechanics of chain molecules[END_REF], that shows the existence of the thermodynamic limit. It does not only provide the existence of W V , but also gives some additional properties, such as quasiconvexity and isotropy. Second, from this Γ-convergence result (and its generalization to the case when Dirichlet boundary conditions are considered, see [START_REF] Alicandro | Convergence analysis of the Böl-Reese discrete model for rubber[END_REF]), we deduce that given a boundary value problem on D and a minimizer u ε of the discrete energy at scale ε > 0, u ε converges in L p (D) as ε → 0 to a minimizer u ∈ W 1,p (D) of the continuum energy with energy density W V , the same boundary conditions, and external loads. This fully justifies the passage to the limit as ε → 0 in the static setting: not only the energy does converge, but also the minimizers.

In the following subsection, we show how to apply Theorem 1 to the discrete model introduced in Subsection 3.1, in order to define the variational model and its energy density W V .

Discrete variational model.

In order to make use of the rigorous results recalled above, we need to put [START_REF] Delone | A local test for the regularity of a system of points[END_REF] in the form of [START_REF] Flory | Network topology and the theory of rubber elasticity[END_REF]. To this aim, we further assume that the network is a tessellation of D into tetrahedra whose edges are the polymer chains, so that we can write [START_REF] Delone | A local test for the regularity of a system of points[END_REF] in the form

F (ξ, D) ≈ inf u∈U {F ε 0 (u ε 0 , D)} ,
for some ε 0 > 0, some function set U related to ξ ∈ M 3 , and with F ε 0 as in [START_REF] Flory | Network topology and the theory of rubber elasticity[END_REF]. The parameter ε 0 is the intrinsic lengthscale of the polymer network, which is the length of a monomer, also denoted by l. An edge e of the tessellation at scale ε 0 is then supposed to be made of

N e ≃ |e| l 2 = |e| ε 0
segments (or monomers). Hence, if the edge e has length L after deformation, its free energy is given by

Wc (L, N e ) = n β N e   L N e l θ L N e l + log θ L Nel sinh θ L Nel   = n β N e   1 √ N e L |e| θ 1 √ N e L |e| + log θ 1 √ Ne L |e| sinh θ 1 √ Ne L |e|   .
This formula allows us to properly define W nn . First we rescale the network so that it is of order ε = 1. The rescaled network will be taken as a reference and denoted by T . The associated parameters in Definition 1 are r = 1 (each chain has at least one monomer), and r 10 (a chain has, say, at most 100 monomers, so that its averaged length is of order 10). Let now e = (e 1 , e 2 ) be an edge of T , and for ε > 0, let u ε ∈ S(T ε ), where T ε is the ε-rescaled version of T . Then, setting λ := |uε(εe 1 )-uε(εe

2 )| ε|e 1 -e 2 |
the edge dilatation and K = n β , the rescaled energy of the rescaled edge εe deformed by u ε is given by

W nn (|e|, λ) = KN e   λ √ N e θ λ √ N e + log θ λ √ Ne sinh θ λ √ Ne   , (17)
corresponding to the unscaled values of the energy of the polymer chains when T is ε 0rescaled. For the volumetric energy, we take W vol (det ξ) = W Helm (ξ), where W Helm is defined in [START_REF] Arruda | Constitutive models of rubber elasticity : a review[END_REF].

It remains to check that W nn and W vol satisfy Hypothesis 1. In fact, the estimate from above is satisfied by W nn provided we consider any order of the Taylor expansion of the inverse of the Langevin function (for instance [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF], in which case p = 10). For the volumetric term, there is a technical difficulty: the Helmholtz energy density (5) does not satisfy Hypothesis 1 since W Helm (ξ) blows up as det ξ → 0 (as it is desirable in nonlinear elasticity). Although the discrete homogenization theorem is expected to hold in this case as well, the current state-of-art does not allow to directly extend the result to this case.

In order to circumvent this difficulty, we have proceeded in two steps in [START_REF] Alicandro | Convergence analysis of the Böl-Reese discrete model for rubber[END_REF]. First we have regularized the behavior of W Helm for det ξ ≤ η for η small so that the new volumetric energy density W η vol satisfies Hypothesis 1. We have then appealed to Theorem 1, which yields the existence of a homogenized energy density W p V,η (with obvious notation for the order 2p of the Taylor expansion (3)). In a second step, we have let the regularization parameter η go to zero, and showed that W p V,η converges to some energy density W p V (in the sense that the associated energy functionals Γ-converge) which satisfies in addition of frame-invariance, quasiconvexity, and weak lower-semicontinuity of the integral functional, the desired behavior

W V (ξ) → +∞, if det ξ → 0.
We refer the reader to [2, Theorem 11] for details. In the rest of this paper, the variational energy density W V is defined as the energy density obtained by the double limit procedure of [2, Theorem 11] :

W p V (ξ) = lim η→0 lim R→∞ 1 R 3 inf{F 1 (u, B R ), u ∈ S(T ), u(x) = ξ • x if dist(x, ∂B R ) ≤ 2r}, with B R = {Rx : |x| ≤ 1}, F 1
given by ( 12), W nn given by the p-order expansion of ( 17), and W vol by the η regularization of W Helm .

Remark 2. As pointed out to us by François Lequeux, the assumption that the polymer network forms a tetrahedral tessellation of R 3 is very strong. In particular it leads to a much too high connectivity of polymer chains. Therefore, the model is expected to overestimate the rigidity of the polymer network. In addition, as noticed by Michael Rubinstein, the density of monomers in the model is not correct: we have imposed local incompressibility whereas our model is not at the packing limit. In order to correct this inconsistency, we let depend the value of K on the number of chains and on the volume of the sample. Last, as remarked by one of the referees, the choice of the energy associated with the network implies that the end-to-end distance of polymer chains coincides with the characteristic spacing between the cross-links, whereas it is known that this distance is much larger than this spacing in most elastomers, see [START_REF] Flory | Network topology and the theory of rubber elasticity[END_REF]. We could also add in the model constrained junctions restrictions as in [START_REF] Flory | Network topology and the theory of rubber elasticity[END_REF] to model diffuse entanglements of cross links. These issues are very important from a physical point of view, and essentially regard the geometry of the network. They will be addressed in detail in a future work [17].

General properties

4.1. Mechanical properties. We now prove that the energy densities W T , W AB and W V of the Treloar, Arruda-Boyce, and variational models all satisfy the following mechanical properties:

• hyperelasticity,

• frame-invariance,

• isotropy,

• blow up of the energy when det ξ → 0 + , • strong ellipticity. Hyperelasticity is obvious since the models only depend on the pointwise value of the strain gradient. Frame-invariance is a consequence of the dependence upon ξ T ξ for the three models (or upon distances at the discrete level).

Isotropy is postulated in (AB) since the representative element aligns itself with the principal directions of ξ. It is a consequence of the integration over all the directions for (T). For (V), one has to assume in addition that the stochastic network is isotropic in the mean (the expectation of the stochastic network is isotropic), which implies the isotropy of the energy density, as proved in [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF]Theorem 6].

The behavior of the energy when det ξ → 0 + is a direct consequence of the definition of W Helm for (T) and (AB), whereas for (V) it is a consequence of the double limit procedure.

The strong ellipticity of the energy density implies two important properties: the material is stable (that is, homogeneous deformations are strict local minimizers), and this property allows one to prove short-time existence results for the system of elastodynamics. We recall that an energy density W :

M 3 → R is elliptic if and only if for all ξ ∈ M 3 and all a, b ∈ R 3 , the function R → R, t → W (ξ + ta ⊗ b) is convex. We say that W is strongly elliptic if t → W (ξ + ta ⊗ b) is in addition α-elliptic for some α > 0 (that is provided W is twice- differentiable: for all ξ ∈ M 3 , there exists α > 0 such that a ⊗ b • D 2 W (ξ)a ⊗ b ≥ α|a ⊗ b| 2 for all a, b ∈ R 3 ).
Lemma 1. The functions WT and WAB are strongly elliptic on M 3 + . Proof. We split the proof into two steps.

Step 1. Convexity of Wc (•, N ). We first claim that for all N > 0, the function r → Wc (r, N ) defined in (1) is convex, increasing, and α-elliptic with constant α = 6n βN l 2 on R + . The argument goes as follows. Since the Langevin function La : t → coth t -1 t is concave and decreasing on R + , its inverse θ is convex and increasing. In addition, the function t → sinh t t is convex and increasing on R + . Hence,

r → sinh θ(r) θ(r)
is convex and increasing on R + as well. Likewise, since t → -log t is convex and increasing on R + * = (0, +∞),

r → log θ(r) sinh θ(r)
is convex and increasing on R + . The same argument holds for r → rθ(r).

In addition, its second derivative, given by r → rθ ′′ (r) + 2θ ′ (r), is bounded by below by 2θ ′ (r) ≥ 6 on R + in view of the convexity of θ and using the series expansion (2) of θ. Taking into account the coefficients in (1) proves the claim.

Step 2. Strong ellipticity of WAB and WT on M 3 + . We first treat WAB . From Step 1 and ( 6), we infer that WAB is convex on M 3 since ξ → I 3 (ξ) is convex as the Frobenius norm on M 3 . For all i, j ∈ {1, 2, 3}, we have

e i ⊗ e j • D 2 ξξ WAB (ξ)e i ⊗ e j = W ′ c ( I 1 (ξ))e i ⊗ e j • D 2 ξξ I 1 (ξ)e i ⊗ e j + W ′′ c ( I 1 (ξ))D ξ I 1 (ξ) ⊗ D ξ I 1 (ξ). Since e i ⊗ e j • D 2 ξξ I 1 (ξ)e i ⊗ e j = I 1 (ξ) -ξ 2 ij I 1 (ξ) 3/2
is positive for all i, j ∈ {1, 2, 3} for det ξ = 0, since W ′ c ( I 1 (ξ)) > 0 for ξ = 0, and since W ′′ c ( I 1 (ξ))D ξ √ I 1 (ξ) ⊗ D ξ √ I 1 (ξ) is an elliptic tensor (although not strongly) for all ξ ∈ M 3 , the energy density WAB is strongly elliptic on M 3 + . We now turn to W T and proceed similarly. For all e ∈ S 2 the function ξ → |ξe| is convex by the triangle inequality. Hence, ξ → Wc (|ξe|) is a convex function as well, so that (4) implies that W T is convex on M 3 as a convex combination of convex functions (for e describing S 2 ). We now prove the strong ellipticity on M 3 + . We have:

S 2 e i ⊗ e j • D 2 ξξ (|ξe(X)|)e i ⊗ e j dσ(X) = S 2 |ξe(X)| 2 -(e i • ξe(X)) 2 |ξe| 3/2 dσ(X)
which is positive provided ξ = 0. As above, this implies the strong ellipticity (of WT ) on M 3 + . In particular, W T and W AB are strongly elliptic as the sum of a strongly elliptic function on M 3 + ( WAB and WT , respectively) and of an elliptic function (the Helmholtz energy density, which is polyconvex and hence elliptic). The strong ellipticity of W V is a much more delicate question, since strong ellipticity can be lost by homogenization in nonlinear elasticity (see in particular the very inspiring contribution of Geymonat, Müller, and Triantafyllidis in [START_REF] Geymonat | Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity[END_REF]). Provided some technical assumptions (which are supported by numerical simulations) and using the specific structure of Wc , the first author has shown in [START_REF] Gloria | Strong ellipticity of nonlinear elastic materials and homogenization of periodic and stochastic discrete systems[END_REF] that W V is strongly elliptic on M 3 + . The following lemma shows that W T and W AB are coercive energy densities: Lemma 2. For all p ≥ 1, the energy densities ξ → W p AB (ξ) and ξ → W p T (ξ) are the sum of a convex function of ξ and of a convex function of det ξ. In addition, they satisfy the following coercivity property: there exists c p > 0 such that for all ξ ∈ M 3

min{W p AB (ξ), W p T (ξ)} + 1 ≥ c p (|ξ| 2p + | det ξ| 2 ). ( 18
)
Proof. The structure of W AB and W T with respect to convexity is a consequence of the convexity of WAB and WT and of the convexity of the function t → t 2 -1 -2 ln t on R + * (for the Helmholtz energy density defined in ( 5)).

The coercivity estimate [START_REF] Gloria | An optimal variance estimate in stochastic homogenization of discrete elliptic equations[END_REF] follows from the Helmholtz energy density, noting that

t 2 -1 -2 ln t ≥ max 0, t 2 2 -2 ,
and from the fact that the Taylor expansion (3) of Wc only has non-negative coefficients for all p ≥ 1. In particular,

W p AB (ξ) ≥ C AB (N, p) I 1 (ξ) 2p = C(N, p)|ξ| 2p ,
where C AB (N, p) is a positive constant depending only on N and p, and |•| is the Frobenius norm on M 3 . For the Treloar model, we slightly anticipate on Section 6, and appeal to [START_REF] Kuhn | Beziehung zwischen elastische Konstanten und Dehnungsdoppelberechnung Eigenschaften hochpolymerer Stoffe[END_REF], that is W T (ξ) ≥ W AB (ξ) for all ξ ∈ M 3 .

We conclude this subsection by the determination of the natural states of W T , W AB , and W V (that is, the absolute minimizers of these energy densities on M 3 + ). In the case when reduced invariants are used, as for the original Arruda-Boyce model W AB defined in [START_REF] Böl | Finite element modelling of rubber-like polymers based on chain statistics[END_REF], the absolute minimizer is the identity. When standard invariants are used, as in the variants of Arruda-Boyce and Treloar W AB and W T considered here, the identity is not the natural state since there is a competition between the energy of the polymer chains (which is minimized at 0) and the Helmholtz energy (which is minimized at identity). This gives rise to: Lemma 3. The energy densities W AB and W T admit the same unique natural state on M 3 + , which is a dilatation αId with 0 < α < 1. Proof. This is a consequence of convexity. The functions W AB and W T are continuous on M 3 + and infinite at infinity and when det ξ → 0 + , so that they attain their minimum. Let α 3 denote the determinant of a minimizer. By frame-invariance it is enough to consider diagonal matrices of the form ξ λ 1 ,λ 2 ,α := diag(λ 1 , λ 2 , α 3 /(λ 1 λ 2 )). The function (λ 1 , λ 2 ) → WAB (ξ λ 1 ,λ 2 ,α ) is strictly convex on R 2 as the composition of an increasing strictly convex function with a strictly convex function. Hence it admits a unique minimum (λ * 1 , λ * 2 ). By symmetry, λ * 1 = λ * 2 , and therefore λ * 1 = λ * 2 = α, and the natural state of determinant α 3 is the dilatation αId. It remains to note that the function t → W AB (tId) is strictly convex, so that there is only one dilatation of minimal energy, and therefore only one natural state. The same arguments hold for W T , and the natural states coincide since W T and W AB coincide on dilatations.

A weaker result holds for W V , based on [26, Theorem A]: Lemma 4. [3, Theorem 8] If W V is isotropic (which follows from the isotropy of the stochastic lattice), it admits a natural state on M 3 + which is a dilatation. Note that uniqueness of the natural state is not ensured by this lemma (its proof only relies on argument which are compatible with the existence of several natural states, see [26, Section 3] e.g.).

Remark 3. All the results of this paragraph hold as well for W p AB , W p T , and W p V for all p ≥ 1 since the Taylor expansion (3) of Wc is strictly convex at any order.

Boundary value problems.

The main interest of constitutive laws is their application to mechanical problems, or in mathematical terms, to the associated boundary value problems. This section is devoted to the study of the well-posedness of boundary value problems in terms of minimization. Within this context, we only have to check classical coercivity and lower-semicontinuity properties ([6], [START_REF] Ciarlet | Mathematical elasticity. Volume I: three-dimensional elasticity[END_REF]).

In particular we have the following (standard) existence result for the minimization problems associated with Treloar, Arruda-Boyce, and variational models: Theorem 2. Let p ≥ 2, Γ 0 be a measurable subset of ∂D of positive measure, and let φ 0 : Γ 0 → R 3 be a measurable function such that Φ = {v ∈ W 1,2p (D), v = φ 0 on Γ 0 , det ∇v ∈ L 2 (D), det ∇v > 0 a. e. on D} is not empty. Let f ∈ L ρ (D) and g ∈ L σ (∂D \ Γ 0 ) be such that the linear form

L : v ∈ W 1,p (D) → L(v) := D f (x) • v(x)dx + ∂D\Γ 0 g(x) • v(x)da(x)
is continuous (da denotes the surface measure), let

F (v) := D W p (∇v(x))dx -L(v),
where W p is either W p T , W p AB , or W p V , and assume that inf v∈Φ F (v) < ∞. Then there exists at least one function u ∈ Φ such that

F (u) = inf v∈Φ F (v).
Proof. We first treat the case of W p T and W p AB using polyconvexity, and then turn to W p V using results of [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF].

Step 1. Proof for W p T and W p AB . This is a direct application of Ball's seminal results in [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]. Since 2p ≥ 3, for all v ∈ W 1,2p (D), det ∇v ∈ L 1 (D), and the set Φ := {v ∈ W 1,2p (D), det ∇v ∈ L 2 (D)} is weaklyclosed (see for instance [9, Section 7.6]). By Lemma 2, F is coercive on Φ and lowersemicontinuous for the weak-convergence in Φ (see [START_REF] Ciarlet | Mathematical elasticity. Volume I: three-dimensional elasticity[END_REF]Section 7.7]). Let u n ∈ Φ be a minimizing sequence of F on Φ. Up to extraction, u n weakly converges to some u ∈ Φ, and

F (u) ≤ inf v∈Φ F (v).
It remains to prove that u ∈ Φ, which is a consequence of the property lim det ξ→0

W p (ξ) = +∞, and of the compactness of the trace operator, see for instance in [9, Section 7.7] for both arguments.

Step 2. Proof for W p V . The structure of the proof is the same as above, noting that the lower-semicontinuity of F for the weak convergence in W 1,p (D) is a direct consequence of Γ-convergence. We refer the reader to [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF]Theorem 11] for details.

Numerical approximation of the variational model

5.1. The numerical solution method. Our aim is to approximate ξ → W V (ξ) with the help of (a variant of) formula ( 16):

W V (ξ) = lim R→∞ 1 |B(R)| inf{F 1 (u, B R ), u ∈ S(T ), u(x) = ξ • x if dist(x, ∂B R ) ≤ 2r}, ( 19
)
with W nn given by ( 17), W vol = W Helm given by ( 5), and where B(R) is the ball of radius R > 0 centered at the origin. For the numerical approximation procedure we take W vol = W Helm , and do not use further regularization as in [START_REF] Alicandro | Convergence analysis of the Böl-Reese discrete model for rubber[END_REF]Theorem 11] -this owes to the fact that this is the case of interest although we are not able to directly deal with W vol = W Helm in the analysis yet.

The physical parameters to be fixed are the typical length e of the edges of the network (that is, the average of the lengths of all the edges of the tessellation), the typical number of monomers N per chain (related to the characteristic length l of a monomer through the identity N = e l 2 ) and K in (17), and K in [START_REF] Arruda | Constitutive models of rubber elasticity : a review[END_REF]. Typical values in the simulations are chosen in function of R so that the ratio of the contributions of the chains and of the volumetric effect is constant (this does not change the value of W V at the limit). We denote by N R and by |B(R)| the number of edges and the volume of B(R), respectively. We will take:

N ∈ [25, 250], K = |B(R)| N R × 26.5 N × 0.27 MPa, K ∈ {5, 50, 100, 500} MPa.
This choice of K makes the contribution to the energy of the N R N monomers be of order |B R |0.27 MPa (such a rescaling makes the number of monomers per unit volume to fixed in the simulations -which is needed to compare polymers with different N since the network is not at the packing limit, see Remark 2). These correspond to the typical values considered in the Arruda-Boyce model to fit Treloar's experiments.

For the numerical approximation procedure one has to pick a large -though finite -R ≫ 1, and minimize F 1 (•, B R ) on the set of continuous and piecewise constant functions on T ∩ B R whose values on the boundary coincide with the affine deformation x → ξ • x.

To this aim, one needs to know T ∩B R . Yet, T is a Delaunay tessellation associated with a point process in the whole space R 3 . We therefore need to approximate the point process itself on the domain B R . The point process we consider on R 3 is the thermodynamic limit of the random parking measure associated with unit balls (that we will refer to as the "parking lattice"), which is rigorously defined as the limit of point processes on finite domains in [START_REF] Penrose | Random parking, sequential adsorption, and the jamming limit[END_REF]. As shown in [START_REF] Gloria | Random parking, Euclidean functionals, and rubber elasticity[END_REF], the parking lattice is stationary, ergodic, isotropic, almost surely general (the associated Delaunay tessellation is unique), and satisfies the hardcore and non-empty space conditions required by Theorem 1. Hence, formula [START_REF] Gloria | An optimal error estimate in stochastic homogenization of discrete elliptic equations[END_REF] makes sense. The approximation of the parking lattice on finite domains B R is as follows.

For all R > 0, we make a uniform mesh of the sphere S R of radius R (with triangles of side ∼ 1) and consider a hardcore Poisson point process (with minimal distance 1) in B R-1 up to the packing limit (that is, until one cannot add any other point: due to the hardcore constraint, there is an easy upper bound on the number of points which can be accepted).

In particular, for such a point process, any two points are at least at distance 1, and any ball of radius 1 has at least one point. We denote by T R a Delaunay tessellation associated with the mesh of S R and the points in B R-1 . As proved in [START_REF] Gloria | Random parking, Euclidean functionals, and rubber elasticity[END_REF], although T R does not coincide with T ∩ B R , we still have

lim R→∞ 1 |B(R)| inf{F 1 (u, B R ), u ∈ S(T R ), u(x) = ξ • x if dist(x, ∂B R ) ≤ 2r} = W V (ξ). ( 20
)
This is the final variant of ( 16) we consider, and which has the advantage to be practically computable.

The numerical approximation of ( 20) is made in two steps:

• We first generate the deterministic set of points on ∂B R and a realization of the stochastic set of points in B R-1 . The latter is generated iteratively. Points are randomly picked in B R-1 . The first point is accepted. When another point is picked, either it is at distance less than 1 from a point which has already been accepted and it is discarded, or it is at distance at least 1 from all the other points and it is accepted. The algorithm stops when B R-1 is packed, that is, when no additional point can be accepted. Given the deterministic set of points on ∂B R and the realization of the random set of points in B R-1 we then construct an associated Delaunay tessellation of B R . Note that this choice sets the relation between the number N R of chains in the volume B(R), and N R /|B(R)| tends to a deterministic limit as R → ∞. • In a second step, we solve the minimization problem associated with [START_REF] Gloria | Random parking, Euclidean functionals, and rubber elasticity[END_REF] for R finite and the Delaunay tessellation of B R (well-defined as the minimization of a smooth coercive function on a finite-dimensional space) by a Newton algorithm, as it is classical in nonlinear elasticity (see for instance [START_REF] Tallec | Numerical methods for nonlinear three-dimensional elasticity[END_REF], and [START_REF] Vidrascu | Solution of non-linear elasticity problems using the continu software[END_REF]) provided the addition of the energy of the edges (which are "non-standard" one-dimensional elements). Continuation methods are also used to ensure the convergence of the Newton algorithm. In practice, we also consider several independent realizations of the stochastic set of point and make an empirical average of the approximations of W V obtained. This enhances the convergence with respect to the randomness. For the analysis of numerical methods to approximate homogenized coefficients for random discrete systems, we refer the reader to the linear case dealt with in [START_REF] Gloria | Numerical approximation of effective coefficients in stochastic homogenization of discrete elliptic equations[END_REF], based on [START_REF] Gloria | An optimal variance estimate in stochastic homogenization of discrete elliptic equations[END_REF][START_REF] Gloria | An optimal error estimate in stochastic homogenization of discrete elliptic equations[END_REF].

Note that with respect to standard three-dimensional finite element methods, the solution of the minimization problem associated with [START_REF] Gloria | Random parking, Euclidean functionals, and rubber elasticity[END_REF] for R finite and the Delaunay tessellation of B R requires to consider one-dimensional elements in addition to three-dimensional elements. These one-dimensional elements are represented by the edges of the tessellation.

For notational convenience, we define

W nn : R + × R d → R + by W nn (h, ζ) := W nn (h, |ζ|),
and we recall that for all ξ ∈ M + d , W vol (det ξ) = W Helm (ξ). At the end of the minimization algorithm, the homogenized energy W V is approximated by the spatial average on B R of the energy density of the minimizer u R which has been numerically obtained:

W V (Λ) ≃ 1 |B R | F 1 (u R , B R ) = 1 |B R | e∈E,e⊂B R W nn |e 1 -e 2 |, |u R (e 1 ) -u R (e 2 )| |e 1 -e 2 | + T ∈T ,T ⊂B R T W vol (det ∇u R ) = 1 |B R | e∈E,e⊂B R W nn |e 1 -e 2 |, u R (e 1 ) -u R (e 2 ) |e 1 -e 2 | + T ∈T ,T ⊂B R T W Helm (∇u R ) .
The Piola-Kirchhoff stress tensor is given by the spatial average on B R of the associated local Piola-Kirchhoff stress tensor (provided the minimizer is isolated, and the local Hessian strongly elliptic, see [START_REF] Gloria | A direct approach to numerical homogenization in finite elasticity[END_REF]Section 4.2] for related arguments in the continuum case):

D Λ W V (Λ) ≃ 1 |B R | e∈E,e⊂B R D ζ W nn |e 1 -e 2 |, u R (e 1 ) -u R (e 2 ) |e 1 -e 2 | ⊗ e 1 -e 2 |e 1 -e 2 | + T ∈T ,T ⊂B R T D ξ W Helm (∇u R ) .
Equivalently, we could compute the stress tensor as a residual in the weak formulation (which is slightly cheaper numerically).

A study of the convergence of the numerical approximation of W V (ξ) in terms of R and the number of realizations is proposed in the following paragraph for some typical values of ξ (both in small and large strains).

NR 130 1600 6000 16000 33000 59500 nR 2160 270 80 34 18 10

Table 1. Number of edges and associated number of realizations in the numerical computations.

Note that the numerical tests are performed using the natural state of W V as the reference configuration in Subsections 5.5 and 5.6, so that the identity has least energy.

Convergence properties.

In this subsection, we investigate the convergence properties of the numerical approximation of W V . The numerical analysis of the numerical approximation procedure is made particularly difficult because • the unknown is a vector and not a scalar,

• the system is nonlinear and nonconvex,

• the randomness is complex (the decorrelation properties of the random parking lattice are not easy to quantify). As seen in [START_REF] Gloria | An optimal variance estimate in stochastic homogenization of discrete elliptic equations[END_REF][START_REF] Gloria | An optimal error estimate in stochastic homogenization of discrete elliptic equations[END_REF][START_REF] Gloria | Numerical approximation of effective coefficients in stochastic homogenization of discrete elliptic equations[END_REF] on a simpler linear scalar case, the error splits into two parts: a random error and a systematic error. The random error is of variance type and measures the fact that the computed quantities are random variables themselves. This error can be made as small as desired by using a sufficiently large number n R of realizations. It is proved to decay at the central limit theorem scaling in terms of n R times the number N R of chains (in this scalar case, chains are viewed as conductances

) in B R , that is (n R N R ) -1/2 .
The systematic error is related to the boundary conditions: we impose linear Dirichlet boundary conditions, which perturbs the energy -no matter the number of realizations. This error also vanishes at the limit N R → ∞. Its scaling is expected to be of the order R -1 since this is a surface effect -although there is no rigorous proof of this fact.

The physical parameters for these tests are as follows

N = 26.5, K = |B(R)| N R × 0.27 MPa, K = 100 MPa.
We have conducted two series of tests with two different strain gradients A 1 and A 2 , the first in moderate deformation (∼ 25%), and the second in large deformation (∼ 300%):

A 1 =   1.1 0 0 0 1.2 0 0 0 25/33   , A 2 =   2 0 0 0 3 0 0 0 1/6   .
Both deformations are isochoric: det

A 1 = det A 2 = 1.
We focus on the first Piola-Kirchhoff stress tensor ∂W V (ξ) ∂ξ . We expect the stress tensors

Π 1 (N R , n R ) = [Π 1 (N R , n R )] ij and Π 2 (N R , n R ) = [Π 2 (N R , n R )] ij ,
associated with A 1 and A 2 respectiveley, to be diagonal. We therefore focus on the principal stresses. Each effective stress tensor is obtained as the empirical average over n R realizations of a stochastic network in B R (with approximately N R edges). These numbers are gathered in Table 1.

The (square root of the) variance of one realization is plotted for the three diagonal terms in function of N R on Figure 2 for A 1 and Figure 3 for A 2 , in log-log scale. The straight lines are linear fittings. Their slopes are approximately -1/2 (between -.45 and -.5), as expected.

To conclude, we show on Figures 4-6 and 7-9 the convergences of the diagonal terms of the Piola stress tensor (with the variance) in function of N R , for A 1 and A 2 , respectively. As can be seen, the systematic error error may largely dominate the random error (see in particular Figure 9). From now on, we shall consider that the approximation has converged for N R ∼ 100000 and n R ∼ 10. 3. Diagonal entries for the rotations of A 2 5.3. Isotropy of W V . The isotropy of the model is directly related to the isotropy of the network. The physical parameters for these tests are as follows:

N = 26.5, K = |B(R)| N R × 0.27 MPa, K = 100 MPa.
To check the isotropy, we consider the two strain gradients A 1 and A 2 that we transform by the following three rotations:

R 1 =   0 0 1 1 0 0 0 1 0   , R 2 =   0 1 0 0 0 1 1 0 0   , R 3 =   0 0 1 √ 2/2 - √ 2/2 0 √ 2/2 √ 2/2 0   .
If W V is isotropic, then the Piola stress tensor should transform as

Π(Rξ) = RΠ(ξ).
For A 1 and A 2 , and each of the three rotations, the diagonal entries of T R i Π(R i A j ) are reported on in Tables 2 and3. The difference with Π(A j ) is a measure of the anisotropy.

As can be seen, the difference on the diagonal terms is of order 0.1%, which shows the network can be considered isotropic, and that the energy density W V is indeed isotropic.

5.4. Dependence on the second invariant. It is rather well admitted that the energy density of rubber materials slightly depends on the second invariant I 2 of the Cauchy-Green strain tensor. By definition, the Arruda-Boyce model only depends on the first and third invariants. As would show an elementary series expansion, the Treloar model actually depends on the second invariant -although the coefficient in front of I 2 is actually rather small compared to the coefficient for I 1 . This is also the case for the variational model as we shall show now. Note that Treloar's experiments do not allow to verify whether a model depends on I 2 or not.

We write the energy as a function of the invariants I 1 , I 2 , and I 3 , so that

W V (ξ) = ω(I 1 , I 2 , I 3 )
for some function ω : R 3 + → R + . This implies for a deformation of the form ξ = diag(ξ 1 , ξ 2 , ξ 3 ) that the Piola stress tensor is diagonal, and that its entries are given for all i ∈ {1, 2, 3} by

Π(ξ) ii = 2ξ i ∂ω ∂I 1 (I 1 , I 2 , I 3 ) + 2ξ i (I 1 -ξ 2 i ) ∂ω ∂I 2 (I 1 , I 2 , I 3 ) + 2I 3 ξ i ∂ω ∂I 3 (I 1 , I 2 , I 3 ).
Let us fix some i 1 > 0 and i 3 = 1, and set α(

I 2 ) = 2 ∂ω ∂I 1 (i 1 , I 2 , i 3 ), β(I 2 ) = 2 ∂ω ∂I 2 (i 1 , I 2 , i 3 ), and γ(I 2 ) = 2 ∂ω ∂I 3 (i 1 , I 2 , i 3 ).
Then, for all ξ = diag(ξ 1 , ξ 2 , ξ 3 ) such that I 1 = i 1 and I 3 = i 3 , the Piola stress tensor is given by

Π(ξ) = α(I 2 )ξ + β(I 2 )diag(ξ 1 (i 1 -ξ 2 1 ), ξ 2 (i 1 -ξ 2 2 ), ξ 3 (i 1 -ξ 2 3 )) + γ(I 2 )diag(ξ -1 1 , ξ -1 2 , ξ -1 3 ). ( 21
)
In particular, W does not depend on I 2 if α and γ are constant functions, and if in addition β ≡ 0.

The physical parameters for these tests are as follows:

N = 26.5, K = |B(R)| N R × 0.27 MPa, K = 100 MPa.
We have conducted numerical tests with i 1 = 9 and i 3 = 1 for the following three deformation gradients ξ: These deformation gradients ξ 1 , ξ 2 , ξ 3 satisfy I 2 = 9, 18.333333, 48.25, respectively. For each deformation gradient, we have computed α, β, and γ by solving a 3 × 3 linear system given by [START_REF] Heinrich | Theoretical and numerical formulation of a molecular based constitutive tube-model of rubber elasticity[END_REF]. Caution should be taken since the systematic and random errors in the computations may be of the order of the coefficients β and γ. We have therefore made explicit a confidence interval, which corresponds to the variance of the coefficients for several realizations (of the polymer network). Tables 4 and5 report on the values of α, β, and γ and of their standard deviations for each deformation gradient (with 100 realizations), for the Treloar and the variational model, respectively. Note that the results on the Treloar model are obtained numerically using the same meshes as for the variational model and within the affine assumption, so that there is also statistical variability due to the randomness of the meshes.

ξ 1 = diag(1,
The dependence of the Treloar model upon the second invariant is clear for ξ 1 , ξ 2 and ξ 3 since the random error on β is of order of 10%. For the variational model, the order of magnitude of β and its sign seem reasonable. Yet, the error is of the order of 50% for the cases ξ 1 and ξ 2 . For the case ξ 3 , the error is only of the order of 10%, which clearly shows that the variational model also slightly depends on the second invariant.

5.5. Validity of the affine assumption. The affine assumption states that the deformation of the sample B R which minimizes the energy for the Dirichlet boundary conditions x → ξ •x is close enough to the linear deformation x → ξ •x in B R so that the minimizer can be replaced by the the linear deformation for the computation of W V (ξ) through [START_REF] Gloria | Random parking, Euclidean functionals, and rubber elasticity[END_REF], and of Π(ξ) (the associated Piola stress tensor). The affine assumption is largely questioned in the litterature, [START_REF] Flory | Network topology and the theory of rubber elasticity[END_REF]. In what follows, we denote by Π(ξ) the Piola stress tensor computed with the linear deformation (that is, within the affine assumption). The numerical tests α ∆α β ∆β γ ∆γ ξ 1 0.3151 ± 0.0015 -0.00342 ± 0.00022 0.00012 ± 0.00017 ξ 2 0.3144 ± 0.0020 -0.00334 ± 0.00041 0.000043 ± 0.00012 ξ 3 0.3290 ± 0.0015 -0.00372 ± 0.00028 0.000040 ± 0.000037 Table 4. Set of parameters for the Treloar model α ∆α β ∆β γ ∆γ ξ 1 0.2898 ± 0.00075 -0.00040 ± 0.00017 0.02177 ± 0.00032 ξ 2 0.2888 ± 0.0015 -0.00057 ± 0.00034 0.01833 ± 0.00021 ξ 3 0.3063 ± 0.00040 -0.00131 ± 0.00013 0.042400 ± 0.00058 Table 5. Set of parameters for the variational model are performed using the natural state of W V as the reference configuration, so that the identity has least energy.

Our numerical experiments show there are two distinct regimes: at moderate deformation, the affine assumption seems to be justified, whereas at large deformation, nonlinear effects play an important role. The threshold above which nonlinear effects are important depends on two parameters:

• the relative strength K of the volumetric term W Helm (the higher, the closer to the affine assumption), • the connectivity of the network (the higher, the more rigid).

To illustrate this, we have conducted numerical simulations in uniaxial traction for the variational model with the following sets of parameters:

N = 26.5, K = |B(R)| N R × 0.27 MPa, K ∈ {5, 50, 500} MPa.
The connectivity of the network influences the parameter |B(R)| N R : for instance, for a fixed random point set in B(R) the number of chains N R is 5 times smaller if the connectivity is set from 20 to 4.

More precisely, the deformation gradient ξ is taken of the form ξ = diag(λ, λ -1/2 , λ -1/2 ) for some λ ≥ 1, and we have plotted the so-called the so-called engineering stress given by

S uniaxial = σ 11 -σ 22 λ .
As announced, the parameter K controls the threshold of non-affinity, as can be seen on Figure 10 which displays the engineering stress for K ∈ {5, 50, 500} MPa and connectivity 20. On Figure 11, we have plotted the engineering stress for two connectivities and K = 50 MPa: the original network obtained from the parking lattice (connectivity around 20), and a modification of this network where edges have been randomly deleted to obtain a connectivity around 4. This illustrates the fact that the connectivity controls the threshold of non-affinity as well. To complete this study, we have plotted on Figure 12 the engineering stress for K ∈ {5, 50, 500} MPa and a network with connectivity 4. Let us conclude this discussion by looking more closely at this nonlinear effect. We turn to the two matrices A 1 and A 2 considered in Section 5. For N R = 59500, n R = 10, the parking lattice (with connectivity 20), and K = 100 MPa (that is of the order of 400 times the characteristic energy of the polymer chains), the first Piola stress tensors are denoted by Π for the Treloar model, and by Π for the variational model. They are given by for A 1 For A 1 , the difference on the diagonal terms is quite small, 0.7%, 1%, and 2%, respectively. This is completely negligible in comparison to the convergence error in N R . Yet, for A 2 , the difference on the diagonal terms can be quite important and cannot be neglected: 0.4%, 5%, and 79%, respectively.

In this subsection we have put in evidence the interplay between non-affinity, the connectivity of the network, and the volumetric term in the energy. Physics should in principle guide us for the connectivity of the network and the volumetric term (recall Remark 2). Note that this interplay is beyond the scope of both the Treloar model (since this is a nonlinear effect due to the non-affinity) and of the Arruda-Boyce model (since the network is reduced there to eight chains -there is no connectivity at all).

5.6. The Rivlin effect. In the experiment literature, models are usually compared using the so-called Mooney plot, which is sensitive to relative errors. This plot only makes sense when the identity is a stress-free configuration. For uniaxial experiments, whose deformation imposed at the boundary is typically given by the linear map

ξ = diag(λ, λ -1/2 , λ -1/2 )
for λ > 0, the Mooney plot consists in dividing the engineering stress

S uniaxial = σ 11 -σ 22 λ
by the universal geometrical factor 2(λ -1/λ 2 ), plotted against 1/λ. In particular, a material exhibits the so-called Rivlin effect if the curve of the Mooney stress

1 λ → σ 11 (λ) -σ 22 (λ) 2(λ 2 -1/λ)
is strictly concave around λ = 1. Rubber materials generically exhibit such a Rivlin effect, see for instance Figure 9 in [START_REF] Rubinstein | Elasticity of polymer networks[END_REF].

It is therefore of interest to check whether this is the case as well for the discrete variational model for rubber studied in this article.

The physical parameters for these tests are as follows:

N ∈ {26.5, 250}, K = |B(R)| N R × 26.5 N × 0.27 MPa, K ∈ {5, 50, 500} MPa,
and the connectivity is around 20 and 4. The numerical tests are performed using the natural state of W V as the reference configuration, so that the identity has least energy. Numerical tests show that the concavity of this curve is driven by the average number of monomers N per polymer chain, the incompressibility factor K, and the connectivity of the network. For small number of monomers (N = 26.5), there is no Rivlin effect (even for small incompressibility factor K = 5 MPa and small connectivity around 4), see Figure 13. For large number of monomers (N = 250), there is a Rivlin effect provided the incompressibility factor is small (K = 5 MPa) or provided the connectivity is small (around 4). A series of tests illustrate these tendencies:

• For N = 250, and connectivity around 20, there is a Rivlin effect for K = 5 MPa, Figure 14 (and not for K = 50 MPa); • For N = 250, and connectivity around 4, there is a strong Rivlin effect for K = 5 MPa, a disctinctive Rivlin effect for K = 50 MPa, and no Rivlin effect for K = 500 MPa, see Figures 15,[START_REF] Gloria | Numerical approximation of effective coefficients in stochastic homogenization of discrete elliptic equations[END_REF]and 17. This dependence upon N is due to the form of the inverse Langevin function, which somehow forbids relaxation for small N . Note that the Arruda-Boyce model never yields the Rivlin effect, as well as the Treloar model.

6.

Comparison of the three models 6.1. General relations for the energy. Although it is not clear a priori that the Arruda-Boyce model is related to the Treloar model, a close look at the defining formulas of both energy densities shows that the Arruda-Boyce model is an elementary relaxation of the Treloar model using convexity. In particular, both the variational and the Arruda-Boyce models are relaxations of the Treloar model: For all deformation gradients ξ,

W V (ξ) W AB (ξ) ≤ W T (ξ). (22)
The first inequality holds because the linear deformation is a test function for the minimization problem defining W V . For the second inequality, we recall that for all N , t → W p c (t, N ) is convex on R + so that Jensen's inequality yields: Although these general relations only hold for the energy density, this suggests that the Arruda-Boyce and variational models are softer than the Treloar model. We thus expect principal stresses to be higher for the Treloar model. This is actually what we observe in the following subsection.

W p T (ξ) = 1 4π π 0 2π 0 W p c (λ 2 1 cos 2 φ + λ 2 2 sin 2 φ cos 2 θ + λ 2 3 sin 2 φ sin 2 θ) sin φdθdφ = 1 4π 1 -1 2π 0 W p c (λ 2 1 s 2 + λ 2 2 (1 -s 2 ) cos 2 θ + λ 2 3 (1 -s 2 ) sin 2 θ)dθds ≥ W p c 1 4π 1 -1 2π 0 (λ 2 1 s 2 + λ 2 2 (1 -s 2 ) cos 2 θ + λ 2 3 (1 -s 2 ) sin 2 θ)dθds = W p c λ 2 1 + λ 2 2 + λ 2 3 3 = W p AB (ξ), where W p c is such that for N fixed, W p c (r 2 /(N l 2 )) = Wc (r, N ) for all r ≥ 0.

Comparison to Treloar's experiments.

In this subsection we compare the three models to Treloar's experiments. The plots 18, 19, and 20 display the engineering stresses (or nominal stresses) associated with the mechanical experiments in uniaxial traction, biaxial traction, and in planar tension, respectively. The plots of the Treloar and Arruda-Boyce models are not new, but the plot of the variational model is. The numerical tests are performed using the natural state of W V as the reference configuration, so that the identity has least energy. Let us quickly recall the definitions of the engineering stress for these three types of sollicitation. Uniaxial traction corresponds to a deformation gradient ξ of the form ξ = diag(λ, λ -1/2 , λ -1/2 ) for some λ ≥ 1. The quantity λ-1 is called the engineering strain. The associated Cauchy stress tensor σ is diagonal, and the associated engineering stress is given by λ -1 → S uniaxial = σ 11 -σ 22 λ .

For biaxial traction the deformation gradient ξ is of the form ξ = diag(λ -2 , λ, λ)

for some λ ≥ 1. The quantity λ -2 -1 is called the engineering strain. The associated Cauchy stress tensor σ is diagonal, and the associated engineering stress is given by λ -2 -1 → S biaxial = σ 11 -σ 22 λ -2 . Finally, the planar tension experiment corresponds to a deformation gradient ξ of the form ξ = diag(1, λ, λ -1 ) for some λ ≥ 1. The quantity λ-1 is called the engineering strain. The associated Cauchy stress tensor σ is diagonal, and the associated engineering stress is given by λ -1 → S planar = σ 22 -σ 11 λ .

For the Arruda-Boyce model, we have used standard coefficients to fit Treloar's data, namely N = 26.5 and n β = 0.27 MPa. For the numerical simulation of the variational model, the stochastic network has been generated by the random parking lattice. We have not tried to optimize the parameters further, and have taken Note that the choice K = 5 MPa is not physical, and this should be optimized together with the network, see Remark 2. As can already be seen, the counterpart to decrease the uniaxial stress is to increase the transversal stresses. This tendency seems to be more pronounced for the variational model. This is even clearer in larger deformation. For the strain gradient A = diag(6.5052, 0.3920755, 0.3920755) Note that the conclusions strongly depend on the strength K of the volumetric energy density and on the connectivity of the network regarding the variational model.

Perspectives

In this article we have compared three discrete models for rubber: the Treloar model, the Arruda-Boyce model, and a variational model based on the statistical physics of reticulated polymer chains. In particular we have seen that both the Arruda-Boyce and the variational models are relaxations of the Treloar model. The three models yield nice mathematical and mechanical properties, and compare rather well to Treloar's experiments. The comparison of the variational model to Treloar's experiments is indeed promising, although the choice of the volumetric stiffness K to obtain these results is not physical. There is definitely room for improvement of the model using more physical input (essentially more information on the polymer chain network), and taking inspiration from more advanced models such as the non-affine microsphere model, the maximal advance path constraint model, the nonaffine tube model, or the constraint junction models. The following program should be addressed to go further:

• generate more realistic networks in terms of connectivity and molecular mass (we expect a better relaxation at moderate strain and still an increase of the stiffness at large strain);

• find explicit formulas for the constitutive law W V by solving an inverse problem in order to turn the model into an attractive alternative to the engineering community; • compare the simulations to physical experiments at the scale of the polymer-chain network (this would challenge the validity of the discrete model itself). These three issues will be addressed in future works.

Figure 1 .

 1 Figure 1. Arruda-Boyce model: undeformed and deformed representative element

Figure 2 .

 2 Figure 2. Variance of the diagonal terms of the Piola stress tensor for A 1 (range [2, 5] × [-3.4, -1.8])

Figure 3 .

 3 Figure 3. Variance of the diagonal terms of the Piola stress tensor for A 2 (range [2, 5] × [-3, -1])

11 Figure 4 . 22 Figure 5 .Figure 6 . 11 Figure 7 . 22 Figure 8 .Figure 9 .

 11422561172289 Figure 4. Convergence of the diagonal part 1 of the Piola stress tensor for A 1 (range [0, 6 10 4 ]×[0.42, 0.46])

Figure 10 .Figure 11 .Figure 12 .

 101112 Figure 10. Uniaxial traction -affine assumption (dashed line), variational model for K = 500, 50, 5 MPa and connectivity 20 (from top to bottom)

Figure 13 .

 13 Figure 13. Mooney plot 1/λ → σ 11 -σ 22 2(λ 2 -1/λ) for N = 26.5, K = 5 MPa and connectivity 20 (logarithmic scale). The curve is convex, indicating that there is no Rivlin effect (range [10 -1 , 15] × [0.5, 0.7]).

Figure 14 .

 14 Figure 14. Mooney plot 1/λ → σ 11 -σ 22 2(λ 2 -1/λ) for N = 250, K = 5 MPa and connectivity 20 (logarithmic scale). The curve is locally concave, indicating that there is a distinct Rivlin effect (range [10 -1 , 15] × [0.055, 0.056]).

Figure 15 .

 15 Figure 15. Mooney plot 1/λ → σ 11 -σ 22 2(λ 2 -1/λ) for N = 250, K = 5 MPa and connectivity 4 (logarithmic scale). The curve is concave, and there is a strong Rivlin effect (range [10 -1 , 15] × [0.045, 0.053]).

Figure 16 .

 16 Figure 16. Mooney plot 1/λ → σ 11 -σ 22 2(λ 2 -1/λ) for N = 250, K = 50 MPa and connectivity 4 (logarithmic scale). The curve is locally concave, and there is a distinct Rivlin effect (range [10 -1 , 15] × [0.054, 0.056]).

Figure 17 .

 17 Figure 17. Mooney plot 1/λ → σ 11 -σ 22 2(λ 2 -1/λ) for N = 250, K = 500 MPa and connectivity 4 (logarithmic scale). The curve is convex, and there is no Rivlin effect (range [10 -1 , 15] × [0.056, 0.057]).

N

  = 26.5, K = |B(R)| N R × 0.27 MPa, K = 5 MPa.

6. 3 .

 3 Comparison of the Cauchy stress tensors for general deformations. Although the comparison to Treloar's experiments seem to reveal that the Arruda-Boyce and the variational models are very close to each other, a closer look at the results shows there are significant differences. In particular, we have only compared the engineering stress, and not the full Cauchy stress tensor (which is unavailable in Treloar's experiments). Let us go back to the case of uniaxial extension, with the parameters N = 26.5, K = |B(R)| N R × 0.27 MPa, K = 5 MPa.

Figure 18 .

 18 Figure 18. Uniaxial traction -Treloar's experiments

2 Figure 19 .

 219 Figure 19. Biaxial tension -Treloar's experiments

Figure 20 .

 20 Figure 20. Planar tension-Treloar's experiments

Table 2 .

 2 Diagonal entries for the rotations of A 1

	Id 0.8201 1.193 0.3246
	R 1 0.8209 1.197 0.3257 R 2 0.8233 1.192 0.3256 R 3 0.8215 1.192 0.3259 Table

Acknowledgments

We gratefully acknowledge the help of Eric Saltel and Paul-Louis Georges from the projectteam GAMMA of INRIA for providing us with a customized Delaunay tessellation software, and the stimulating discussions on the subject with François Lequeux and Michael Rubinstein. We also thank the anonymous referees for their very insightful comments and suggestions on the manuscript. The research of the authors was supported by INRIA, under the grant "Action de Recherche Collaborative" DISCO.


Dhanjal Charanpal
email: charanpal.dhanjal@telecom-paristech.fr

Ltci

Romaric Gaudel
email: romaric.gaudel@univ-lille3.fr

Stéphan Clémençon
email: clemenco@telecom-paristech.fr

Incremental Spectral Clustering with the Normalised Laplacian

. Beyond a precise description of its practical implementation and an evaluation of its complexity, its performance is illustrated through numerical experiments, based on datasets modelling the evolution of a HIV epidemic and the purchase history graph of an e-commerce website.

Introduction

Graph-mining has recently received increasing attention in the machine-learning literature, motivated by application domains such as the Internet, social networks, epidemics of transmissible infectious diseases, sensor and biological networks. A crucial task in exploratory analysis and data visualization is graph clustering [START_REF] Schaeffer | Graph clustering[END_REF][START_REF] Fortunato | Community detection in graphs[END_REF], which aims to partition the vertices in a graph into groups or clusters, with dense internal connections and few connections between each other. There is a large body of work on graph clustering. A possible approach is to consider a certain measure that quantifies community structure and formulate the clustering issue as an optimization problem (which is generally NP-hard), for which fairly good solutions can be obtained recursively or by using ade-1 quate metaheuristics, see [START_REF] Newman | Finding and evaluating community structure in networks[END_REF][START_REF] Newman | Fast algorithm for detecting community structure in networks[END_REF][START_REF] Flake | Graph clustering and minimum cut trees[END_REF][START_REF] Satuluri | Scalable graph clustering using stochastic flows: applications to community discovery[END_REF] for example. In this paper, focus is on the spectral clustering approach [START_REF] Luxburg | A tutorial on spectral clustering[END_REF], which performs well empirically, is often simple to implement and benefits computationally from fast linear algebra libraries. The general principle of spectral clustering is to compute the smallest eigenvectors of some particular matrix L (refer to Section 2 for further details) and then cluster the vertices based on their representation in the corresponding eigen-space. The popularity of this approach arises from the fact that the obtained clustering of vertices is closely connected to the spectral relaxation of the minimization of the normalized cut criterion, see [START_REF] Shi | Normalized cuts and image segmentation[END_REF].

In many applications such as communications networks (e.g. the Web and Internet), biological networks (of proteins, metabolic reactions, etc.), social networks or networks of scientific citations for instance, the graphs of interest slowly change with time. A naive approach to this incremental problem is to cluster each graph in the sequence separately, however this is computationally expensive for spectral clustering as the cost of solving the eigenvalue problem is O(n 3 ) at each iteration, where n is the number of vertices. There has been some previous work on the incremental spectral clustering problem, for example [START_REF] Valgren | Incremental spectral clustering and its application to topological mapping[END_REF][START_REF] Ning | Incremental spectral clustering by efficiently updating the eigen-system[END_REF][START_REF] Ning | Incremental spectral clustering with application to monitoring of evolving blog communities[END_REF] however only [START_REF] Ning | Incremental spectral clustering by efficiently updating the eigen-system[END_REF][START_REF] Ning | Incremental spectral clustering with application to monitoring of evolving blog communities[END_REF] update the eigen-system. In this paper we propose an efficient method for clustering a sequence of graphs which leverages the eigen-decomposition and the clustering on the previous graph to find the new clustering. Firstly, a fast approximation of a rank-k eigen-decomposition of L t+1 knowing that of L t is derived from the Singular Value Decomposition (SVD) updating approach used for Latent Semantic Indexing in [START_REF] Zha | On updating problems in latent semantic indexing[END_REF]. Here, the update is efficient to compute when the change (that will be described more precisely in the subsequent analysis) between L t and L t+1 is a low-rank matrix. Secondly the clustering of vertices is updated accordingly to the new eigen-space. The efficiency of the complete approach, in term of clustering accuracy, is demonstrated on real data.

The paper is organised as follows. A Laplacian-based graph clustering approach is recalled in Section 2. Then Section 3 details the proposed eigen-decomposition update and Section 4 describes its use in the framework of spectral clustering. Numerical results are gathered in Section 5, and the paper ends with Section 6 discussing results and planned future work.

Notation: A bold uppercase letter represents a matrix, e.g. X, and a column vector is displayed using a bold lowercase letter, e. g. x. The transpose of a matrix or vector is written X T . The concatenation of the columns of two matrices A and B is written

[A B].

Laplacian Based Graph Clustering

Consider an undirected graph G = (V, E), composed of a set of vertices {v 1 , . . . , v n } = V and edges E ⊆ V × V such that for every edge (v i , v j ) there is also an edge (v j , v i ). One way of representing the edges of G is using an adjacency matrix A which has A ij = 1 if there is an edge from the i-th to the j-th vertex otherwise A ij = 0. More generally, the weight matrix W allows one to assign nonzero numerical values on edges, and thus W ij = 0 when there is an edge from the i-th to the j-th vertex.

In the perspective of spectral clustering, a useful way of representing G is through its Laplacian matrix [START_REF] Chung | Spectral graph theory[END_REF]. There are several definitions of the Laplacian matrix, however we are interested in the normalised Laplacian matrix, which is symmetric and positive semi-definite:

Definition 2.1. The unnormalised Laplacian matrix of a graph G is defined as L = D -W where D is the degree matrix with diagonal entries D ii = deg(v i ), denoting by deg(v i ) = j W ij
the degree of the i-th vertex, and zeros elsewhere, and W is the weight matrix. The normalised Laplacian matrix of a graph G is then defined as

L = D -1 2 LD -1 2 .
The normalised Laplacian matrix is used in the spectral clustering approach of Ng et al. [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF].

The algorithm computes the Laplacian and then finds the k smallest eigenvectors which form the columns of a matrix V. This matrix is then normalised so that its rows have unit norm, which are then clustered with the k-means algorithm. In general the Laplacian matrix is often sparse and one can use power or Krylov subspace methods such as the Lanczos method to find the eigenvectors.

There are several variants of this algorithm, such as that of [START_REF] Shi | Normalized cuts and image segmentation[END_REF] which uses the so-called random walk Laplacian and clusters the smallest eigenvectors in a similar way. One of the motivations for these clustering methods is from the spectral relaxation of the minimisation of the normalised cut criterion.

A naive approach to spectral clustering on a sequence of graphs has a large update cost due to the computation of the eigen-decomposition of Lt at each iteration t. A more efficient approach is to update relevant eigenvectors from one iteration to the next. In [START_REF] Valgren | Incremental spectral clustering and its application to topological mapping[END_REF] the authors use the spectral clustering of [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] however they do not update the eigenvectors incrementally but instead the clustering directly.

The iterative clustering approach of Ning et al. [START_REF] Ning | Incremental spectral clustering by efficiently updating the eigen-system[END_REF][START_REF] Ning | Incremental spectral clustering with application to monitoring of evolving blog communities[END_REF] uses the spectral clustering method given in [START_REF] Shi | Normalized cuts and image segmentation[END_REF] and updates eigenvectors incrementally. The algorithm incrementally updates the generalised eigenproblem Lv = λDv by finding the derivatives on the eigenvalues/vectors with perturbations in all of the quantities involved. An iterative refinement algorithm is given for the eigenvalues and eigenvectors given a change in the edges or vertices of a graph. One then clusters the resulting k smallest eigenvectors using k-means clustering. In order to limit errors which can build up cumulatively, the authors recompute the eigenvectors after every R-th graph in the sequence.

A disadvantage of the approach of [START_REF] Ning | Incremental spectral clustering by efficiently updating the eigen-system[END_REF][START_REF] Ning | Incremental spectral clustering with application to monitoring of evolving blog communities[END_REF] is that to update an eigenvector one must invert a small matrix for each weight change which makes updates costly. The size of this matrix is proportional to the number of neighbours of the vertices incident to the changed edge.

Incremental Eigen-approximation

In this section we address the updating problem upon the rank-k eigen-decomposition of a symmetric matrix, i.e. we are interested only in the largest k eigenvalues/vectors of the corresponding matrix: 1 This update is a general operation, however it will be explained in the context of spectral clustering later in Section 4. We assume that Y 1 , Y 2 ∈ R n×p and one does not have direct access to A ∈ R m×n but only the matrice A T A.

(A T A) k + U, where U = Y 1 Y T 2 + Y 2 Y T 1 and (A T A) k is the rank-k eigen-decomposition of A T A.
To compute this update, we first project the columns of Y 1 into the space orthogonal to the eigenvectors Q k of (A T A) k , a process known as deflation. Hence the matrix Y 1 is deflated as follows:

Ȳ1 = (I -Q k Q T k )Y 1 . Note that Ȳ1 Θ 1 for some Θ 1 , is orthogonal to Q k since Q T k Ȳ1 Θ 1 = (Q T k Y 1 -Q T k Y 1 )Θ 1 = 0. If we take the SVD Ȳ1 = P1 Σ1 QT 1 then P1 is orthog- onal to Q k since P1 = Ȳ1 Q1 Σ-T 1 assuming Σ1 has nonzero diagonal entries.
At the next stage we would like to orthogonalise the columns of Y 2 with respect to both Q k and P1 . Hence we deflate Y 2 in the following way:

Ȳ2 = (I -P1 PT 1 -Q k Q T k )Y 2
where we have used the fact that P1 is orthogonal to Q k . Proved in a similar way to the step used earlier, the matrix in the column space of Ȳ2 , Ȳ2 Θ 2 for some Θ 2 , is orthogonal to Q k and P1 . Hence, we compute the SVD Ȳ2 = P2 Σ2 QT 2 and note that the matrices P1 , P2 and Q k are mutually orthogonal and span the space spanned by (A T A) k + U. This allows one to write

(A T A) k + U = Q∆ QT with Q = [Q k P1 P2 ] ∈ R n×(k+2p) and ∆ =    Ω k + Q T k UQ k Q T k U P1 Q T k Y 1 Q2 Σ2 PT 1 UQ k PT 1 U P1 Σ1 QT 1 Q2 Σ2 Σ2 QT 2 Y T 1 Q k Σ2 QT 2 Q1 Σ1 0    ,
where, the above expression of ∆ is simply an expansion of the matrix QT ((A T A) k + U) Q since Q is a basis of the space occupied by ((A T A) k + U). We take the rank-k eigen-decomposition

∆ k = H k Π k H T
k and then the final eigen-approximation is given by ( QH k )Π k ( QH k ) T . The deflation and SVD of Ȳ1 and Ȳ2 cost O(npk) and O(np 2 ) respectively and the eigendecomposition of ∆ is O((k + 2p) 3 ) as ∆ ∈ R (k+2p)×(k+2p) . In order to compute ∆ one can reuse the computations of

Q T k Y 1 , Q T k Y 2 , PT 1
Y 2 which are used for deflations and also the matrices used for the SVD decompositions. Thus ∆ is found in O(p 3 + p 2 k + pk 2 ), and the overall complexity of this algorithm is O((k 2 + p 2 )(p + k) + np(p + k)). Of note here is that the dimension n of A T A scales the complexity in a linear fashion, however costs are cubically related to k and p.

Under the specific circumstances of working with sparse matrices and updating graph clustering this cost is later shown to reduce.

Eigenupdates corresponding to the addition and deletion of rows/columns are a particular case of the considered update, in which some rows/columns have zero values before or after the update. However, these are not detailed due to space limitations.

Incremental Cluster Membership

We return to the eigen-problem Lv = λv in which we are interested in the eigenvectors with the smallest eigenvalues. We define the shifted Laplacian as L = 2I -L = I + D -1 2 WD -1 2 , which is positive semi-definite since L is positive semi-definite with largest eigenvalue 2. Since we are interested in the smallest eigenvectors of L they correspond exactly to the maximum eigenvectors of L and we can use the eigen-update methods described above.

Putting the ingredients together allows us to outline an efficient incremental method for performing graph clustering called Incremental Approximate Spectral Clustering (IASC), see Algorithm 1. At a high level the algorithm is quite simple: at each iteration t one mostly recovers the steps of standard spectral clustering. In this case, we use the shifted Laplacian matrix Lt instead of the normalised Laplacian matrix, which does not affect the clustering as both matrices have same eigen-vectors. Furthermore we use the eigen-approximation method above to approximate the largest eigenvectors of Lt using the eigenvectors computed at the previous iteration, instead of (costly) recomputing the exact eigen-decomposition at each iteration. Similarly, we initialise k-means with centers at previous iteration. Finally, we store the first k 2 eigenvectors of Lt , with k 2 ≥ k 1 , to reduce cumulative errors introduced in the loop at the expense of increased computation.

Algorithm 1 Incremental Approximate Spectral Clustering using the Normalised Laplacian Require: Graphs G 1 , . . . , G T of sizes n 1 , . . . , n T , no. of clusters k 1 , matrix approximation rank

k 2 ≥ k 1 , eigen-decomposition recomputation step R 1: for t = 1 → T do 2:
Compute the shifted Laplacian for G t , Lt 3:

if i % R == 1 then 4: Compute rank-k 2 eigen-decomposition of Lt , Q (t) k2 Ω (t) k2 Q (t) k2 T 5: else 6:
Use rank-k 2 eigen-decomposition of Lt-1 and approximated eigen-update of Section 3 to compute rank-k 2 eigen-decomposition of Lt , Q

(t) k2 Ω (t) k2 Q (t) k2 T 7: end if 8: Let Q (t) k1 be the matrix of the first k 1 columns of Q (t) k2 9:
Normalise the rows of

Q (t) k1 , i.e. Q (t) k1 ← diag(Q (t) k1 (Q (t) k1 ) T ) -1 2 Q (t) k1 10:
Use k-means on rows of

Q (t)
k1 , using c t-1 to find initial centroids (if not at first iteration), and store indicators c t ∈ {1, . . . , k 1 } nt 11: end for 12: return Cluster membership vectors c 1 ∈ {1, . . . , k} n1 , . . . , c T ∈ {1, . . . , k} n T The complexity of Algorithm 1 is dictated by the sparsity of the graphs and the extent of the change between successive graphs in the input sequence. At iteration t and step 2 one can compute Lt from the weight and degree matrix at a cost of O(n t + |E t |). In the following step, if there is a change between edges incident to vertices S = {v I1 , . . . , v I ℓ } then the rows and columns corresponding to the union of the neighbours of S, n(S), will change in the corresponding shifted Laplacian. In this case p = |n(S)| in Y 1 , Y 2 ∈ R n×p and if n(S) is small, this update can be efficiently computed.

Computational Results

In this section we compare the clustering quality of: standard spectral clustering (denoted Exact) Ning et al.'s incremental strategy, and IASC. We use the clustering methods on two real datasets. The first one represents individuals in Cuba who have been detected as having HIV between 1986 and 2004, see [START_REF] Auvert | The HIV/AIDS epidemic in Cuba: description and tentative explanation of its low HIV prevalence[END_REF] for precise details of the dataset. Edges in the graph represent a sexual contact between two individuals as determined using contact tracing, which is a technique used in epidemiology to limit the spread of an infection whereby contacts of an infected person are identified and diagnosed. The full sexual contact graph at the end of 2004 consists of 5389 people however we consider the largest connected component of size 2387. We examine the epidemic graph at 3 month steps and cluster using the 3 approaches starting at the point in which the graph has 500 vertices. The experiment is run using k 1 = 25, k 2 = 100 and R = 10 for IASC. For the method of Ning et al. we also use 25 clusters and recompute exact eigenvectors after every 10 iterations. Note that in our implementation of k-means clustering, k is an upper bound on the number of clusters.

The next dataset was introduced in [START_REF] Richard | Link discovery using graph feature tracking[END_REF] and corresponds to the purchase history of an e-commerce website over a period of almost two years. The dataset is a bipartite graph between users and product with more than 700,000 users and 1,200,000 products. In the current experiment we focus on the first 5000 users. A graph is constructed between users with edge weights equal to the number of commonly purchased products between two users.

To evaluate the learned clusters we measure the k-way normalised cut, defined as follows:

N = 1 k k ℓ=1 ij W ij δ(c i , c ℓ )(1 -δ(c j , c ℓ )) ij W ij δ(c i , c ℓ )
.

A cut between two clusters A and B is the sum of the weights between the clusters and the normalised cut is this sum divided by the sum of the weights of all edges incident to vertices in cluster A. Hence the k-way normalised cut is the mean normalised cut between each cluster and its complementary vertices, and the lower it is the better.

Figure 1 shows the k-way normalised cuts for both datasets. In both IASC and the Ning et al. method, since eigenvectors are recomputed every 10 iterations, this manifests itself as sudden changes in the k-way normalised cuts. These changes are more pronounced in the Ning et al. case. Except for the first few iterations, IASC matches or improves results on the exact approach, while keeping only 10% or fewer of the eigenvectors. The Ning et al. approximation leads to cuts worse than those obtained with the exact approach, except for the first few iterations immediately after exact recomputations of the eigenvectors. For the Bemol graphs Ning et al. failed to complete in a reasonable time frame, however on a smaller dataset of 1000 vertices, it resulted in cuts considerably worse than either IASC or the exact spectral approach.

Discussion

We have presented a novel incremental method for spectral graph clustering which updates the eigenvectors of the normalised Laplacian in a computationally efficient way. Such an algorithm is useful for finding clusterings in time evolving graphs such as biological and social networks, and the Internet. We compare our approach empirically to a naive spectral clustering approach which recomputes eigenvectors for each graph, and the method of Ning et al. and demonstrate improvements in the clustering on HIV and product information datasets as well as good cluster stability.

Recall that the approximation error of our proposed approach is kept under control by recomputing the exact eigen-decomposition after every R iterations. One promising extension to this work is to estimate this approximation error. This would would allow one to limit the frequency of this computationally expensive step and offer a way to detect abrupt changes in the evolution of a graph sequence.

Figure 1 :

 1 Figure 1: The k-way normalised cuts using the HIV and e-commerce graphs.

We have written the update in terms of a symmetric matrix A T A to improve notation. Note that any positive semi-definite symmetric matrix M can be decomposed into the form M = A T A, where A is real, for example by using an eigen-decomposition.