E Tikarrouchine

G Chatzigeorgiou

F Praud

B Piotrowski

Y Chemisky

F Meraghni
email: fodil.meraghni@ensam.eu

Three-dimensional FE 2 method for the simulation of non-linear, rate-dependent response of composite structures

Keywords: Multi-scale finite element computation, FE 2 method, periodic homogenization, composite materials, elastoviscoplastic behavior, ductile damage

In this paper, a two scale Finite Element method (FE 2 ), is presented to predict the non-linear macroscopic response of 3D composite structures with periodic microstructure that exhibit a timedependent response. The sensitivity to the strain rate requires an homogenization scheme to bridge the scales between the macroscopic boundary conditions applied and the local evaluation of the strain rate. In the present work, the effective response of composite materials where the matrix has a local elasto-viscoplastic behavior with ductile damage are analyzed using periodic homogenization, solving simultaneously finite element problems at the microscopic scale (unit cell) and at the macroscopic scale. This approach can integrate any kind of periodic microstructure with any type of non-linear behavior for the constituents (without the consideration of non-linear geometric effects), allowing to treat complex mechanisms that can occur in every phase and at their interface. The numerical implementation of this simulation strategy has been performed with a parallel computational technique in ABAQUS/Standard,with the implementation of a set of dedicated scripts. The homogenization process is performed using a user-defined constitutive law that solve a set full-field non-linear simulations of a Unit Cell and perform the necessary homogenization of the mechanical quantities. The effectiveness of the method is demonstrated with three examples of 3D composite structures with plastic or viscoplastic and ductile damage matrix. In the first example, the numerical results obtained by this full field approach are compared with a semi-analytical solution on elastoplastic multilayer composite structure. The second example investigates the macroscopic response of a complex viscoplastic composite structure with ductile damage and is compared with the mean field Mori-Tanaka method. Finally, 3D corner structure consisting of periodically aligned short fibres composite is analysed under complex loading path. These numerical simulations illustrate the capabilities of the FE 2 strategy under non-linear regime.

when time dependent constitutive models describe the response of the constituents

Introduction

Polymer based composite materials are considered to be a good technological solution for automotive and aeronautic industries, thanks to their structural durability and their lightness. A major preoccupation of these industries is to predict the response of such structures with in-service loadings. This requires the development of predictive models that are able to capture the microstructure impact on the mechanical response, and the proper identification of the mechanical properties of the constituents. In this purpose, advanced modelling and simulation methods that integrate the effect of the microstructure is an active area of research. According to the bibliography, several numerical approaches have been proposed for the numerical simulation of the non-linear response of polymer based composite structures including: i) Phenomenological models, which predict the overall response of the composite materials without taking into account the effect of the different constituents observed at the microscopic scale. Several authors have proposed constitutive models that integrate various rheologies and deformation mechanisms, i.e. viscoelasticity [START_REF] Moreau | Analysis of thermoelastic effects accompanying the deformation of pmma and pc polymers[END_REF][START_REF] Akhtar | Thermo-mechanical large deformation response and constitutive modeling of viscoelastic polymers over a wide range of strain rates and temperatures[END_REF], viscoplasticity [START_REF] Duan | A uniform phenomenological constitutive model for glassy and semi crystalline polymers[END_REF][START_REF] Achour | Implicit implementation and consistent tangent modulus of a viscoplastic model for polymers[END_REF][START_REF] Drozdov | Cyclic viscoplasticity of solid polymers: The effects of strain rate and amplitude of deformation[END_REF], coupled viscoelasticity and viscoplasticity [START_REF] Miled | Coupled viscoelastic-viscoplastic modeling of homogeneous and isotropic polymers: Numerical algorithm and analytical solutions[END_REF][START_REF] Miled | Micromechanical modeling of coupled viscoelastic-viscoplastic composites based on an incrementally affine formulation[END_REF], or even both coupled viscoelasticity, viscoplasticity and damage [START_REF] Launay | Cyclic behaviour of short glass fibre reinforced polyamide: Experimental study and constitutive equations[END_REF][START_REF] Launay | Multiaxial fatigue models for short glass fiber reinforced polyamide -part i: Nonlinear anisotropic constitutive behavior for cyclic response[END_REF][START_REF] Krairi | A thermodynamically-based constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage[END_REF][START_REF] Praud | Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF];

ii) Multi-scale methods, that can be classified into two main categories: mean-field and full field approaches. The mean-field approaches are used to describe the behavior of composites for certain categories of microstructures through the Mori-Tanaka Method [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Doghri | Homogenization of two-phase elasto-plastic composite materials and structures: Study of tangent operators, cyclic plasticity and numerical algorithms[END_REF] or the self-consistent scheme [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF]15,[START_REF] Walpole | On bounds for the overall elastic moduli of inhomogeneous systems-i[END_REF][START_REF] Milton | Variational bounds on the effective moduli of anisotropic composites[END_REF]. These methodologies have been developed in order to estimate the overall behavior of the composite using average stress and strain quantities for each material phase [START_REF] Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF][START_REF] Meraghni | Micromechanical modelling of matrix degradation in randomly oriented discontinuous-fibre composites[END_REF][START_REF] Meraghni | Implementation of a constitutive micromechanical model for damage analysis in glass mat reinforced composite structures[END_REF]. These methods have been proved to be accurate for the linear cases. However, for non-linear constitutive laws, especially when the matrix phase exhibits a non-linear behavior, the response of these approaches is inaccurate. It is commonly observed in the literature that the response of the composite obtained by mean-field methods appears to be stiffer than the reality especially when the matrix is ductile and the reinforcements are stiffer [START_REF] Gavazzi | On the numerical evaluation of eshelby's tensor and its application to elastoplastic fibrous composites[END_REF][START_REF] Lagoudas | Elastoplastic behavior of metal matrix composites based on incremental plasticity and the mori-tanaka averaging scheme[END_REF][START_REF] Chaboche | On the capabilities of mean-field approaches for the description of plasticity in metal matrix composites[END_REF]. The numerical simulation of these composite systems has necessitated the development of full-field approaches. To determine the response of a composite structure, accounting for the description of the microstructure, the so-called FE 2 method, appear to be an adequate solution. The major benefit of the FE 2 method is the ability to analyse complex mechanical problems with heterogeneous phases that present a variety of behavior at different scales. This idea was originally introduced by Feyel [START_REF] Frédéric | Multiscale fe2 elastoviscoplastic analysis of composite structures[END_REF], then this method was used and developed by several authors, for example [START_REF] Feyel | Fe2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials[END_REF][START_REF] Nezamabadi | A multilevel computational strategy for handling microscopic and macroscopic instabilities[END_REF][START_REF] Nezamabadi | A multiscale finite element approach for buckling analysis of elastoplastic long fiber composites[END_REF][START_REF] Asada | Fully implicit formulation of elastoplastic homogenization problem for two-scale analysis[END_REF][START_REF] Tchalla | An abaqus toolbox for multiscale finite element compu-tation[END_REF][START_REF] Schröder | Algorithmic two-scale transition for magneto-electro-mechanically coupled problems: Fe2-scheme: Localization and homogenization[END_REF][START_REF] Papadopoulos | The impact of interfacial properties on the macroscopic performance of carbon nanotube composites. a fe2-based multiscale study[END_REF].

The majority of these works consider two-dimensional structures, which if they provide a good study case for the analysis of the capabilities of the method, is of limited interest for practical use for the prediction of the overall response of heterogeneous materials and composites, since the spatial arrangement of the phases is mostly three-dimensional.

In this paper, a two-level FE 2 method, based on the concept of periodic homogenization under the small strain assumption is implemented in a commercial FE code (ABAQUS/Standard). The method predicts the 3D non-linear macroscopic behavior of a composite with periodic microstructure by considering that each macroscopic integration point is a material point where the characteristics at the macroscopic scale are represented by its own unit cell, which includes the material and geometrical characteristics of the constituents (fibre, matrix) in the microstructure. Therefore, a multilevel finite element analysis has been developed using an implicit resolution scheme, with the use of a Newton-Raphson algorithm to solve simultaneously the non-linear system of equations on the two scales (macroscopic and microscopic).

The main advantage of this methodology is that it can account for any type of non-linear behavior of the constituents (plasticity, viscoelasticity, viscoplasticity and damage), as well as any type of periodic microstructure. The proposed FE 2 approach is implemented through a parallelization technique, leading to a significant reduction of the computational time.

The layout of this paper is as follows: in section 2, the theoretical formulation of the homogenization theory is described as well as the principle of scale transition between the local and the global fields. The section also presents the rate dependent constitutive law considered for the matrix phase. In section 3, details of the numerical implementation of the FE 2 method is given for a 3D non-linear problem in ABAQUS/Standard with the parallel implementation. In section 4, the approach is validated by comparing the FE 2 results with semi-analytical method on 3D multilayer composite structure. Afterwards, an example of 3D composite structure exhibiting non-uniform strain fields, in which the microstructure consists of an elastoviscoplastic polymer matrix with ductile damage, reinforced by short glass fibres is presented. The numerical results of the simulation are compared with the Mori-Tanaka method. Finally, the capabilities of this method are shown by simulating the mechanical response of a more complex structure under complex loading path with different strain rate.

Theoretical background and Scale transition

In this section, the periodic homogenization principle, as well as the transition between the two scales (microscopic and macroscopic) are presented. The principal objective is to determine the macroscopic quantities (stress and tangent modulus) that are obtained through periodic homogenization by accounting for the different mechanisms that exist in the microscopic level, as non-linear plastic/viscoplastic behavior with ductile damage of the matrix. After that, the local constitutive law of each constituent is presented, where a linear elastic law is chosen for the reinforcement and a constitutive model that incorporate elastoviscoplasticity coupled with ductile damage for the matrix.

Theoretical background for periodic homogenization

The objective of the periodic homogenization theory is to define a fictitious homogenized medium having an equivalent response of the heterogeneous medium that is representative of the microstructure. A periodic medium is characterized by a repeated unit cell in the three spatial directions, which forms an unit cell. The theory of periodic homogenization is valid as long as the separation between the scales exists, i.e. the sizes of the unit cell are much smaller than the macroscopic sizes of the medium (x >> x) (Fig. 1). In this paper, the notation (•) will be used to denote macroscopic quantities. The motion of any macroscopic and microscopic material points M(x) and M(x, x) respectively, are governed by the macroscopic and the microscopic equations (Tab.1).

In Tab.1 σ, ε, σ and ε represent the microscopic and the macroscopic stress and strain tensors respectively, b v is the body forces, V and V are the volumes of the micro and the macro structure.

Equations

Macro-scale Micro-scale

∀ x ∈ V ∀ x ∈ V, ∀ x ∈ V Equilibrium div x σ + b v = 0 div x (σ) = 0 Kinematics ε = 1 2 (Grad x (u) + Grad x (u)) ε = 1 2 (Grad x (u) + Grad x (u)) Constitutive law σ = F(x, ε) σ = F(x, x, ε)
Strain energy rate Ẇε = σ : ε Ẇε = σ : ε Moreover, x, x, u and u are the microscopic and the macroscopic positions and displacement vectors, while F and F are operators that define the micro and macro relationships between the stress and strain. Both F and F are considered non-linear operators in this work.

The homogenization theory attempts to define the F operator, which characterizes the macroscopic behavior, from the local behaviors defined by the F operator. In order to make this possible, it is necessary to introduce the concept of scale transition between the macro and the micro scales.

According to the average stress and strain theorems, it can be demonstrated that the stress and strain averages within the unit cell are equal to the stress and strains corresponding to uniform tractions and linear displacements respectively that are applied at its boundaries. These averages represent the macroscopic stress and strain tensors respectively. The relationships between the two scales are given by the following equations:

σ = σ = 1 V V σ dV = 1 V ∂V σ.n ⊗ x dS (1)
ε = ε = 1 V V ε dV = 1 2V ∂V (u ⊗ n + n ⊗ u) dS ( 2
)
where n is the outgoing normal of the unit cell boundary ∂V. • is the mean operator and ⊗ the dyadic product.

Non-linear scale transition: incremental approach

Since the homogenization is based on the separation between the different scales, the connection between these scales (microscopic and macroscopic problems) should be defined in order to be able to predict the overall behavior of the structure.

Microscopic problem

The periodicity condition implies that, the displacement field u of any material point located in x can be described by an affine part, in which a periodic fluctuation u is added as is presented in Fig. 2: The periodic fluctuating quantity u takes the same value on each pair of opposite parallel sides of the unit cell and the strain average produced by u is null [Eq. 5]. Therefore, the full strain average is well equal to the macroscopic strain [Eq. 6].

u(x, x, t) = ε(x, t) • x + u (x, x, t) (3)
ε (u) = ε + ε (u ) (4)
ε (u ) = 1 V V ε (u ) dV = 0, (5)
ε (u) = ε + ε (u ) = ε (6)
The traction vector σ. n is anti periodic and satisfies the conditions of equilibrium within the unit cell. The micro problem is formulated as follows:

                     σ = F (x, ε (u (x))) ∀x ∈ V, div x (σ (x)) = 0 ∀x ∈ V, u i -u j = ε . x i -x j ∀x ∈ V (7)
where u i , u j , x i and x j are the displacements and the positions of each pair of opposite parallel material point of the unit cell boundary respectively, while ε is the macroscopic strain. The relationship between the microscopic stress and the microscopic strain in incremental approaches is provided by the linearised expression [Eq. 8]:

∆σ(x) = C t (x) : ∆ε (x) ∀x ∈ V, (8)
where C t is the local tangent operator tensor defined as the numerical differentiation of the stress with respect to the total strain.

Macroscopic problem

The relationship between macroscopic stress and strain cannot be explicitly provided by a stiffness tensor. Nevertheless, for a given macroscopic strain, the macroscopic stress response can be computed using an implicit resolution scheme, where the local behavior is linearized and corrected at each strain increment [Eq. 8]. Then, using the same incremental methodology, the macroscopic behavior can also be linearized in order to predict the next increment.This linearization requires to write the macroscopic constitutive law in the non-linear form of Eq. 9. The equilibrium at the macroscopic level in the absence of body forces is given by the Eq. 10.

∆σ x = C t (x) : ∆ε (x) ∀x ∈ V, (9)
div x ∆σ x = 0 ∀x ∈ V, (10)
where ∆σ(x) is the macroscopic stress tensor associated with the point x of the macrostructure at each macroscopic strain increment.

The relationship between the macroscopic stress and strain is given in Voigt notation in Eq. 11.

The macroscopic tangent operator C t is recovered by computing the macroscopic stress resulting from the six elementary strain states written in Eq. 12 (also in Voigt notation) at each macroscopic strain increment:

                                         ∆σ 1 ∆σ 2 ∆σ 3 ∆σ 4 ∆σ 5 ∆σ 6                                          =                                          C t,
                                         ×                                          ∆ε 1 ∆ε 2 ∆ε 3 2 ∆ε 4 2 ∆ε 5 2 ∆ε 6                                          (11)
                                         ∆ε (1) = (K 0 0 0 0 0) T ∆ε (2) = (0 K 0 0 0 0) T ∆ε (3) = (0 0 K 0 0 0) T ∆ε (4) = (0 0 0 K 0 0) T ∆ε (5) = (0 0 0 0 K 0) T ∆ε (6) = (0 0 0 0 0 K) T (12)
Then, the ij component of the tangent operator is given by the i th component of the stress vector calculated with the j th elementary strain state, divided by the j th component of the strain vector of the elementary strain state:

C t, i j = ∆σ ( j) i K , i, j = 1, 2, 3, 4, 5, 6. (13)
Usually, K is chosen to be equal to 1.

Local elasto-viscoplastic behavior with ductile damage for the matrix

The constitutive law of the matrix material is defined through a thermodynamically based phenomenological model for viscoplasticity and ductile damage in semi-crystalline polymers [START_REF] Lemaitre | Mechanics of solid materials[END_REF][START_REF] Krairi | A thermodynamically-based constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage[END_REF][START_REF] Praud | Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF]. These materials exhibit a dissipative behavior that combines solid and fluid properties with some apparent stiffness reduction. The model is described by the rheological scheme given in the Fig. 3. It is composed of: one single linear spring, subjected to an elastic strain ε e , and a viscoplastic branch, subjected to a viscoplastic strain ε p which consists of a frictional element, a non-linear spring and a non-linear dash-pot. The linear spring and the viscoplastic branch are positioned in series. The model is formulated within the thermodynamics framework [START_REF] Lemaitre | Mechanics of solid materials[END_REF][START_REF] Praud | Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF]. The state laws are obtained by differentiation of the Helmholtz potential with respect to the state variables. This potential is formulated as the sum of the stored energies of the spring and the viscoplastic branch.

ρψ ε, r, ε p , D = 1 2 ε -ε p : (1 -D) C e : ε -ε p + r 0 R(ξ) dξ (14)
The internal state variables ε p , r and D represent the viscoplastic strain, effective equivalent viscoplastic strain variable and the damage variable respectively. C e is the initial fourth order stiffness tensor of the single spring, classically defined for bulk isotropic materials. R is the hardening function, chosen under the form of the power law function, that must be increasing, positive and vanishes at r = 0:

R (r) = Kr n , (15)
where K and n are the viscoplastic material parameters. According to the second principle of thermodynamics, dissipation is always positive or null (Clausius Duhem inequality). Assuming that the mechanical and thermal dissipations are uncoupled, the rate of the mechanical dissipated energy Φ is positive or zero and is given by the difference between the strain energy rate Ẇε and the stored energy rate ρ ψ (Eq. 16).

Φ = Ẇε -ρ ψ = σ : ε -ρ ∂ψ ∂ε : ε + ∂ψ ∂ε p : εp + ∂ψ ∂r : ṙ + ∂ψ ∂D : Ḋ = σ : εp -Rṙ + Y Ḋ ≥ 0. ( 16
)
The viscoplasticity and damage are considered to be coupled phenomena [START_REF] Lemaitre | Coupled elasto-plasticity and damage constitutive equations[END_REF][START_REF] Krairi | A thermodynamically-based constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage[END_REF]. Consequently, the evolution of ε p , r and D are described by the normality of a convex indicative function that satisfies the above inequality:

F (σ, R, Y; D) = eq(σ) 1 -D -R -R 0 f (σ, R; D) + S (β + 1) (1 -D) Y S β+1 f D (Y;D) (17)
In the last expression, the term f (σ, R; D) denotes the yield criterion function which activates the mechanism (ṙ > 0 if f > 0, else ṙ = 0). The function f is expressed in the effective stress space. f D is an additive term that takes into account the evolution of the damage at the same time as the viscoplasticity. eq(σ) denotes the von Mises equivalent stress, R 0 denotes the yield threshold, while S and β are damage related material parameters. The viscous effect is introduced by considering a relation between the positive part of f and ṙ through a function Q. This function is chosen under the form of the power law:

f + = Q (ṙ) , Q (ṙ) = Hṙ m (18)
where H and m are the material parameters. The function Q (ṙ) must be increasing, positive and

null at ṙ = 0.
This type of model allows to capture some well known effects of thermoplastic polymers, namely the rate effect through the creep and relaxation phenomena, as well as the stiffness reduction due to the ductile damage. Tab. 2 summarizes the thermodynamic variables, the evolution laws and the von Mises type viscoplastic criterion of the model. In the table, Dev(σ) denotes the deviatoric part of the stress.

ε σ = ρ ∂ψ ∂ε = (1 -D) C e : ε -ε p State variables Associated variables Evolution laws r R = ρ ∂ψ ∂r = R (r) ṙ = - ∂F ∂R λ = λ ε p -σ = ρ ∂ψ ∂ε p εp = ∂F ∂σ λ = 3 2 Dev(σ) eq(σ) ṙ 1 -D D Y = ρ ∂ψ ∂D Ḋ = ∂F ∂Y λ = Y S β ṙ 1 -D Multiplier Criterion Active ( λ > 0 if f > 0) λ = r f (σ, R; D) = eq(σ) 1 -D -R -R 0 f + = Q (ṙ)

Multi scale FE computation and numerical implementation

To predict the macroscopic behavior of a composite structure, taking into account the effect of the microstructure, homogenization scheme within the framework of a FE 2 is an accurate solution. According to Feyel [START_REF] Frédéric | Multiscale fe2 elastoviscoplastic analysis of composite structures[END_REF], this approach considers that the macroscopic problem and the microscopic heterogeneous unit cell are solved simultaneously. On the macroscopic scale, the material is assumed as a homogenized medium with non-linear behavior. The macroscopic response is calculated by solving an appropriate periodic boundary value problem at the microscopic level within a homogenization scheme. The important macroscopic information (strain) passes to the unit cell through the constraint drivers. The concept of constraint drivers is explained in the next subsection.

It is pointed out that the response at the macroscopic scale is obtained by the homogenization process and is frequently called "homogenized". The macroscopic fields and tangent moduli depend on the microscopic response at each unit cell. Since the macroscopic strains are heterogeneous in the structure, the homogenized response varies at every macroscopic point, providing a type of spatial heterogeneity.

Unit cell computations for periodic homogenization using the concept of constraint drivers

The method of constraint drivers is a numerical technique which allows to apply any state of macroscopic stress, strain or even mixed stress/strain on a periodic finite element unit cell. More detailed exposition about this concept is given in [START_REF] Li | On the unit cell for micromechanical analysis of fibre-reinforced composites[END_REF][START_REF] Li | General unit cells for micromechanical analyses of unidirectional composites[END_REF][START_REF] Shuguang | Unit cells for micromechanical analyses of particle-reinforced composites[END_REF].

In the finite element framework a unit cell for periodic media should be associated with periodic mesh. This means that for each border node, there must be another node at the same relative position on the opposite side of the unit cell. The aim of the constraint drivers in a periodic homogenization approach is to apply a macroscopic strain ε on the unit cell, taking into account the periodic boundary conditions. In practice, a displacement gradient is applied between each pair of opposite parallel border nodes (denoted by the indices i and j). This gradient is directly related to the macroscopic strain tensor ε i j by the following general kinematic relationship:

u i = u j ⇐⇒ u i -u j = ε . (x i -x j ) ∀x ∈ V (19)
The proposed method introduces the six components of the macroscopic strain tensor as additional degrees of freedom (constraint drivers) that are linked to the mesh of the unit cell using the kinematic equation 19. The displacements of these additional degrees of freedom, noted as , respectively, and they permit to recover directly the corresponding components of the macroscopic stress tensor (Fig. 4) at the end of the unit cell calculations. Dividing the dual force by the unit cell volume leads to the corresponding macroscopic stresses.

u

Concept and numerical algorithm of FE 2 method

After defining the concept of constraint drivers, the implementation of a two-scale finite element approach is the next step in the computational homogenization framework. The proposed method lies within the general category of multi-scale models. In this method the macroscopic constitutive behavior is calculated directly from the unit cell, providing the geometry and the phenomenological constitutive equations of each constituent. The FE 2 method consists of three main steps according to [START_REF] Frédéric | Multiscale fe2 elastoviscoplastic analysis of composite structures[END_REF]:

(1) A geometrical description and a FE model of the unit cell.

(2) The local constitutive laws expressing the response of each component of the composite within the unit cell.

(3) Scale transition relationships that define the connection between the microscopic and the macroscopic fields (stress and strain).

The scale transition is provided by the concept of homogenization theory, using volume averaging of microscopic quantities of the unit cell, which is solved thanks to periodic boundary conditions. The macroscopic fields (stress and strain) are introduced in a unit cell by using the six additional degrees of freedom (constraint drivers), that are linked with the boundaries through the kinematic equations [Eq. 19]. The macroscopic behavior of a 3D composite structure is computed by considering that the material response of each macroscopic integration point is established from the homogenization of a unit cell that is connected to each macroscopic integration point. Each unit cell contains the local constitutive laws of different phases and the geometrical characteristics of the microstructure.

The FE 2 approach presented here has been developed using an implicit resolution scheme, with the use of a Newton-Raphson algorithm, that solves the non-linear problems at the two scales. At each macroscopic integration point, the macroscopic stress and the macroscopic tangent operator are computed for the calculated macroscopic strain at each time increment, by solving iteratively a FE problem at the microscopic scale.

Concept of transition between scales in FE 2 computations

In the framework of FE 2 modelling the global resolution step is performed at each time increment by solving a local equilibrium problem at each macroscopic integration point. At each step, the microscopic problem is solved by applying the macroscopic strain increment to the unit cell through the periodic boundary conditions. The system of equations in the linearized incremental form is given as follows:

                     ∆σ (x) = C t (x) : ∆ε (x) ∀x ∈ V, div x (∆σ (x)) = 0 ∀x ∈ V, ∆u i -∆u j = ∆ε . x i -x j ∀x ∈ V (20)
By using the developed user subroutine at the microscopic scale which contains the non-linear local behavior of the constituents, the microscopic stress, tangent operator and internal state variables V k are computed at every microscopic point. The macroscopic stress σ is then computed through volume averaging of the microscopic stresses, and the local tangent operators of all microscopic points are utilized to obtain the macroscopic tangent operator C t by solving six elastic-type loading cases with the elementary strain states described in subsection 2.1.2. The internal state variables and the local stress are saved as initial conditions for the next time increment. Once the macroscopic quantities σ and C t are computed, the analysis at the macroscopic level is then performed and the macroscopic strain increment ∆ε is provided by the Finite Element Analyses Package ABAQUS at every macroscopic point through the global equilibrium resolution. This information is passed to the macroscopic scale by using a user defined constitutive model (denoted here as Meta-UMAT) that represent the behavior of a macroscopic material point and contains the unit cell equations and hence the process returns to the local problem. The iterative procedure inside the Meta-UMAT is depicted in Fig. 5. The loop is repeated until numerical convergence is achieved in both micro and macro-scales numerical problems. After the convergence, the analysis proceeds to the next time step. Both the Meta-UMAT and the structural analysis in the macroscopic level define the FE 2 approach.

Algorithm of FE 2 and parallel calculation

The algorithm of the FE 2 computational strategy for the non-linear case in ABAQUS/Standard is presented in Fig. 6.

As shown in Figures 5 and6, the macroscopic problem is solved at each increment in a linearized manner, considering the homogenized tangentmod modulus C t . The elastic prediction -inelastic correction is performed at the scale of the constituents laws (Micro-UMAT) using the well-known "return mapping algorithm -convex cutting plane" scheme [START_REF] Simo | Computational Inelasticity[END_REF].

The aim of the FE 2 approach is to perform structural numerical simulations, thus reduction of the computational time is of outmost importance. Since the FE 2 homogenization requires very costly computations, parallel calculation procedures for running the analysis on multiple CPUs are unavoidable.

Parallel implementation of the FE 2 code in ABAQUS/Standard

It is known that the FE 2 computation is expensive in terms of CPU time, caused by the transition between the two scales and the degree of freedom number of the microscopic and the macroscopic models. To reduce this computational time, a parallel implementation of the FE 2 procedure is set-up in ABAQUS/Standard. All the Finite Element Analyses of the unit cells within a sin-initialisation -Apply the PBCs on the unit cell.

-Compute the initial macroscopic tangent modulus C t .

Macro-level -Solve the macro problem.

-Get the macroscopic strain increment ∆ε n+1 .

Micro-level

-Python script for the micro problem.

-Compute the local fields σ, ε, C t .

-Compute the macroscopic stress σ.

-Compute the macroscopic tangent modulus C t .

Global check convergence

Next increment n=n+1 Update all fields: gle macroscopic element (one per integration point) is sent on a single computation node (a set of processors) and their are solved iteratively. The computations that correspond to each macroscopic element are solved in different computations nodes. Thus, theoretically the parallelization can be performed simultaneously on every element. In practice the parallel computation is limited to the number of available calculation nodes. Note that the computations of every microscopic Finite Element Analysis can also be computed in parallel within the computation node if it possess several processors (which is often the case) and this parallelization process is governed by the Finite Element Analyses Package ABAQUS. In practice, the Meta-UMAT calls an appropriate python script that solves the local problem (including the computation of the macroscopic tangent modulus) at each macroscopic integration point, with the use of the microscopic UMAT, which contains the local non-linear behavior of the constituents (Fig. 7). Afterwards, the global solver of ABAQUS checks that all calculations at different processors are completed before proceeding to the resolution of the macroscopic problem, before passing to the next time increment, or the next macroscopic iteration.

V k , σ, ε, C t , σ, C t . C t ∆ε n+1 C t , σ non yes

Implementation of the microscopic problem in ABAQUS

With regard to the microscopic problem, as mentioned previously, the Meta-UMAT executes a properly designed python script in each macroscopic integration point of the composite structure with the first macroscopic strain increment given by ABAQUS.

The periodic boundary conditions (PBCs) and the macroscopic strain are applied on the unit cell by means of the python script at each time increment, since this last information is given at each integration point from the prediction of the strain increment that should satisfy the global equilibrium. The script also calls the solver ABAQUS to solve the microscopic Finite Element Analyses, which utilize the microscopic user subroutines that contains the local constitutive laws of the constituents, in order to obtain the microscopic response through a return mapping iterative process.

Once the local equilibrium is achieved, the local response (σ and C t ) are computed. Then, the macroscopic stress is recovered as a reaction force divided by the unit cell volume on the constraint drivers (section 3.1). The macroscopic tangent modulus is calculated by mapping the local tangent moduli on the unit cell through the six elementary strain states. Through the python script, the macroscopic quantities (σ and C t ) are calculated and transferred to in the Meta-UMAT. At this point, the global equilibrium is checked, if the convergence is reached, we proceed to the next time increment n+1.

Applications and Capabilities of the FE 2 framework

In order to validate the two-scale computational approach within the framework of 3D nonlinear composite structures, two test cases have been addressed: the first one is a periodic multilayer composite structure with non-linear, elastoplastic phases. It has been demonstrated that the use of an incremental linearized temporal integration approach, there exist a semi-analytical solution for this problem [START_REF] Chatzigeorgiou | Computational micro to macro transitions for shape memory alloy composites using periodic homogenization[END_REF]. This test case is utilized as a validation of the implementation of the FE 2 framework. The second one is the simulation of three-dimensional composite structure, with a two-phase microstructure: A matrix phase that exhibits a coupled elastoviscoplastic with ductile damage response, reinforced by short glass fibres. The results of such multi-scale simulation are compared with a legacy modelling approach, i.e. the use of an incremental Mori-Tanaka scheme [START_REF] Lagoudas | Elastoplastic behavior of metal matrix composites based on incremental plasticity and the mori-tanaka averaging scheme[END_REF].

Comparison with semi-analytical homogenization method for elastoplastic multilayer composites

The multi-scale structure simulated is presented in Fig. 8 and is composed at the microscopic scale of a periodic stack of two different layers, one with an elastic response (superscript e) and the second one with an elastic-plastic response (superscript p). The volume fraction of the two phases is equal, i.e. c e = c p = 0.5. The macroscopic shape of the structure is a cuboid. For the elastic-plastic phase, the plastic yield criterion is given by:

f p (σ, p) = eq(σ) -R(p) -R 0 ≤ 0. ( 21
)
where eq(σ) is the equivalent Von Mises stress and R 0 is the yield threshold. The hardening function R(p) is chosen under the form of a power law [START_REF] Lemaitre | Mechanics of solid materials[END_REF]:

R(p) = K.p n (22)
where K and n are material parameters. p is the accumulated plastic strain. The material parameters of the two phases are given in Tab. 3. As discussed in the Section 3.1, periodic boundary conditions are applied to the unit cell of the multilayer material. The macroscopic boundary conditions imposed correspond to a pure shear loading and are such that the relationship between the displacement at the boundary is u 0 = ε 0 n, n being the outward normal of the surface and all the components of the tensor ε 0 are zero except ε 12 = ε 21 , see Fig. 9-b. Note that under such conditions, the numerical results of the two Finite Element Analyses should be mesh-independent, since homogeneous fields are considered in all the phases. The results of the two approaches, the FE 2 and the semi-analytical are identical (Fig. 10), which demonstrates the capability of the computational method to predict the response of 3D non-linear multi-scale composite structures.

3D structure (Meuwissen) with short fibre reinforced composite

To demonstrate the capabilities of the FE 2 approach to identify the overall behavior of 3D composite structures close to parts that are commonly manufactured, the second test case is performed on structure where heterogeneous strain and stress fields are observed during a tensile load field. the composite material is considered as an elastoviscoplastic polymer matrix with ductile damage, reinforced by aligned glass short fibres arranged in a periodic hexagonal array (Fig. 11). The volume fractions of the matrix and the fibres are V m = 0.925 and V f = 0.075 respectively, while the aspect ratio for the elliptic fibre is (4, 1, 1). The fibres elastic properties are the following: a Young's modulus E f = 72000 MPa and a Poisson's ratio ν f = 0.26. The material properties of the matrix phase are listed in Tab. 4. It should be mentioned that these material parameters are motivated by the work of [START_REF] Praud | Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation[END_REF], but they do not consider the viscoelastic response, which is taken into account in that article. Thus, the material properties are related to viscoplastic behavior coupled to damage in polymeric media. The structure presented in the Fig. 12-a is clamped at the left side and subjected to the loading path of Fig. 12-b at the right side. The displacement controlled path consists in three loading steps with different velocities (u (1) x = 1 mm.s -1 , u(2)

x = 0.2 mm.s -1 , u (3)
x = 0.8 mm.s -1 ) followed by an unloading stage at a displacement rate of (u (4) x = 2 mm.s -1 ). The results of the full-field FE 2 method are compared with those obtained by using the incremental mean-field Mori-Tanaka method. This method has been widely utilized for the simulation of composites [START_REF] Doghri | Homogenization of two-phase elasto-plastic composite materials and structures: Study of tangent operators, cyclic plasticity and numerical algorithms[END_REF] as well as smart structures [START_REF] Piotrowski | Modeling of niobium precipitates effect on the ni47ti44nb9 shape memory alloy behavior[END_REF]. Such homogenization scheme is however considered valid under specific cases [START_REF] Lagoudas | Elastoplastic behavior of metal matrix composites based on incremental plasticity and the mori-tanaka averaging scheme[END_REF], and some specific corrections might be required [START_REF] Chaboche | On the capabilities of mean-field approaches for the description of plasticity in metal matrix composites[END_REF][START_REF] Brassart | Homogenization of elasto-plastic composites coupled with a nonlinear finite element analysis of the equivalent inclusion problem[END_REF]. Since the proposed corrections are not unique and depends on the type of composites, the regular incremental method is employed, where the linearized problem is written in term of the anisotropic algorithmic tangent modulus of the non-linear phases [START_REF] Lagoudas | Elastoplastic behavior of metal matrix composites based on incremental plasticity and the mori-tanaka averaging scheme[END_REF]. The advantage of the Mori-Tanaka scheme relies in its computational efficiency, since it is a semi-analytical method and accounts for the material non-linearities only on an average sense and not at every local microscopic point in the unit cell. The overall load-displacement response computed using FE 2 approach is shown in Fig. 13 and compared to the global response predicted by the mean-field Mori-Tanaka approach. As expected the both approaches predict comparable responses notably in the elastic regime. However, for the viscoplastic regime, the mean-field based simulation does not capture well the strain rate change of the applied loading path and provides stiffer response than the full-field FE 2 . This aspect is known and occurs when one phase exhibits a non-linear behavior. Similar observations have been reported to the literature especially when the matrix phase behaves as viscoelastic-viscoplastic media [START_REF] Miled | Micromechanical modeling of coupled viscoelastic-viscoplastic composites based on an incrementally affine formulation[END_REF][START_REF] Chaboche | On the capabilities of mean-field approaches for the description of plasticity in metal matrix composites[END_REF].

The authors proposed specific numerical formulations to address this limit of mean-field based methods. Fig. 14 demonstrates stress-strain curves at the macroscopic point A (Fig. 15-a). Due to the semi-analytical form of the Mori-Tanaka method, the computations are faster than in the FE 2 method but it requires a smaller time increment. The results indicate that the response of the two approaches describe the changing of the rate loading caused by the viscous behavior of polymer matrix, but it is clearly shown that the Mori-Tanaka response misdescribed this phenomena because it is more rigid with a considerable loss of plasticity as expected, compared to the FE 2 .

The results illustrate that the response of the composite is highly influenced by the presence of the matrix, exhibiting both viscoplastic response through relaxation phenomena, as well as stiffness reduction during unloading due to the ductile damage.

It is worth noticing that the inelastic characteristics of the different phases are mainly taken into account in the microscale and, accordingly, the unit cell is adequately meshed (6857 elements).

The authors have performed several analyses at different meshes of the macroscopic structure and have confirmed that the chosen meshing of 100 elements was sufficient for the purposes of the manuscript.

At a characteristic critical point of the structure (centre of one notch), the deformed macro-scale structure and the microscopic stress response (component 11) of the unit cell that represent a macroscopic integration point A are shown in Fig. 15. It is clear that at such critical material point, the adopted incremental Mori-Tanaka scheme do not predict the local response with a sufficient accuracy to be able to utilize such results for the computation of damage evolution of fatigue life predictions, which are unavoidable in the case of most load-bearing application of composite structures. Even if the mesh convergence is difficult to reach to obtain exact results, the FE 2 framework could provide a standard of predictability that is much higher than the mean-fields methods when a composite is simulated, where the matrix present a strongly non-linear response. It deserves to be mentioned that the periodic homogenization gives excellent results for 3D structures, and the numerical accuracy depends on that of the FE calculations. However, when addressing plate or shell structures, the periodic homogenization requires proper modifications, as described in [START_REF] Kalamkarov | Analysis, design and optimization of composite structures[END_REF], due to the loss of periodicity in the thickness direction (out of plane). This results in less accurate prediction for the out of plane Poisson ratio. Nevertheless, the out of plane periodicity can be reasonably assumed when the microstructure contains high number of fibers or layers in the thickness direction.

Complex 3D structure with corner shape

In this section, a second 3D composite structure is simulated in order to illustrate the capability and the flexibility of the approach, when more complex boundary conditions are applied to the macroscopic structure. The modelled structure consists of a 3D part having a corner shape (Fig. 16-a). It is made of a thermoplastic aligned short fibre reinforced composite in which the matrix and reinforcement phases exhibit the same behavior as in Section 4.2. The structure is clamped at the bottom side and subjected to a normal uniform displacement path at left side (Fig. 16-b).

The displacement controlled path consists in two loading steps with different displacement rates (u (1) x = 2.1875 mm.s -1 , u(2) x = 0.15625 mm.s -1 ) and an unloading step at a displacement rate of (u (3) x = 0.9375 mm.s -1 ). In Figs. 17 and 18, the whole response of the composite in terms of macroscopic stress vs strain are depicted at two distinct points A and B (Fig. 19-c). The approach is able to reproduce the effect of such microstructure on the overall response of the composite, as on the most stressed point shown in Fig. 19. Indeed, on the clamped part (point A), it is clear that the structure is subjected to a tensile load according to the 22 direction, and a shear load in the direction 12.

These results are attributable primarily by the macroscopic boundary conditions. Furthermore, for the point B, a high stress value in the direction of loading was noticed. Fig. 19-c For the parallelization procedure and with the same number of increments (42 increments), the computation becomes 18 times faster than the non parallel solution. The actual computational time of the analysis performed on 18 processors was approximately 72h for a macroscopic structure containing 90 elements of type C3D8 with 6857 microscopic elements of type C3D4.

Conclusions and further work

This work presents a non-linear three-dimensional two-scale finite element (FE 2 ) framework fully integrated in the Finite Element Analysis Package ABAQUS/Standard, using parallel computation. The main advantage of the method is that it does not require an analytical form for the constitutive law at the macro-scale, while accounting for the microstructural effects and the local behaviors. It can integrate any kind of periodic microstucture with any type of non-linear behavior of the reinforcement (fibres and/or particles) and the matrix (plastic, viscoelastic, viscoplastic and damage).

The multi-scale strategy has been tested on three independent numerical examples: In the first example, a 3D multilayer composite structure with elastoplastic phases is simulated and compared with semi analytical solution, to validate the numerical implementation. In the second example, a short glass fibre reinforced composite with elastoviscoplastic-damageable matrix under complex loading is examined through the FE 2 strategy and the results are compared to those obtained by the Mori-Tanaka method. The obtained responses were in agreement with those presented in the literature in similar cases, and highlight the importance of utilizing full-field method for a generic modelling strategy with high predictability capabilities. In the third example, 3D a complex composite structure with corner shape is simulated in which the microstucture is made of an elastoviscoplastic matrix with ductile damage reinforced by short glass fibre. The capabilities of such approach to reproduce the effect of such microstructure on the macrostructure response at each macroscopic integration point has been demonstrated. The response of the structure clearly highlights creep and relaxation phenomena, which are characteristic for rate dependent responses. This viscous behavior and the stiffness reduction observed during unloading have been induced by the viscoplastic nature of the polymer matrix. It worth noticing that for composites where the matrix is viscoplastic material, the Mori-Tanaka method under proper modifications can provide quite accurate results [START_REF] Mercier | Homogenization of elastic-viscoplastic heterogeneous materials: Self-consistent and mori-tanaka schemes[END_REF][START_REF] Miled | Micromechanical modeling of coupled viscoelastic-viscoplastic composites based on an incrementally affine formulation[END_REF] compared to the full-field based approach.

A last advantage of this approach is that it can be extended to predict the overall fully coupled thermomechanical response of 3D composite structures [START_REF] Berthelsen | Computational homogenisation for thermoviscoplasticity: application to thermally sprayed coatings[END_REF][START_REF] Chatzigeorgiou | Computational micro to macro transitions for shape memory alloy composites using periodic homogenization[END_REF] with more complex mechanisms between fibres-matrix as interfacial damage mechanisms. Such fully-coupled analyses on multiscale structures should be of a high interest for industrial applications that are usually computed with commercial finite element analyses packages.

Figure 1 :

 1 Figure 1: Schematic representation of the homogenization computational.

Figure 2 :

 2 Figure 2: Definition of the displacement field as the sum of an affine part and a periodic fluctuation

Figure 3 :

 3 Figure 3: Rheological scheme of the viscoplastic behavior and ductile damage [11].

Figure 4 :

 4 Figure 4: Connection of the constraint drivers with the unit cell.

Figure 5 :

 5 Figure 5: Meta-UMAT for the overall response computation of the composite using FE 2 approach at time increment n+1.

Figure 6 :

 6 Figure 6: The flow chart of the two scales FE 2 algorithm in ABAQUS/Standard for non-linear case.

Figure 7 :

 7 Figure 7: Parallelization steps of the FE 2 code.

Figure 8 :

 8 Figure 8: Multilayer composite structure with their microstructure associated with each macroscopic integration point.

Figure 9 :

 9 Figure 9: Multilayer composite structure under shear loading path.

Figure 10 :

 10 Figure 10: Comparison of the numerical result of FE 2 approach with semi-analytical solution on multilayer witch elastoplastic phases in term of macroscopic stress-strain response.

  (a) Mesh of the entire unit cell.(b) Short fibres reinforcement.

Figure 11 :

 11 Figure 11: Composite microstructure.

  (a) Mesh of the entire 3D composite structure. (b) Applied loading path.

Figure 12 :

 12 Figure 12: Tensile and compression test on the 3D Meuwissen test tube.

Figure 13 :Figure 14 :

 1314 Figure 13: The overall load-displacement response of the structure in the directions 11. Comparison between the FE 2 and the Mori-Tanaka solution.

( a )

 a Macroscopic stress field of the composite structure. (b) Microscopic stress field of the microstructure at point A.

Figure 15 :

 15 Figure 15: FE 2 solution with ABAQUS/Standard (component 11).

  shows the stress response (component 11) of the macroscopic structure and the resulting microscopic stress in the two unit cells situated at two different macroscopic integration points A and B (Figs. 19-a and 19-b respectively). The response of the composite is highly affected by the matrix behavior through the relaxation phenomena caused by the change of the loading rate. The apparent stiffness reduction during the unloading caused by the ductile damage in the matrix is clearly observed.

  (a) Macroscopic mesh of the 3D composite structure. (b) Applied loading path.

Figure 16 :

 16 Figure 16: Tensile and compression test on the 3D composite structure with corner.

- 1 .

 1 75 -1.50 -1.25 -1.00 -0.75 -0.50 -0.25 0

Figure 17 :

 17 Figure 17: Macroscopic response of the composite at point A in term of stress-strain in the directions 11, 22, 33 and shear 12.

Figure 18 :

 18 Figure 18: Macroscopic response of the composite at point B in term of stress-strain in the directions 11, 22, 33 and shear 12.

  (a) Microscopic stress field of the microstructure at point A. (b) Microscopic stress field of the microstructure at point B. (c) Macroscopic stress field of the 3D composite structure (component 11).

Figure 19 :

 19 Figure 19: FE 2 solution with ABAQUS/Standard (component 11).

Table 1 :

 1 Macroscopic and microscopic scale transition[START_REF] Praud | Modélisation multi-échelle des composites tissés à matrice thermoplastique sous chargements cycliques non proportionnels[END_REF]

Table 2 :

 2 State and evolution laws

	Observable state variable	Associated variable

  cd 11 , u cd 22 , u cd 33 , u cd 12 , u cd 13 and u cd 23 , take the values of each component of the macroscopic strain tensor ε 11 , ε 22 , ε 33 , 2ε 12 , 2ε 13 and 2ε 23 , respectively. The dual forces of the constrain drivers are noted as F cd 11 , F cd 22 , F cd 33 , F cd 12 , F cd 13 and F cd 23

Table 3 :

 3 Material parameters for the two phases

		Elastic-plastic phase
	Parameter	value	unit
	E p	2000	MPa
	ν p	0.3	-
	R 0	10	MPa
	K	60.0	MPa
	n	0.15	-
		Elastic phase
	E e	6000	MPa
	ν e	0.2	-

Table 4 :

 4 Material parameters for polymer matrix

	Parameter	value	unit
	E m	1680	MPa
	ν m	0.3	-
	R 0	10	MPa
	K	365.0	MPa
	n	0.39	-
	H	180.0	MPa.s
	m	0.3	-
	S	6.0	MPa
	β	-1.70	-


Paulin Jacquot
email: paulin.jacquot@polytechnique.edu

Cheng Wan
email: cheng.wan.2005@polytechnique.org

† Cheng

Routing Game on Parallel Networks: the Convergence of Atomic to Nonatomic

We consider an instance of a nonatomic routing game. We assume that the network is parallel, that is, constituted of only two nodes, an origin and a destination. We consider infinitesimal players that have a symmetric network cost, but are heterogeneous through their set of feasible strategies and their individual utilities. We show that if an atomic routing game instance is correctly defined to approximate the nonatomic instance, then an atomic Nash Equilibrium will approximate the nonatomic Wardrop Equilibrium. We give explicit bounds on the distance between the equilibria according to the parameters of the atomic instance. This approximation gives a method to compute the Wardrop equilibrium at an arbitrary precision.

Introduction

Motivation. Network routing games were first considered by Rosenthal [START_REF] Rosenthal | The network equilibrium problem in integers[END_REF] in their "atomic unsplittable" version, where a finite set of players share a network subject to congestion. Routing games found later on many practical applications not only in transport [START_REF] Marcotte | Equilibria with infinitely many differentiated classes of customers[END_REF][START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF], but also in communications [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF], distributed computing [START_REF] Altman | Nash equilibria in load balancing in distributed computer systems[END_REF] or energy [START_REF] Atzeni | Demand-side management via distributed energy generation and storage optimization[END_REF]. The different models studied are of three main categories: nonatomic games (where there is a continuum of infinitesimal players), atomic unsplittable games (with a finite number of players, each one choosing a path to her destination), and atomic splittable games (where there is a finite number of players, each one choosing how to split her weight on the set of available paths).

The concept of equilibrium is central in game theory, for it corresponds to a "stable" situation, where no player has interest to deviate. With a finite number of players-an atomic unsplittable game-it is captured by the concept of Nash Equilibrium [START_REF] Nash | Equilibrium points in n-person games[END_REF]. With an infinite number of infinitesimal players-the nonatomic case-the problem is different: deviations from a finite number of players have no impact, which led Wardrop to its definition of equilibria for nonatomic games [START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF]. A typical illustration of the fundamental difference between the nonatomic and atomic splittable routing games is the existence of an exact potential function in the former case, as opposed to the latter [START_REF] Nisan | Algorithmic game theory[END_REF]. However, when one considers the limit game of an atomic splittable game where players become infinitely many, one obtains a nonatomic instance with infinitesimal players, and expects a relationship between the atomic splittable Nash equilibria and the Wardrop equilibrium of the limit nonatomic game. This is the question we address in this paper.

Main results. We propose a quantitative analysis of the link between a nonatomic routing game and a family of related atomic splittable routing games, in which the number of players grows. A novelty from the existing literature is that, for nonatomic instances, we consider a very general setting where players in the continuum [0, 1] have specific convex strategy-sets, the profile of which being given as a mapping from [0, 1] to R T . In addition to the conventional network (congestion) cost, we consider individual utility function which is also heterogeneous among the continuum of players. For a nonatomic game of this form, we formulate the notion of an atomic splittable approximating sequence, composed of instances of atomic splittable games closer and closer to the nonatomic instance. Our main results state the convergence of Nash equilibria (NE) associated to an approximating sequence to the Wardrop equilibrium of the nonatomic instance. In particular, Thm. 11 gives the convergence of aggregate NE flows to the aggregate WE flow in R T in the case of convex and strictly increasing price (or congestion cost) functions without individual utility; Thm. 14 states the convergence of NE to the Wardrop equilibrium in ((R T ) [0,1] , . 2 ) in the case of player-specific strongly concave utility functions. For each result we provide an upper bound on the convergence rate, given from the atomic splittable instances parameters. An implication of these new results concerns the computation of an equilibrium of a nonatomic instance. Although computing an NE is a hard problem in general [START_REF] Koutsoupias | Worst-case equilibria[END_REF], there exists several algorithms to compute an NE through its formulation with finite-dimensional variational inequalities [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF]. For a Wardrop Equilibrium, a similar formulation with infinite-dimensional variational inequalities can be written, but finding a solution is much harder.

Related work. Some results have already been given to quantify the relation between Nash and Wardrop equilibria. Haurie and Marcotte [START_REF] Haurie | On the relationship between nashcournot and wardrop equilibria[END_REF] show that in a sequence of atomic splittable games where atomic splittable players replace themselves smaller and smaller equal-size players with constant total weight, the Nash equilibria converge to the Wardrop equilibrium of a nonatomic game. Their proof is based on the convergence of variational inequalities corresponding to the sequence of Nash equilibria, a technique similar to the one used in this paper. Wan [START_REF] Wan | Coalitions in nonatomic network congestion games[END_REF] generalizes this result to composite games where nonatomic players and atomic splittable players coexist, by allowing the atomic players to replace themselves by players with heterogeneous sizes.

In [START_REF] Gentile | Nash and wardrop equilibria in aggregative games with coupling constraints[END_REF], the authors consider an aggregative game with linear coupling constraints (generalized Nash Equilibria) and show that the Nash Variational equilibrium can be approximated with the Wardrop Variational equilibrium. However, they consider a Wardrop-type equilibrium for a finite number of players: an atomic player considers that her action has no impact on the aggregated profile. They do not study the relation between atomic and nonatomic equilibria, as done in this paper. Finally, Milchtaich [START_REF] Milchtaich | Generic uniqueness of equilibrium in large crowding games[END_REF] studies atomic unsplittable and nonatomic crowding games, where players are of equal weight and each player's payoff depends on her own action and on the number of players choosing the same action. He shows that, if each atomic unsplittable player in an n-person finite game is replaced by m identical replicas with constant total weight, the equilibria generically converge to the unique equilibrium of the corresponding nonatomic game as m goes to infinity. Last, Marcotte and Zhu [START_REF] Marcotte | Equilibria with infinitely many differentiated classes of customers[END_REF] consider nonatomic players with continuous types (leading to a characterization of the Wardrop equilibrium as a infinite-dimensional variational inequality) and studied the equilibrium in an aggregative game with an infinity of nonatomic players, differentiated through a linear parameter in their cost function and their feasibility sets assumed to be convex polyhedra.

Structure. The remaining of the paper is organized as follows: in Sec. 2, we give the definitions of atomic splittable and nonatomic routing games. We recall the associated concepts of Nash and Wardrop equilibria, their characterization via variational inequalities, and sufficient conditions of existence. Then, in Sec. 3, we give the definition of an approximating sequence of a nonatomic game, and we give our two main theorems on the convergence of the sequence of Nash equilibria to a Wardrop equilibrium of the nonatomic game. Last, in Sec. 4 we provide a numerical example of an approximation of a particular nonatomic routing game.

Notation. We use a bold font to denote vectors (e.g. x) as opposed to scalars (e.g. x).

2 Splittable Routing: Atomic and Nonatomic

Atomic Splittable Routing Game

An atomic splittable routing game on parallel arcs is defined with a network constituted of a finite number of parallel links (cf Fig. 1) on which players can load some weight. Each "link" can be thought as a road, a communication channel or a time slot on which each user can put a load or a task. Associated to each link is a cost or "latency" function that depends only of the total load put on this link.

O D t = 1, c 1 t = 2, c 2 • • • t = T, c T

Definition 1. Atomic Splittable Routing Game

An instance G of an atomic splittable routing game is defined by:

• a finite set of players I = {1, . . . , I},

• a finite set of arcs T = {1, . . . , T },

• for each i ∈ I, a feasibility set X i ⊂ R T + , • for each i ∈ I, a utility function u i : X i → R,

• for each t ∈ T , a cost or latency function c t (.) : R → R .

Each atomic player i ∈ I chooses a profile (x i,t ) t∈T in her feasible set X i and minimizes her cost function:

f i (x i , x -i ) := t∈T x i,t c t j∈I x j,t -u i (x i ).
(

) 1
composed of the network cost and her utility, where x -i := (x j ) j =i . The instance G can be written as the tuple:

G = (I, T , X , c, (u i ) i ) , (2)
where

X := X 1 × • • • × X I and c = (c t ) t∈T .
In the remaining of this paper, the notation G will be used for an instance of an atomic game (Def. 1).

Owing to the network cost structure (1), the aggregated load plays a central role. We denote it by X t := i∈I x i,t on each arc t, and denote the associated feasibility set by:

X := X ∈ R T : ∃x ∈ X s.t. i∈I x i = X . (3)
As seen in ( 1), atomic splittable routing games are particular cases of aggregative games: each player's cost function depends on the actions of the others only through the aggregated profile X.

For technical simplification, we make the following assumptions:

Assumption 1. Convex costs Each cost function (c t ) is differentiable, convex and increasing.

Assumption 2. Compact strategy sets For each i ∈ I, the set X i is assumed to be nonempty, convex and compact.

Assumption 3. Concave utilities Each utility function u i is differentiable and concave.

Note that under Asms. 1 and 3, each function f i is convex in x i .

An example that has drawn a particular attention is the class of atomic splittable routing games considered in [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF]. We add player-specific constraints on individual loads on each link, so that the model becomes the following. Example 1. Each player i has a weight E i to split over T . In this case, X i is given as the simplex:

X i = { x i ∈ R T + : t x i,t = E i and x i,t ≤ x i,t ≤ x i,t
} . E i can be the mass of data to be sent over different canals, or an energy to be consumed over a set of time periods [START_REF] Jacquot | Analysis and implementation of an hourly billing mechanism for demand response management[END_REF]. In the energy applications, more complex models include for instance "ramping" constraints r i,t ≤ x i,t+1 -x i,t ≤ r i,t .

Example 2. An important example of utility function is the distance to a preferred profile y i = (y i,t ) t∈T , that is:

u i (x i ) = -ω i x i -y i 2 2 = -ω i t (x i,t -y i,t ) 2 , (4)
where ω i > 0 is the value of player i's preference. Another type of utility function which has found many applications is :

u i (x i ) = -ω i log (1 + t x i,t ) , (5)
which increases with the weight player i can load on T .

Below we recall the central notion of Nash Equilibrium in atomic non-cooperative games.

Definition 2. Nash Equilibrium (NE) An NE of the atomic game G = (I, X , (f i ) i ) is a profile x ∈ X such that for each player i ∈ I: f i ( xi , x-i ) ≤ f i (x i , x-i ), ∀x i ∈ X i . Proposition 1. Variational Formulation of an NE Under Asms. 1 to 3, x ∈ X is an NE of G if and only if: ∀x ∈ X , ∀i ∈ I, ∇ i f i ( xi , x-i ), x i -xi ≥ 0 , (6)
where 6) is the necessary and sufficient first order condition for xi to be a minimum of f i (., x-i ).

∇ i f i ( xi , x-i ) = ∇f i (•, x-i )| •= xi = c t ( Xt ) + xi,t c t ( Xt ) t∈T -∇u i ( xi ). An equivalent condition is: ∀x ∈ X , i∈I ∇ i f i ( xi , x-i ), x i -xi ≥ 0 . Proof. Since x i → f i (x i , x -i ) is convex, (
Def. 1 defines a convex minimization game so that the existence of an NE is a corollary of Rosen's results [START_REF] Rosen | Existence and uniqueness of equilibrium points for concave n-person games[END_REF]:

Theorem 2 (Cor. of [START_REF] Rosen | Existence and uniqueness of equilibrium points for concave n-person games[END_REF]. Existence of an NE If G is an atomic routing congestion game (Def. 1) satisfying Asms. 1 to 3, then there exists an NE of G.

Rosen [START_REF] Rosen | Existence and uniqueness of equilibrium points for concave n-person games[END_REF] gave a uniqueness theorem applying to any convex compact strategy sets, relying on a strong monotonicity condition of the operator (∇ xi f i ) i . For atomic splittable routing games [START_REF] Orda | Competitive routing in multiuser communication networks[END_REF], an NE is not unique in general [START_REF] Bhaskar | Equilibria of atomic flow games are not unique[END_REF]. To our knowledge, for atomic parallel routing games (Def. 1) under Asms. 1 to 3, neither the uniqueness of NE nor a counter example of its uniqueness has been found. However, there are some particular cases where uniqueness has been shown, e.g. [START_REF] Jacquot | Analysis and implementation of an hourly billing mechanism for demand response management[END_REF] for the case of Ex. 1.

However, as we will see in the convergence theorems of Sec. 3, uniqueness of NE is not necessary to ensure the convergence of NE of a sequence of atomic unsplittable games, as any sequence of NE will converge to the unique Wardrop Equilibrium of the nonatomic game considered.

Infinity of Players: the Nonatomic Framework

If there is an infinity of players, the structure of the game changes: the action of a single player has a negligible impact on the aggregated load on each link. To measure the impact of infinitesimal players, we equip real coordinate spaces R k with the usual Lebesgue measure µ.

The set of players is now represented by a continuum Θ = [0, 1]. Each player is of Lebesgue measure 0.

Definition 3. Nonatomic Routing Game

An instance G of a nonatomic routing game is defined by:

• a continuum of players Θ = [0, 1],
• a finite set of arcs T = {1, . . . , T },

• a point-to-set mapping of feasibility sets X . : Θ ⇒ R T + , • for each θ ∈ Θ, a utility function u θ (.) : X θ → R,

• for each t ∈ T , a cost or latency function c t (.) : R → R.

Each nonatomic player θ chooses a profile x θ = (x θ,t ) t∈T in her feasible set X θ and minimizes her cost function:

F θ (x θ , X) := t∈T x θ,t c t X t -u θ (x θ ), (7)
where X t := Θ x θ,t dθ denotes the aggregated load. The nonatomic instance G can be written as the tuple:

G = (Θ, T , (X θ ) θ∈Θ , c, (u θ ) θ∈Θ ) . (8)
For the nonatomic case, we need assumptions stronger than Asms. 2 and 3 for the mappings X . and u., given below: Assumption 4. Nonatomic strategy sets There exists M > 0 such that, for any θ ∈ Θ, X θ is convex, compact and X θ ⊂ B 0 (M ), where B 0 (M ) is the ball of radius M centered at the origin. Moreover, the mapping θ → X θ has a measurable graph Γ

X := {(θ, x) : θ ∈ Θ, x ∈ X θ } ⊂ R T +1 .
Assumption 5. Nonatomic utilities There exists Γ > 0 s.t. for each θ, u θ is differentiable, concave and ∇u θ ∞ < Γ. The function Γ X (θ, x θ ) → u θ (x θ ) is measurable.

Def. 3 and Asms. 4 and 5 give a very general framework. In many models of nonatomic games that have been considered, players are considered homogeneous or with a finite number of classes [START_REF] Nisan | Algorithmic game theory[END_REF]Chapter 18]. Here, players can be heterogeneous through X θ and u θ . Games with heterogeneous players can find many applications, an example being the nonatomic equivalent of Ex. 1: Example 3. Let θ → E θ be a density function which designates the total demand E θ for each player θ ∈ Θ. Consider the nonatomic splittable routing game with feasibility sets

X θ := {x θ ∈ R T + : t x θ,t = E θ }.
As in Ex. 1, one can consider some upper bound x θ,t and lower bound x θ,t for each θ ∈ Θ and each t ∈ T , and add the bounding constraints ∀t ∈ T , x θ,t ≤ x θ,t ≤ x θ,t in the definition of X θ .

Heterogeneity of utility functions can also appear in many practical cases: if we consider the case of preferred profiles given in Ex. 2, members of a population can attribute different values to their cost and their preferences.

Since each player is infinitesimal, her action has a negligible impact on the other players' costs. Wardrop [START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF] extended the notion of equilibrium to the nonatomic case.

Definition 4. Wardrop Equilibrium (WE)

x * ∈ (X θ ) θ is a Wardrop equilibrium of the game G if it is a measurable function from θ to X and for almost all θ ∈ Θ,

F θ (x * θ , X * ) ≤ F θ (x θ , X * ), ∀x θ ∈ X θ ,
where

X * = θ∈Θ x * θ dθ ∈ R T . Proposition 3. Variational formulation of a WE Under Asms. 1, 4 and 5, x * ∈ X is a WE of G iff for almost all θ ∈ Θ: c(X * ) -∇u θ (x * θ ), x θ -x * θ ≥ 0, ∀x θ ∈ X θ . (9)
Proof. Given X * , ( 9) is the necessary and sufficient first order condition for x * θ to be a minimum point of the convex function F θ (., X * ).

According to [START_REF] Jacquot | Analysis and implementation of an hourly billing mechanism for demand response management[END_REF], the monotonicity of c is sufficient to have the VI characterization of the equilibrium in the nonatomic case, as opposed to the atomic case in [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF] where monotonicity and convexity of c are needed.

Theorem 4 (Cor. of Rath, 1992 [16]). Existence of a WE If G is a nonatomic routing congestion game (Def. 3) satisfying Asms. 1, 4 and 5, then G admits a WE.

Proof. The conditions required in [START_REF] Rath | A direct proof of the existence of pure strategy equilibria in games with a continuum of players[END_REF] are satisfied. Note that we only need (c t ) t and (u θ ) θ∈Θ to be continuous functions.

The variational formulation of a WE given in Thm. 3 can be written in the closed form: Theorem 5. Under Asms. 1, 4 and 5, x * ∈ X is a WE of G iff:

θ∈Θ c(X * ) -∇u θ (x * θ ), x θ -x * θ dθ ≥ 0, ∀x ∈ X . (10)
Proof. This follows from Thm. 3. If x * ∈ X is a Wardrop equilibrium so that (9) holds for almost all θ ∈ Θ, then [START_REF] Koutsoupias | Worst-case equilibria[END_REF] follows straightforwardly. Conversely, suppose that x * ∈ X satisfies condition [START_REF] Koutsoupias | Worst-case equilibria[END_REF] but is not a WE of G. Then there must be a subset S of Θ with strictly positive measure such that for each θ ∈ S, (9) does not hold: for each θ ∈ S, there exists

y θ ∈ X θ such that c(X * ) -∇u θ (x * θ ), y θ -x * θ < 0 For each θ ∈ Θ \ S, let y θ := x * θ . Then y = (y θ ) θ∈Θ ∈ X , and θ∈Θ c(X * ) -∇u θ (x * θ ), y θ -x * θ dθ = θ∈S c(X * ) -∇u θ (x * θ ), y θ -x * θ dθ < 0 contradicting (10).
Corrolary 6. In the case where u θ ≡ 0 for all θ ∈ Θ, under Asms. 1 and 4,

x * ∈ X is a WE of G iff: c(X * ), X -X * ≥ 0, ∀X ∈ X . (11)
From the characterization of the WE in Thm. 5 and Thm. 6, we derive Thms. 7 and 8 that state simple conditions ensuring the uniqueness of WE in G.

Theorem 7. Under Asms. 1, 4 and 5, if u θ is strictly concave for each θ ∈ Θ, then G admits a unique WE.

Proof. Suppose that x ∈ X and y ∈ X are both WE of the game. Let X = θ∈Θ x θ dθ and Y = θ∈Θ y θ dθ. Then, according to Theorem 5,

θ∈Θ c(X) -∇u θ (x θ ), y θ -x θ dθ ≥ 0 ( 12
) θ∈Θ c(Y ) -∇u θ (y θ ), x θ -y θ dθ ≥ 0 (13)
By adding ( 12) and ( 13), one has

θ∈Θ c(X) -c(Y ) -∇u θ (x θ ) + ∇u θ (y θ ), y θ -x θ dθ ≥ 0 ⇒ c(X) -c(Y ), θ∈Θ (y θ -x θ )dθ + θ∈Θ -∇u θ (x θ ) + ∇u θ (y θ ), y θ -x θ dθ ≥ 0 ⇒ c(X) -c(Y ), X -Y + θ∈Θ -∇u θ (x θ ) + ∇u θ (y θ ), x θ -y θ dθ ≤ 0
Since for each θ, u θ is strictly concave, ∇u θ is thus strictly monotone. Therefore, for each θ ∈ Θ, -∇u θ (x θ ) + ∇u θ (y θ ), x θ -y θ ≥ 0 and equality holds if and only if

x θ = y θ . Besides, c is monotone, hence c(X) -c(Y ), X -Y ≥ 0. Consequently, c(X) -c(Y ), X -Y + θ∈Θ -∇u θ (x θ ) + ∇u θ (y θ ),
x θ -y θ dθ ≥ 0, and equality holds if and only if for almost all θ ∈ Θ, x θ = y θ . (In this case, X = Y .)

Theorem 8. In the case where u θ ≡ 0 for all θ ∈ Θ, under Asms. 1 and 4, if c = (c t ) T t=1 : [0, M ] T → R T is a strictly monotone operator, then all the WE of G have the same aggregate profile X * ∈ X .

Proof. Suppose that x ∈ X and y ∈ X are both WE of the game. Let X = θ∈Θ x θ dθ and Y = θ∈Θ y θ dθ. Then, according to Corollary 6,

c(X), Y -X ≥ 0 (14) c(Y ), X -Y ≥ 0 ( 15
)
By adding ( 14) and ( 15), one has

c(X) -c(Y ), Y -X ≥ 0
Since c is strictly monotone, c(X) -c(Y ), X -Y ≥ 0 and equality holds if and only X = Y .

Consequently, X = Y .

Remark 1. If for each t ∈ T , c t (.) is (strictly) increasing, then c is a (strictly) monotone operator from [0, M ] T → R T .

One expects that, when the number of players grows very large in an atomic splittable game, the game gets close to a nonatomic game in some sense. We confirm this intuition by showing that, considering a sequence of equilibria of approximating atomic games of a nonatomic instance, the sequence will converge to an equilibrium of the nonatomic instance.

Approximating Nonatomic Games

To approximate the nonatomic game G, the idea consists in finding a sequence of atomic games (G (ν) ) with an increasing number of players, each player representing a "class" of nonatomic players, similar in their parameters.

As the players θ ∈ Θ are differentiated through X θ and u θ , we need to formulate the convergence of feasibility sets and utilities of atomic instances to the nonatomic parameters.

Approximating the nonatomic instance

Definition 5. Atomic Approximating Sequence (AAS) A sequence of atomic games G (ν) = I (ν) , T , X (ν) , c, (u (ν) i ) i is an approximating sequence (AAS) for the nonatomic instance G = Θ, T , (X θ ) θ , c, (u θ ) θ if for each ν ∈ N, there exists a partition of cardinal I (ν) of set Θ, denoted by (Θ

(ν) i ) i∈I (ν) , such that: • I (ν) -→ +∞, • µ (ν) := max i∈I (ν) µ (ν) i -→ 0 where µ (ν) i := µ(Θ (ν) i ) is the Lebesgue measure of subset Θ (ν) i , • δ (ν) := max i∈I (ν) δ (ν) i
-→ 0 where δ i is the Hausdorff distance (denoted by d H ) between nonatomic feasibility sets and the scaled atomic feasibility set:

δ (ν) i := max θ∈Θi d H X θ , 1 µ (ν) i X (ν) i , (16)
• d (ν) := max i∈I (ν) d (ν) i -→ 0 where d i is the L ∞ -distance (in B 0 (M ) → R)
between the gradient of nonatomic utility functions and the scaled atomic utility functions:

d (ν) i = max θ∈Θi max x∈B0(M ) ∇u (ν) i µ (ν) i x -∇u θ (x) 2 . ( 17
)
From Def. 5 it is not trivial to build an AAS of a given nonatomic game G, one can even be unsure that such a sequence exists. However, we will give practical examples in Secs. 3.4.1 and 3.4.2.

A direct result from the assumptions in Def. 5 is that the players become infinitesimal, as stated in Thm. 9.

Lemma 9. If (G (ν) ) ν is an AAS of a nonatomic instance G, then considering the maximal diameter M of X θ , we have:

∀i ∈ I (ν) , ∀x i ∈ X (ν) i , x i 2 ≤ µ (ν) i (M + δ (ν) i ) . ( 18
)
Proof. Let x i ∈ X (ν) i . Let θ ∈ Θ (ν) i
and denote by P X θ the projection on X θ . By definition of δ (ν) i , we get:

x i µ (ν) i -P X θ x i µ (ν) i 2 ≤ δ (ν) i (19) ⇐⇒ x i 2 ≤ µ (ν) i δ (ν) i + P X θ x i µ (ν) i 2 ≤ µ (ν) i (δ (ν) i + M ) . (20)
Lemma 10. If (G (ν) ) ν is an AAS of a nonatomic instance G, then the Hausdorff distance between the aggregated sets X = Θ X . and X (ν) = i∈I (ν) X (ν) i is bounded by:

d H X (ν) , X ≤ δ (ν) . (21)
Proof. Let (x θ ) θ ∈ X be a nonatomic profile. Let P i denote the Euclidean projection on X (ν) i for i ∈ I (ν) and consider

y i := P i Θ (ν) i x θ dθ ∈ X (ν)
i . From ( 16) we have:

Θ x θ dθ - i∈I (ν) y i 2 = i∈I (ν) Θ (ν) i x θ dθ -y i 2 (22) = i∈I (ν) Θ (ν) i x θ - 1 µ (ν) i y i dθ 2 (23) ≤ i∈I (ν) Θ (ν) i x θ - 1 µ (ν) i y i 2 dθ (24) ≤ i∈I (ν) Θ (ν) i δ (ν) i dθ = i∈I (ν) µ (ν) i δ (ν) i ≤ δ (ν) , (25)
which shows that d X, X (ν) ≤ δ (ν) for all X ∈ X . On the other hand, if i∈I (ν) x i ∈ X (ν) , then let us denote by Π θ the Euclidean projection on X θ for θ ∈ Θ, and

y θ = Π θ 1 µ (ν) i x i ∈ X θ for θ ∈ Θ (ν) i . Then we have for all θ ∈ Θ (ν) i , 1 µ (ν) i x i -y θ 2 ≤ δ (ν) i
and we get:

i∈I (ν) x i - Θ y θ dθ 2 ≤ i∈I (ν) Θ (ν) i 1 µ (ν) i x i -y θ dθ 2 (26) ≤ i∈I (ν) Θ (ν) i 1 µ (ν) i x i -y θ 2 dθ (27) ≤ i∈I (ν) µ (ν) i δ (ν) i ≤ δ (ν) , (28)
which shows that d X, X ≤ δ (ν) for all X ∈ X (ν) and concludes the proof.

To ensure the convergence of an AAS, we make the following additional assumptions on costs functions (c t ) t : Assumption 6. Lipschitz continuous costs For each t ∈ T , c t is a Lipschitz continuous function on [0, M ]. There exists C > 0 such that for each t ∈ T , |c t (•)| ≤ C.

Assumption 7.

Strong monotonicity There exists c 0 > 0 such that, for each t ∈ {1, . . . , T },

c t (•) ≥ c 0 on [0, M ]
In the following sections, we differentiate the cases with and without utilities, because we found different convergence results in the two cases.

Players without Utility Functions: Convergence of the Aggregated Equilibrium Profiles

In this section, we assume that u θ ≡ 0 for each θ ∈ Θ.

We give a first result on the approximation of WE by a sequence of NE in Thm. 11.

Theorem 11. Let (G (ν) ) ν be an AAS of a nonatomic instance G, satisfying Asms. 1, 2, 4, 6 and 7. Let ( x(ν) ) a sequence of NE associated to (G (ν) ), and (x * θ ) θ a WE of G. Then:

X(ν) -X * 2 2 ≤ 2 c 0 × B c × δ (ν) + C(M + 1) 2 × µ (ν) ,
where

B c := max x∈B0(M ) c(x) 2 .
Proof. Let P i denote the Euclidean projection onto X (ν) i and Π the projection onto X . We omit the index ν for simplicity. From [START_REF] Marcotte | Equilibria with infinitely many differentiated classes of customers[END_REF], we get:

c(X * ), Π( X) -X * ≥ 0 . ( 29
)
On the other hand, with x * i := Θi x θ dθ, we get from (1):

0 ≤ i∈I c t ( Xt ) + xi,t c t ( Xt ) t∈T , P i (x * i ) -xi (30) = c( X), i P i (x * i ) -X + R( x, x * ) (31) with R( x, x * ) = i xi,t c t ( Xt ) t , P i (x * i ) -xi .
From the Cauchy-Schwartz inequality and Thm. 9, we get:

|R( x, x * )| ≤ i∈I (ν) xi,t c t ( Xt ) t 2 × P i (x * i ) -xi 2 (32) ≤ i∈I (ν) (µ (ν) i (M + δ (ν) i )C × 2(µ (ν) i (M + δ (ν) i )) (33) ≤ 2C(M + 1) 2 max i µ (ν) i . ( 34
)
Besides, with the strong monotonicity of c and from ( 29) and (30):

c 0 X -X * 2 ≤ c( X) -c(X * ), X -X * = c( X), X -X * + c(X * ), X * -X ≤ c( X), X - i P i (x * i ) + c(X * ), X * -Π( X) + c( X), i P i (x * i ) -X * + c(X * ), Π( X) -X ≤ |R( x, x * )| + 0 + 2B c × max i δ i ,
which concludes the proof.

Players with Utility Functions: Convergence of the Individual Equilibrium Profiles

In order to establish a convergence theorem in the presence of utility functions, we make an additional assumption of strong monotonicity on the utility functions stated in Asm. 8. Note that this assumption holds for the utility functions given in Ex. 2.

Assumption 8. Strongly concave utilities For all θ ∈ Θ, u θ is strongly concave on B 0 (M ), uniformly in θ: there exists α > 0 such that for all x, y ∈ B 0 (M ) 2 and any τ ∈]0, 1[ :

u θ ((1 -τ )x + τ y) ≥ (1 -τ )u θ (x) + τ u(y) + α 2 τ (1 -τ ) x -y 2 .
Remark 2. If u θ (x θ ) is α θ -strongly concave, then the negative of its gradient is a strongly monotone operator:

-∇u θ (x θ ) -∇u θ (y θ ), x θ -y θ ≥ α θ x θ -y θ 2 . (35)
We start by showing that, under the additional Asm. 8 on the utility functions, the WE profiles of two nonatomic users within the same subset Θ (ν) i are roughly the same. Proposition 12. Let (G (ν) ) ν be an AAS of a nonatomic instance G and (x * θ ) θ the WE of G satisfying Asms. 1, 4, 5 and 8. Then, if θ, ξ ∈ Θ (ν) i , we have:

x * θ -x * ξ 2 2 ≤ 2 α M d (ν) i + (B c + Γ)δ (ν) i v ξ (x * ξ ) -c(X) ,
which gives the desired result when combined with (40).

This result reveals the role of the strong concavity of utility functions: when α goes to 0, the right hand side of the inequality diverges. This is coherent with the fact that, without utilities, only the aggregated profile matters, so that we cannot have a result such as Thm. 12.

According to Thm. 12, we can obtain a continuity property of the Wardrop equilibrium if we introduce the notion of continuity for the nonatomic game G, relatively to its parameters: Definition 6. Continuity of a nonatomic game The nonatomic instance G = Θ, T , (X θ ) θ , c, (u θ ) θ is said to be continuous at θ ∈ Θ if, for all ε > 0, there exists η > 0 such that:

∀θ ∈ Θ, θ -θ ≤ η ⇒ d H (X θ , X θ ) ≤ ε max x∈X θ ∪X θ ∇u θ (x) -∇u θ (x) 2 ≤ ε . (41)
Then the proof of Thm. 12 shows the following intuitive property:

Proposition 13. Let G = Θ, T , (X θ ) θ , c, (u θ ) θ be a nonatomic instance. If G is continuous at θ 0 ∈ Θ and (x * θ ) θ is a WE of G, then θ → x * θ is continuous at θ 0 .
The next theorem is one of the main results of this paper. It shows that a WE can be approximated by the NE of an atomic approximating sequence. and(x * θ ) θ the WE of G. Under Asms. 1 to 6 and 8, the approximating solution defined by

Theorem 14. Let (G (ν) ) ν be an AAS of a nonatomic instance G. Let ( x(ν) ) a sequence of NE associated to (G (ν) ),
x(ν)

θ := 1 µ (ν) i x(ν) i for θ ∈ Θ (ν) i satisfies: θ∈Θ x(ν) θ -x * θ 2 2 dθ ≤ 2 α B c + Γ)δ (ν) + C(M + 1) 2 µ (ν) + M d (ν) .
Proof. Let ( xi ) i be an NE of G (ν) , and x * ∈ X the WE of G. For the remaining of the proof we ommit the index (ν) for simplicity.

Let us consider the nonatomic profile defined by xθ := 1 µi xi for θ ∈ Θ i , and its projection on the feasibility set ŷθ := P X θ ( xθ ). Similarly, let us consider the atomic profile given by x * i := Θi x * θ dθ for i ∈ I (ν) , and its projection y * i := P Xi (x * i ). For notation simplicity, we denote ∇u θ by v θ . From the strong concavity of v θ and the strong monotonicity of c, we have:

α θ∈Θ x(ν) θ -x * θ 2 2 + c 0 X(ν) -X * 2 2 (42) ≤ Θ c( X) -v θ ( xθ ) -(c(X * ) -v θ (x * θ )) , xθ -x * θ dθ (43) = Θ c( X) -v θ ( xθ ), xθ -x * θ dθ + Θ c(X * ) -v θ (x * θ ), x * θ -xθ dθ . ( 44
)
To bound the second term, we use the characterization of a WE given in Thm. 3, with ŷθ ∈ X θ :

Θ c(X * ) -v θ (x * θ ), x * θ -xθ dθ (45) = Θ c(X * ) -v θ (x * θ ), x * θ -ŷθ dθ + Θ c(X * ) -v θ (x * θ ), ŷθ -xθ dθ (46) ≤ 0 + i∈I (ν) Θi c(X * ) -v θ (x * θ ) 2 × ŷθ -xθ 2 dθ (47) ≤ i∈I (ν) Θi (B c + Γ) × δ i ≤ (B c + Γ) × δ . (48)
To bound the first term of (44), we divide it into two integral terms:

Θ c( X) -v θ ( xθ ), xθ -x * θ dθ (49) = i∈I (ν) Θi c( X) -v i ( xi ), xθ -x * θ dθ + Θi v i ( xi ) -v θ ( xθ ), xθ -x * θ dθ . ( 50
)
The first integral term is bounded using the characterization of a NE given in Thm. 1:

i∈I (ν) Θi c( X) -v i ( xi ), xθ -x * θ dθ (51) = i∈I (ν) c( X) -v i ( xi ), xi -x * i (52) ≤ i∈I (ν) c( X) -v i ( xi ), xi -y * i + i∈I (ν) c( X) -v i ( xi ), y * i -x * i (53) ≤ -R( x, x * ) + i∈I (ν) c( X) -v i ( xi ) 2 × y * i -x * i 2 (54) ≤ 2C(M + 1) 2 µ + (B c + Γ) × 2M δ i∈I (ν) µ i (55) = 2C(M + 1) 2 µ + (B c + Γ) × 2M δ . (56)
For the second integral term, we use the distance between utilities (17):

i∈I (ν) Θi v i ( xi ) -v θ ( xθ ), xθ -x * θ dθ (57) ≤ i∈I (ν) µ i v i ( xi ) -v θ ( xθ ) 2 × xθ -x * θ 2 (58) ≤ i∈I (ν) µ i d i × 2M ≤ d2M . (59)
We conclude the proof by combining (48),( 56) and (59).

As in Thm. 12, the uniform strong concavity of the utility functions plays a key role in the convergence of disaggregated profiles ( x(ν) θ ) ν to the nonatomic WE profile x * .

Construction of an Approximating Sequence

In this section, we give examples of the construction of an AAS for a nonatomic game G, under two particular cases: the case of piecewise continuous functions and, next, the case of finitedimensional parameters.

Piecewise continuous parameters, uniform splitting

In this case, we assume that the parameters of the nonatomic game are piecewise continuous functions of θ ∈ Θ: there exists a finite set of K discontinuity points 0

≤ σ 1 < σ 2 < • • • < σ K ≤ 1,
and the game is uniformly continuous (Def. 6) on (σ k , σ k+1 ), for each k ∈ {0, . . . , K + 1} with the convention σ 0 = 0 and σ K = 1.

For ν ∈ N * , consider the ordered set of I ν cutting points (υ

(ν) i ) Iν i=0 := k ν 0≤k≤ν ∪ {σ k } 1≤k≤K and define the partition (Θ (ν) i ) i∈I (ν) of Θ by: ∀i ∈ {1, . . . , I ν } , Θ (ν) i = [υ (ν) i-1 , υ (ν) i ). ( 60
)
Proposition 15. For ν ∈ N * , consider the atomic game G (ν) defined with I (ν) := {1 . . . I ν }, and for each i ∈ I (ν) :

X (ν) i := µ (ν) i X ῡ(ν) i and u (ν) i := x → µ (ν) i u ῡ(ν) i 1 µ (ν) i x , with ῡ(ν) i = υ (ν) i-1 +υ (ν) i 2 . Then G (ν) ν = I (ν) , T , X (ν) , c, u (ν)
ν is an AAS of the nonatomic game G = (Θ, T , X . , c, (u θ ) θ ).

Proof. We have I (ν) > ν -→ ∞ and for each i ∈ I (ν) , µ(Θ

(ν) i ) ≤ 1
ν -→ 0. The conditions on the feasibility sets and the utility functions are obtained with the piecewise uniform continuity conditions. If we consider a common modulus of uniform continuity η associated to an arbitrary ε > 0, then, for ν large enough, we have, for each i ∈< I (ν) , µ

(ν) i < η. Thus, for all θ ∈ Θ (ν) i , |ῡ (ν) i
-θ| < η, so that from the continuity conditions, we have:

d H X θ , 1 µ (ν) i X (ν) i = d H (X θ , X ῡ(ν) i ) < ε (61)
and max

x∈B0(M ) ∇u (ν) i µ (ν) i x -∇u θ (x) 2 = µ (ν) i µ (ν) i ∇u ῡ(ν) i 1 µ (ν) i µ (ν) i x -∇u θ (x) 2 < ε , (62)
which concludes the proof.

Proof. The proof follows [START_REF] Batson | Combinatorial behavior of extreme points of perturbed polyhedra[END_REF] in several parts, but we extend the result on the compact set B, and drop the irredundancy assumption made in [START_REF] Batson | Combinatorial behavior of extreme points of perturbed polyhedra[END_REF].

For each b, we denote by V (b) the set of vertex of the polyhedron Λ b . Under Assumption 4, V (b) is nonempty for any b ∈ B.

First, as Λ b is a polyhedra, we have Λ b = conv(V (b)) where conv(X) is the convex hull of a set X. As the function x → d(x, Λ b ) defined over Λ b is continuous and convex, by the maximum principle, its maximum over the polyhedron Λ b is achieved on V (b). Thus, we have:

d H (Λ b , Λ b ) = max[ max x∈Λ b d(x, Λ b ) , max x∈Λ b d(Λ b , x)] (64) = max[ max x∈V (b) d(x, Λ b ) , max x∈V (b ) d(Λ b , x)] (65)
≤ max[ max

x∈V (b) d(x, V (b )) , max x∈V (b ) d(V (b), x)] (66) =d H (V (b), V (b )) . (67)
Let's denote by H i (b) the hyperplane {x :

A i x = b i } and by H - i (b) = {x : A i x ≤ b i } and H + i (b) = {x : A i x ≥ b i } the associated half-spaces. Then Λ b = i∈[1,m] H - i (b). Now fix b 0 ∈ B and consider v ∈ V (b 0 )
. By definition, v is the intersection of hyperplanes i∈K H i (b 0 ) where K ⊂ {1, . . . , m} is maximal (note that k := card(K) ≥ n otherwise v can not be a vertex).

For J ∈ {1, . . . , m}, let A J denote the submatrix of A obtained by considering the rows A j for j ∈ J. Let us introduce the sets of derived points (points of the arangement) of the set K, for each b ∈ B:

V K (b) := {x ∈ R n ; ∃J ⊂ K ; A J is invertible and x = A -1 J b} . By definition, V K (b 0 ) = {v} and, for each b ∈ B, V K (b) is a set of at most k n elements. First, note that for each b ∈ B and v := A -1 J b ∈ V K (b), one has: v -v = A -1 J b 0 -A -1 J b ≤ A -1 J × b 0 -b ≤ α b 0 -b (68)
where α := max with y θ = (0, E θ ) the preference of user θ for period P and ω θ := θ the preference weight of player θ.

We consider approximating atomic games by splitting Θ uniformly (Sec. 3.4.1) in 5, 20, 40 and 100 segments (players). We compute the NE for each atomic game using the best-response dynamics (each best-response is computed as a QP using algorithm [START_REF] Brucker | An o(n) algorithm for quadratic knapsack problems[END_REF], see [START_REF] Jacquot | Analysis and implementation of an hourly billing mechanism for demand response management[END_REF] for convergence properties) and until the KKT optimality conditions for each player are satisfied up to an absolute error of 10 -3 . Fig. 2 shows, for each NE associated to the atomic games with 5, 20, 40 and 100 players, the linear interpolation of the load on the peak period x θ,P (red filled area), while the load on the offpeak period can be observed as x θ,O = E θ -x θ,P . We observe the convergence Figure 2: Convergence of the Nash Equilibrium profiles to a Wardrop Equilibrium profile to the limit WE of the nonatomic game as stated in Thm. 14. We also observe that the only discontinuity point of θ → x * θ,P comes from the discontinuity of θ → E θ at θ = 0.7, as stated in Thm. 13.

Conclusion

This paper gives quantitative results on the convergence of Nash equilibria, associated to atomic games approximating a nonatomic routing game, to the Wardrop equilibrium of this nonatomic game. These results are obtained under different differentiability and monotonicity assumptions. Several directions can be explored to continue this work: first, we could analyze how the given theorems could be modified to apply in case of nonmonotone and nondifferentiable functions. Another natural extension would be to consider routing games on nonparallel networks or even general aggregative games: in that case, the separable costs structure is lost and the extension is therefore not trivial.

Figure 1 :

 1 Figure 1: A parallel network with T links.

  η := min j∈{1...m}\K d v, H + j (b) . By the maximality of K, η > 0. As x → d(x, H + j ) is continuous for each j, and from (68) , there exists δ > 0 such that: b0 -b ≤ δ =⇒ ∀v ∈ V K (b), min j∈{1...m}\K d v , H + j (b) > 0. Next, we show that, for b such that b 0 -b ≤ δ, there exists v ∈ V K (b) ∩ V (b). We proceed by induction on k -n. If k = n, then v = A -1 K b 0 and for any b in the ball S δ (b 0 ), V K (b) = {A -1 K b}. Thus v := A -1 K b verifies A K v = b K , and A j v < b j for all j / ∈ K, thus v belongs to V (b). If k = n + t with t ≥ 1, there exists j 0 ∈ K such that with K = K \ {j 0 }, V K (b 0 ) = {v}. Consider the polyhedron P = i∈K H - i (b 0 ). By induction, there exists J ⊂ K such that A -1 J b is a vertex of P . If it satisfies also A j0 x ≤ b j0 then it is an element of V (b). Else, consider a vertex v of the polyhedron P ∩ H - j0 (b) on the facet associated with H j0 (b). Then, v ∈ V K (b) and, as b ∈ S δ (b 0 ), it verifies A j v < b j for all j / ∈ K, thus v ∈ V (b). Thus,in any case and for b ∈ S δ (b 0 ), d(v, V (b)) ≤ vv ≤ α b 0 -b and finally d (V (b 0 ), V (b)) ≤ α b 0 -b . The collection S δ b (b) b∈B is an open covering of the compact set B, thus there exists a finite subcollection of cardinal r that also covers B, from which we deduce that there exists D ≤ max(rα) such that: ∀b, b ∈ B, d H (V (b ), V (b)) ≤ D bb .

Acknowledgments

We thank Stéphane Gaubert, Marco Mazzola, Olivier Beaude and Nadia Oudjane for their insightful comments.

This work was supported in part by the PGMO foundation.

Finite dimension, meshgrid approximation

Consider a nonatomic routing game G = (Θ, X , F ) (Def. 3) satisfying the following two hypothesis:

• The feasibility sets are K-dimensional polytopes: there exist A ∈ M K,T (R) and b : Θ → R K bounded, such that for any θ, X θ := {x ∈ R T ; Ax ≤ b θ }, with X θ nonempty and bounded (as a polytope, X θ is closed and convex).

• There exist a bounded function s : Θ → R q and a function u : R q × B 0 (M ) → R such that for any θ ∈ Θ, u θ = u(s θ , .). Furthermore, u is Lipschitz-continuous in s. For ν ∈ N * , we consider the uniform meshgrid of ν K+q classes of

Let

[s k , s k ] which will give us a set of I (ν) = ν K+q subsets. More explicitly, if we define:

the set of indices for the meshgrid, and with the cutting points b k,

n of Θ as:

Since some of the subsets Θ

n can be of Lebesgue measure 0, we define the set of players I (ν) as the elements n of Γ (ν) for which µ(Θ

Remark 3. If there is a set of players of positive measure that have equal parameters b and s, then the condition max i∈I (ν) µ i → 0 will not be satisfied. In that case, adding another dimension in the meshgrid by cutting Θ = [0, 1] in ν uniform segments solves the problem. Proposition 16. For ν ∈ N * , consider the atomic game G (ν) defined by:

Before giving the proof of Thm. 16, we show the following nontrivial Thm. 17, from which the convergence of the feasibility sets is easily derived.

Proof of Thm. 16. First, to show the divergence of the number of players and their infinitesimal weight, we have to follow Remark 3 where we consider an additional splitting of the segment Θ = [0, 1]. In that case, we will have I (ν) ≥ ν hence goes to positive infinity and for each n ∈ I (ν) , µ(Θ

ν hence goes to 0. Then, the convergence of the strategy sets follows from the fact that, for each n ∈ I (ν) :

and from Thm. 17 which implies that, for each θ ∈ Θ (ν)

n :

Finally, the convergence of utility functions comes from the Lipschitz continuity in s. For each n ∈ I (ν) and each θ ∈ Θ (ν) n , we have:

which terminates the proof. Of course, the number of players considered in Thm. 16 is exponential in the dimensions of the parameters K + q, which can be large in practice. As a result, the number of players in the approximating atomic games considered can be very large, which will make the NE computation really long. Taking advantage of the continuity of the parametering functions and following the approach of Thm. 15 gives a smaller (in terms of number of players) approximating atomic instance.

Numerical Application

We consider a population of consumers Θ = [0, 1] with an energy demand distribution θ → E θ . Each consumer θ splits her demand over T := {O, P }, so that her feasibility set is X θ := {(x θ,O , x θ,P ) ∈ R 2 + x θ,O + x θ,P = E θ }. The index O stands for offpeak-hours with a lower price c O (X) = X and P are peak-hours with a higher price c P (X) = 1 + 2X. The energy demand and the utility function in the nonatomic game are chosen as the piecewise continuous functions: