Xiaoxin Lu

Fabrice Detrez
email: fabrice.detrez@univ-eiffel.fr

Julien Yvonnet

Jinbo Bai

Identification of elastic properties of interphase and interface in graphene-polymer nanocomposites by atomistic simulations

Keywords: Polymer nanocomposites, Graphene, Interface, Interphase, Mechanical properties

come

Identification of elastic properties of interphase and interface in graphene-polymer nanocomposites by atomistic simulations

Introduction

Graphene, as one of the stiffest known materials, whose Young's modulus is around 1 TPa [START_REF] Van Lier | Ab initio study of the elastic properties of single-walled carbon nanotubes and graphene[END_REF]2], is an attractive candidate for use in high-performance polymer-based nanocomposites. From 2000s, a significant amount of research has been carried out on graphene nano-platelets (GNP) or graphene-derived materials based polymer nanocomposites due to their potential for large increases in toughness, strength and stiffness [3,4,5,[START_REF] Park | Epoxy toughening with low graphene loading[END_REF]. The enhanced mechanical properties of the nanocomposites have been found to depend on various factors including the GNP-polymer bonding, the concentration and dispersion of GNP, the defects in GNP [START_REF] Zhao | Enhanced mechanical properties of graphene-based poly (vinyl alcohol) composites[END_REF][START_REF] Terrones | Interphases in graphene polymer-based nanocomposites: achievements and challenges[END_REF][START_REF] Tang | The effect of graphene dispersion on the mechanical properties of graphene/epoxy composites[END_REF]. Many experimental studies have shown that the GNP polymer interface and the interphase region greatly influence the effective mechanical properties of the nanocomposites (see, e.g. [START_REF] Zaman | Epoxy/graphene platelets nanocomposites with two levels of interface strength[END_REF]11,[START_REF] Ma | Covalently bonded interfaces for polymer/graphene composites[END_REF]). In this paper, the polymer close to GNP filler is called interphase and the term interface indicates a discontinuity surface, which is generally located between two materials.

The existence and the influence of interfaces and interphases have been experimentally observed in GNP polymer nanocomposites [START_REF] Gong | Interfacial stress transfer in a graphene monolayer nanocomposite[END_REF][START_REF] Zhang | The effect of surface chemistry of graphene on rheological and electrical properties of polymethylmethacrylate composites[END_REF][START_REF] Liu | Quantitative study of interface/interphase in epoxy/graphene-based nanocomposites by combining STEM and EELS[END_REF]. They have been theoretically demonstrated in polymer nanocomposites by atomistic simulations [START_REF] Brown | A molecular dynamics study of a model nanoparticle embedded in a polymer matrix[END_REF][START_REF] Li | Atomistic simulations on multilayer graphene reinforced epoxy composites[END_REF][START_REF] Rahman | Molecular modeling of crosslinked graphene-epoxy nanocomposites for characterization of elastic constants and interfacial properties[END_REF][START_REF] Rissanou | Dynamics of various polymer-graphene interfacial systems through atomistic molecular dynamics simulations[END_REF][START_REF] Moon | Interfacial strengthening between graphene and polymer through Stone-Thrower-Wales defects: Ab initio and molecular dynamics simulations[END_REF]. However, their respective roles on elastic behavior are not yet well understood. So, it is necessary to have an identification procedure which is able to decouple the elastic contributions of interfaces and interphases. This is hardly accessible through an inverse analysis based on atomistic simulations or experimental data because the effects of interfaces and interphases can be hardly distinguished. The objective of this paper is to propose an identification methodology of elastic models with interphases and/or imperfect interfaces in GNP polymer nanocomposites. Possible applications of these continuum models could be the definition of a finite element framework taking into account the size effects to study at mesoscale the concentration and dispersion of GNP on elastic properties (see, e.g. [START_REF] Mortazavi | Interphase effect on the elastic and thermal conductivity response of polymer nanocomposite materials: 3D finite element study[END_REF][START_REF] Mortazavi | Multiscale modeling of thermal and mechanical properties of nanostructured materials and polymer nanocomposites[END_REF][START_REF] Lu | Multiscale modeling of nonlinear electric conductivity in graphene-reinforced nanocomposites taking into account tunnelling effect[END_REF][START_REF] Lu | Multiscale electro-mechanical modeling of graphene/polymer nanocomposites[END_REF][START_REF] Lu | Low electrical percolation thresholds and nonlinear effects in graphene-reinforced nanocomposites: a numerical analysis[END_REF] for similar studies on electric and thermal conductivity properties). We will restrict our study to a sandwich structure model of nanocomposite (see Figure 2) for reasons of clarity.

The strong influence of interfaces and interphases induces size effects on the mechanical behavior. These size effects have been observed in several polymer nanocomposites (see, e.g. [START_REF] Tjong | Structural and mechanical properties of polymer nanocomposites[END_REF][START_REF] Crosby | Polymer nanocomposites: the "nano" effect on mechanical properties[END_REF] for a review). There are mainly two families of continuum models in the literature able to capture size effects: the models with an interphase surrounding inclusions and the imperfect interface models. The first class of models with interphases have been used on polymer nanocomposites for mean-field homogenization (see, e.g. [START_REF] Peng | Modeling of nano-reinforced polymer composites: Microstructure effect on Young's modulus[END_REF][START_REF] Liu | An extended micromechanics method for probing interphase properties in polymer nanocomposites[END_REF]) and for numerical homogenization [START_REF] Mortazavi | Interphase effect on the elastic and thermal conductivity response of polymer nanocomposite materials: 3D finite element study[END_REF]. The mechanical studies using an interphase region in polymer carbon nanocomposites were mostly performed on the polymer reinforced with carbon nanotubes [START_REF] Wan | Effect of CNT length and CNT-matrix interphase in carbon nanotube (CNT) reinforced composites[END_REF][START_REF] Seidel | Micromechanical analysis of the effective elastic properties of carbon nanotube reinforced composites[END_REF][START_REF] Hernández-Pérez | Modeling the influence of interphase on the elastic properties of carbon nanotube composites[END_REF][START_REF] Yang | Nonlinear multiscale modeling approach to characterize elastoplastic behavior of CNT/polymer nanocomposites considering the interphase and interfacial imperfection[END_REF][START_REF] Srivastava | A continuum model to study interphase effects on elastic properties of CNT/GSnanocomposite[END_REF][START_REF] Tserpes | Prediction of yield strength of MWCNT/PP nanocomposite considering the interphase and agglomeration[END_REF]. Atomistic simulations are also performed to identify the interphase elastic parameter by inverse analysis to reproduce the size effect [START_REF] Odegard | Modeling of the mechanical properties of nanoparticle/polymer composites[END_REF][START_REF] Marcadon | Confrontation between molecular dynamics and micromechanical approaches to investigate particle size effects on the mechanical behaviour of polymer nanocomposites[END_REF][START_REF] Le | Stochastic continuum modeling of random interphases from atomistic simulations. Application to a polymer nanocomposite[END_REF][START_REF] Kim | Multiscale modeling of interphase in crosslinked epoxy nanocomposites[END_REF]. In these studies, the interphase is assumed isotropic and the elastic parameters are fited to reproduce the experimental or atomistic simulation data.

The study of interphase anisotropy is very delicate through inverse analysis. The second class of models that incorporates zero-thickness imperfect interfaces [START_REF] Bövik | On the modelling of thin interface layers in elastic and acoustic scattering problems[END_REF][START_REF] Gurtin | A general theory of curved deformable interfaces in solids at equilibrium[END_REF][START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces[END_REF][START_REF] Benveniste | Models of thin interphases with variable moduli in plane-strain elasticity[END_REF] is also able to capture the size effects [START_REF] Gu | Size-dependent effective elastic moduli of particulate composites with interfacial displacement and traction discontinuities[END_REF][START_REF] Chatzigeorgiou | Generalized interfacial energy and size effects in composites[END_REF]. The general imperfect interfaces combine the discontinuity of displacement of cohesive zone model [START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture. general ideas and hypotheses. Axiallysymmetric cracks[END_REF][START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and the discontinuity of the traction vector of the membrane-type interface model [START_REF] Gurtin | A continuum theory of elastic material surfaces[END_REF][START_REF] Murdoch | A thermodynamical theory of elastic material interfaces[END_REF]. To the best of our knowledge, although the atomistic simulations are widely used to identity the physical effective parameters of the materials such as Young's modulus, glass transition parameter, etc. [START_REF] Skountzos | Determination of the mechanical properties of a poly(methyl methacrylate) nanocomposite with functionalized graphene sheets through detailed atomistic simulations[END_REF][START_REF] Fan | Molecular dynamics predictions of thermomechanical properties of an epoxy thermosetting polymer[END_REF], the atomistic identification of the anisotropic local materials parameters remains an open question. Nevertheless, special cases of imperfect interfaces have identified by atomistic simulations, such as membrane-type interface in thin material [START_REF] Shenoy | Atomistic calculations of elastic properties of metallic fcc crystal surfaces[END_REF][START_REF] Park | A surface Cauchy-Born model for nanoscale materials[END_REF][START_REF] Yvonnet | Finite element model of ionic nanowires with size-dependent mechanical properties determined by ab initio calculations[END_REF][START_REF] Yvonnet | Characterization of surface and nonlinear elasticity in wurtzite ZnO nanowires[END_REF][START_REF] Davydov | On molecular statics and surface-enhanced continuum modeling of nanostructures[END_REF][START_REF] Davydov | A comparison of atomistic and surface enhanced continuum approaches at finite temperature[END_REF][START_REF] Hoang | Size-dependent mechanical properties of axial and radial mixed AlN/GaN nanostructure[END_REF][START_REF] Elsner | Surface excess elasticity of gold: Ab initio coefficients and impact on the effective elastic response of nanowires[END_REF]; and cohesive zone model to the study of the decohesion between filler and matrix in polymer nanocompoistes [START_REF] Lu | Multiscale electro-mechanical modeling of graphene/polymer nanocomposites[END_REF][START_REF] Namilae | Multiscale model to study the effect of interfaces in carbon nanotube-based composites[END_REF][START_REF] Lu | A cohesive law for multi-wall carbon nanotubes[END_REF][START_REF] Awasthi | Modeling of graphene-polymer interfacial mechanical behavior using molecular dynamics[END_REF][START_REF] Liu | A boundary element method for the analysis of CNT/polymer composites with a cohesive interface model based on molecular dynamics[END_REF][START_REF] Chen | Mechanical characterization of interfaces in epoxy-clay nanocomposites by molecular simulations[END_REF][START_REF] Song | Multiscale modeling of damage progression in nylon 6/clay nanocomposites[END_REF][START_REF] Paliwal | Nanomechanical modeling of interfaces of polyvinyl alcohol (PVA)/clay nanocomposite[END_REF]. Note that some authors performed atomistic simulations to study the cohesion between carbon filler and polymer matrix without identifying the cohesive zone model [START_REF] Li | Atomistic simulations on multilayer graphene reinforced epoxy composites[END_REF][START_REF] Rahman | Molecular modeling of crosslinked graphene-epoxy nanocomposites for characterization of elastic constants and interfacial properties[END_REF][START_REF] Liao | Interfacial characteristics of a carbon nanotube-polystyrene composite system[END_REF][START_REF] Gou | Computational and experimental study of interfacial bonding of single-walled nanotube reinforced composites[END_REF][START_REF] Jiang | A cohesive law for carbon nanotube/polymer interfaces based on the van der Waals force[END_REF][START_REF] Chowdhury | Computer simulation of carbon nanotube pull-out from polymer by the molecular dynamics method[END_REF].

We introduce the methodology of Atomistic Local IdentificAtion of Stiffness, called in this paper ALIAS methodology, and which allows to establish an equivalent continuum model from an atomistic model. It consists of measuring the local stress and strain fields under different deformations of the simulation box and then deducing the elastic tensor fields. The novel possibilities offered by the ALIAS methodology are the identification of a continuum model allowing to take into account size effects from an unique size box of atomistic model, and also the characterization of interphase anisotropy. The ALIAS methodology is a sequential multiscale approach following the classification proposed in the comprehensive book of [START_REF] Tadmor | Modeling materials: continuum, atomistic and multiscale techniques[END_REF]. The other category of atomistic-continuum coupling approaches is the concurrent method (see, e.g. [START_REF] Tadmor | Modeling materials: continuum, atomistic and multiscale techniques[END_REF][START_REF] Li | Handbook of micromechanics and nanomechanics[END_REF] for an overview) and the application to polymer material [START_REF] Pfaller | An Arlequin-based method to couple molecular dynamics and finite element simulations of amorphous polymers and nanocomposites[END_REF][START_REF] Pfaller | Investigation of interphase effects in silica-polystyrene nanocomposites based on a hybrid molecular-dynamics-finite-element simulation framework[END_REF][START_REF] Pfaller | Optimisation of the Capriccio method to couple particle-and continuum-based simulations of polymers[END_REF]. Although they contain more physical ingredients, the main drawback for GNP nanocomposite applications is the GNP number in a Representative Volume Element (RVE) that induces a large number of atomistic domains, leading to untractable computational times. However, the ALIAS method is inspired by concurrent methods in the sense that they are based on the same two key features: the collection of continuum information from atomic level and the constraint of atomistic simulations to impose the continuum field. For the first point, the Murdoch-Hardy procedure [START_REF] Murdoch | Physical Foundations of Continuum Mechanics[END_REF] is used because it allows to define both the velocity field and the stress field at all points.

Regarding the second point, the Cauchy-Born rule [START_REF] Born | On the stability of crystal lattices I[END_REF][START_REF] Weiner | Hellmann-Feynman theorem, elastic moduli, and the Cauchy relations[END_REF] is the fundamental assumption to link the continuum strain to the atom displacement. It is suitable for the multiscale coupling of bulk materials such as crystalline materials [see, e.g. [START_REF] Arroyo | An atomistic-based finite deformation membrane for single layer crystalline films[END_REF][START_REF] Arroyo | Finite crystal elasticity of carbon nanotubes based on the exponential Cauchy-Born rule[END_REF]. However, for the amorphous materials such as glass or polymer, it is no longer applicable [START_REF] Tanguy | Continuum limit of amorphous elastic bodies: a finite-size study of low-frequency harmonic vibrations[END_REF][START_REF] Wittmer | Vibrations of amorphous, nanometric structures: When does continuum theory apply?[END_REF][START_REF] Maloney | Universal breakdown of elasticity at the onset of material failure[END_REF][START_REF] Maloney | Amorphous systems in athermal, quasistatic shear[END_REF][START_REF] Lemaître | Sum rules for the quasi-static and visco-elastic response of disordered solids at zero temperature[END_REF]. To overcome this difficulty, we impose the deformation to the atomic simulation box and we compute the atoms displacement by the Athermal, Quasistatic Simulations algorithm introduced by [START_REF] Maeda | Computer simulation of deformation in two-dimensional amorphous structures[END_REF][START_REF] Kobayashi | Computer simulation of deformation of amorphous Cu57Zr43[END_REF][START_REF] Maeda | Atomistic process of plastic deformation in a model amorphous metal[END_REF] and used to determined the elastic constants of amorphous glasses materials by [START_REF] Tanguy | Continuum limit of amorphous elastic bodies: a finite-size study of low-frequency harmonic vibrations[END_REF][START_REF] Wittmer | Vibrations of amorphous, nanometric structures: When does continuum theory apply?[END_REF][START_REF] Maloney | Universal breakdown of elasticity at the onset of material failure[END_REF][START_REF] Maloney | Amorphous systems in athermal, quasistatic shear[END_REF][START_REF] Lemaître | Sum rules for the quasi-static and visco-elastic response of disordered solids at zero temperature[END_REF]. The novelty of ALIAS method is the estimation of the local displacement field during simulations for a given macroscopic strain imposed to the atomistic simulation box. This, allows us to deduce the strain fields and the stiffness tensor at each point.

The paper is organized as follows. First, in section 2, a short overview of the general imperfect interface framework is established. In section 3, we describe the atomistic model, the ALIAS methodology and their results. Finally, in section 4, we discuss the influence of interphase and interface on elastic properties of sandwich graphene/polymer nano-composite.

Notations and definitions

Notations

A : B = (A) ij (B) ij .
Quantities defined on the interface are distinguished from those in the bulk by a superscript (•) s . The continuum fields identified by the ALIAS methodology are designated by a sub-script (•) w Moreover, macroscale quantities and effective quantities are differentiated from microscale quantities by a bar placed above the quantity ( •). The jump of a quantity (•) over the interface I are defined by

• = •| + I -•| - I .
The list of symbols used in this paper is provided in Table 1. Atomistic mass density at point x u

(E) w (x), v w (x)
Atomistic Eulerian displacement and atomistic velocity fields at point x σ w (x), w (x) Atomistic Cauchy stress tensor and atomistic strain tensor at point x C w (x)

Atomistic stiffness tensor at point x s w (x)

Nematic order parameter at point x

Generalized imperfect interface framework

The purpose of this section is to establish the equations governing continua embedding general interface elastic model. A comprehensive description of this model can be found for the case of infinitesimal displacements in [START_REF] Chatzigeorgiou | Generalized interfacial energy and size effects in composites[END_REF] and for its extension to finite deformation in [START_REF] Javili | Micro-to-macro transition accounting for general imperfect interfaces[END_REF][START_REF] Javili | A note on traction continuity across an interface in a geometrically non-linear framework[END_REF]. In this section, we assume infinitesimal displacements and linear behavior for the sake of simplicity. The interface I splits the body Ω into two regions (see Figure 1). The unit normal of the interface I is denoted by n I .

The displacement in the body is denoted by u (x) and the displacement jump, u , across the interface

I is defined by u = u| + I -u| - I , (1)
where u| The displacement field, u s , inside the interface I is defined as the median value according to:

u s (x) = 1 2 u| + I + u| - I , ∀x ∈ I. (2)
The infinitesimal strain tensor in the bulk is defined by

(x) = 1 2 ∇u (x) + ∇u T (x) , ∀x ∈ Ω, (3)
and the infinitesimal strain tensor on the interface is defined by

s (x) = 1 2 ∇ s u s (x) + (∇ s u s (x)) T , ∀x ∈ I, (4)
where ∇ is the gradient operator, ∇ s = P • ∇ is the interface gradient operator and P = In I ⊗ n I is the projector onto the tangent plane of the interface I.

In absence of external force densities in bulk and in the interface, the balance equations are given by:

∇ • σ = 0 ∀x ∈ Ω, (5)
∇ s • σ s + t = 0 ∀x ∈ I, (6)
where ∇ • is the divergence operator;

∇ s • { • } = ∇ • { • } : P is the interface divergence operator; σ s
the interfacial stress and t is the traction vector jump across the interface defined by

t = σ • n I = σ| + I -σ| - I • n I . ( 7
)
We assume the existence of a bulk free energy density φ( ) and an interfacial free energy density φ s ( s , u ) in the bulk and the interface respectively, such as the behavior laws derive from them

σ = ∂ φ ∀x ∈ Ω, (8)
σ s = ∂ s φ s ∀x ∈ I, (9)
t s = ∂ u φ s ∀x ∈ I, (10)
where t s denotes the average traction across the interface, defined by

t s = 1 2 σ| + I + σ| - I • n I . ( 11
)
Under the assumptions of standard linear elastic behavior with internal stresses, both free energies are given by

φ ( ) = 1 2 : C (x) : + τ (x) : + φ 0 ∀x ∈ Ω, ( 12
)
φ s ( s , u ) = 1 2 s : C s (x) : s + τ s (x) : s + 1 2 u .K s (x) . u + φ s 0 ∀x ∈ I, (13)
where C (x) is the fourth-order symmetric bulk stiffness tensor; C s (x) is the fourth-order symmetric interfacial stiffness tensors; τ is the internal stresses in the bulk and τ s is the internal surface stress; φ 0 and φ s 0 are the free energies of bulk and interface under zero strain and zero jump displacement; and K s is the second-order symmetric tensor which corresponds to the cohesive stiffness of interface. c). The initial system is prepared from the off-lattice self-avoiding random walk combining the simultaneously molecular dynamics relaxation [START_REF] Lu | Multiscale study of influence of interfacial decohesion on piezoresistivity of graphene/polymer nanocomposites[END_REF]. We assume that the particles interact through the DREIDING potential energy function [START_REF] Mayo | DREIDING: a generic force field for molecular simulations[END_REF], where the electrostatic contributions are neglected for simplicity reasons. The potential energy U ({r α } , ¯ ), depends on the current position of particles, {r α }, via the DREIDING potential as well as on the shape of simulation box via the effective infinitesimal strain tensor ¯ of the box defined by

Identification of continuum elastic parameters by atomistic simulations

¯ = 1 2 F + F T -I, where F = h • h (0) -1 , (14)
I is the unit second order tensor and h (0) defines the shape of the parallelepiped simulation box in the initial reference configuration (¯ = 0).

In our study, we suppose that the tensor h (t) is parameterized by a scalar t called pseudo time.

We assume that the system follows an equilibrium trajectory, even though the actual potential energy

U ({r α } , ¯ (t)
) is in a local minimum. The Athermal, Quasistatic Simulations algorithm is used to ensure that the system follows an equilibrium trajectory. This simulation technics have been introduced by [START_REF] Maeda | Computer simulation of deformation in two-dimensional amorphous structures[END_REF][START_REF] Kobayashi | Computer simulation of deformation of amorphous Cu57Zr43[END_REF][START_REF] Maeda | Atomistic process of plastic deformation in a model amorphous metal[END_REF] and used to determine the elastic constants of amorphous glasses materials by [START_REF] Tanguy | Continuum limit of amorphous elastic bodies: a finite-size study of low-frequency harmonic vibrations[END_REF][START_REF] Wittmer | Vibrations of amorphous, nanometric structures: When does continuum theory apply?[END_REF][START_REF] Maloney | Universal breakdown of elasticity at the onset of material failure[END_REF][START_REF] Maloney | Amorphous systems in athermal, quasistatic shear[END_REF][START_REF] Lemaître | Sum rules for the quasi-static and visco-elastic response of disordered solids at zero temperature[END_REF].

For each small pseudo time increment, δt, this algorithm consists of moving affinely the atoms from previous state, and then to perform a minimization to find the nearest local minimum.

The velocity of α-th atom is estimated by backward finite difference approximation

v α (t) ≈ r α (t) -r α (t -δt) δt . ( 15
)
It is interesting to note that the velocity tends to zero in the quasistatic limit (i.e. v α → 0 and δt → ∞), therefore we could neglect it in most of cases. However, the atom velocity atomic is necessary to define the displacement field.

The equilibrium trajectory imposes that the force on particle α is always zero

f α = -∂ rα U ({r α } , ¯ (t)) = β α f αβ = 0. ( 16
)
Here f αβ is the force of β-th atom on the α-th atom. Since the potentials used are conservative, they can be expressed in terms of set of inter-atomic distances r αβ between the α-th and β-th atoms [START_REF] Tadmor | Modeling materials: continuum, atomistic and multiscale techniques[END_REF][START_REF] Admal | A unified interpretation of stress in molecular systems[END_REF].

Consequently, a set of central forces between the pairs of atoms f αβ can be defined by

f αβ = -∂ r αβ U ({r α } , ¯ (t)) r α -r β r αβ . ( 17
)
The expressions of DREIDING potential in term of interatomic distances is given in supplementary information. It is interesting to note that this definition ensures that

f αβ = -f βα , and f α = β α f αβ . ( 18
)

Atomistic Local IdentificAtion of Stiffness: ALIAS methodology

The ALIAS methodology is based on the measurement of continuum stress and strain tensors for several imposed effective infinitesimal strain tensors ¯ . The number of imposed effective infinitesimal strain tensors is chosen to explore a necessary set of strain-stress combinations to defined the stiffness fields. The Murdoch-

Hardy procedure [START_REF] Murdoch | Physical Foundations of Continuum Mechanics[END_REF] is used to measure the continuum velocity and Cauchy's stress tensor fields from atomistic quantities. In this paper, the continuum displacement field that allows to compute the strain field is obtained by integration of the velocity field.

Murdoch-Hardy procedure

The main point here is to identify the Eulerian specification of continuum fields, such as mass density, Cauchy's stress tensor and velocity, with local space averages of atomistic quantities, such as atom mass, atom velocity and interatomic forces. A physical interpretation of the Murdoch-Hardy procedure is that a device with a finite size measures the conservative quantities of interest at a point x. This device is modeled by a scalar-valued weighting function w(x) which defines a spatial averaging and satisfies the normalization condition

R 3 w (x) dx = 1. ( 19
)
The mass density ρ w (x, t) is defined as

ρ w (x, t) = α m α w (r α (t) -x) . ( 20
)
Here, the mass density ρ w (x, t) must verify the continuity equation (conservation of mass) in the current configuration

∂ t ρ w + ∇ x • (ρ w v w ) = 0, (21)
where ∂ t and ∇ x • denote the partial derivative with respect to the time and the divergence operator with respect to the spatial coordinate x, respectively.

To satisfy the continuity equation ( 21), the Eulerian specification of velocity fields v w (x, t) is defined by

[see e.g. 76] v w (x, t) = 1 ρw α m α v α (t) w (r α (t) -x) if ρ w (x, t) 0 0 otherwise . ( 22
)
Following [START_REF] Admal | A unified interpretation of stress in molecular systems[END_REF][START_REF] Admal | The non-uniqueness of the atomistic stress tensor and its relationship to the generalized beltrami representation[END_REF][START_REF] Admal | Material fields in atomistics as pull-backs of spatial distributions[END_REF], the atomistic Cauchy stress tensor is decomposed in two terms: the kinematics contribution, σ w,k , and the potential contribution, σ w,v , such as

σ w (x, t) = σ w,k (x, t) + σ w,v (x, t) , ( 23
)
where

σ w,k (x, t) = - α (v α (t) -v w (x, t)) ⊗ (v α (t) -v w (x, t)) w (r α (t) -x) . ( 24
)
Here, the kinematics contribution of the atomistic Cauchy stress tensor and the velocity fields are assumed to be negligible under the quasistatic assumption (δt → ∞), i.e. v α ∝ (δt) -1 → 0 for all α, t implying

v w (x, t) ∝ (δt) -1 → 0 and σ w,k (x, t) ∝ (δt) -2 → 0 for all x, t.
The potential contribution of atomistic Cauchy stress, σ w,v satisfies the balance equation

∇ x .σ w,v (x, t) = α f α (t) w (r α (t) -x) . ( 25
)
The expression of σ w,v is not unique, the most common one is

σ w,v (x, t) = 1 2 α,β f αβ (t) ⊗ (r β (t) -r α (t)) b w (x; r α (t) , r β (t)) , ( 26
)
where b w (x, r α (t) , r β (t)) is the bond function expressed by

b w (x; r α (t) , r β (t)) = 1 0 w((1 -s)r α + sr β -x)ds. (27)
Note that the atomisitic definition of stress is not unique [START_REF] Tadmor | Modeling materials: continuum, atomistic and multiscale techniques[END_REF][START_REF] Admal | A unified interpretation of stress in molecular systems[END_REF][START_REF] Admal | The non-uniqueness of the atomistic stress tensor and its relationship to the generalized beltrami representation[END_REF] and is a controversial topic in the literature [START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF].

Methotology of stiffness identification

First, we introduce a definition of the infinitesimal strain tensor field from the Mudoch-Hardy velocity field v w (x, t) through the concept of the Eulerian displacement field. Following [START_REF] Gremaud | Eulerian theory of newtonian deformable lattices-dislocation and disclination charges in solids[END_REF], we introduce the Eulerian displacement field, u

(E)
w (x, t), which allows to reconstruct the solid body as it was at the initial instant in an absolute referential before all transformations. The relation between the Eulerian displacement field and velocity fields is

v w (x, t) = - du (E) w dt , or ∂ t u (E) w = -v w . I + ∇ x u (E) w , ( 28
)
where ∇ x is the gradient operator with respect to the spatial coordinate x. The Eulerian displacement field can be estimated by a finite difference approximation

u (E) w (x, t + δt) ≈ u (E) w (x, t) -δtv w (x, t + δt) • I + ∇ x u (E) w (x, t) . ( 29
)
By introducing the definitions of the Eulerian specification of velocity fields v w (x, t) (Eq. 22) combined with the definition of atomic velocities v α (t) (Eq. 15), we can write the increment of the Eulerian displacement field :

∆u (E) w (x, t) = u (E) w (x, t + δt) -u (E) w (x, t) = - δt ρ w α m α v α (t) w (r α (t) -x) • I + ∇ x u (E) w (x, t) (30) = - 1 ρ w α m α (r α (t) -r α (t -δt)) w (r α (t) -x) • I + ∇ x u (E) w (x, t) .
Note that the increment of Eulerian displacement field is independent of speudo time δt and that it is well defined in the quasistatic limit (δt → ∞).

Under infinitesimal displacement assumption, the initial and final configuration can be considered as the same, therefore the classical displacement field u w can be defined as the opposite of the Eulerian displacement

field u w = -u (E)
w (x). The infinitesimal strain tensor field w is defined by:

( w ) ij (x) = 1 2 ∂ xi (u w ) j (x) + ∂ xj (u w ) i (x) . ( 31
)
Here, we use a centered finite difference approximation of the spatial derivative of displacement field

∂ xi (u w ) j (x) ≈ (u w ) j (x + δxe i ) -(u w ) j (x -δxe i ) 2δx . ( 32
)
The atomistic stiffness tensor C w (x) at each point where the atomistic Cauchy's stress σ w and the atomistic infinitesimal strain w are measured is defined by the least squares method, for a given set of K effective strain tensors ¯ (k) , such as

C w (x) = argmin Cw(x)∈Ela K k=1 σ (k) w (x) -C w (x) : (k) w (x) -τ w (x) 2 F ( 33
)
where Ela is the set of fourth order elastic tensor with minor and major symmetry and

a 2 F = a ij a ij is the
Frobenius norm associated to 2nd order tensor. Here σ

(k) w (x) (resp. (k) w (x))
is the atomistic Cauchy's stress tensor (resp. the atomistic infinitesimal strain) that is measured at x point for the effective strain tensors of the simulation box ¯ = ¯ (k) . τ w (x) is the internal Cauchy's stress tensor at x point that is computed on reference configuration (¯ = 0). Note that

(k)
w (x) is null by definition for all x on the reference configuration.

Application of sandwich structure of graphene polymer nanocomposite

The origin of absolute referential, R Euler (O, e 1 , e 2 , e 3 ), in the Eulerian specification is the center of the graphene sheet and its orientation is defined such that e 3 is normal to the graphene and the triangles in the graphene hexagonal lattice pointing along the e 1 direction (see Figure 2). x, ( Å) Here, we want to probe the mechanical properties as a function of the z distance from the graphene sheet. For that reason an anisotropic measuring device with rectangular cuboid shape is used such with lengths l w 1 = L x , l w 2 = L y and l w 3 = 4Å. The length l w 3 is chosen to be greater than the C-C bond and small enough to avoid considering graphene atoms when we estimate the interphase properties. We assume that the measure is uniform for all points in the cuboid and that the atoms occupy a small sphere of radius r w = 0.77Å to regularize the weighting function w(x). The atomic radius r w is chosen as half of the C-C bond of the DREIDING potential [START_REF] Mayo | DREIDING: a generic force field for molecular simulations[END_REF]. Under these assumptions, the weighting function is defined by a convolution product

ω (x, 0, 0) × ω 0 -4 -2 0 2 4 0 0.5 1 z, ( Å) ω (0, 0, z) × ω 0 (a) (b)
w (x) = 1 w 0 Ω f (Cub) (x)f (Sph) (x -x )dx , ( 34
)
where w 0 is the normalization constant given by equation [START_REF] Rissanou | Dynamics of various polymer-graphene interfacial systems through atomistic molecular dynamics simulations[END_REF] and f (Cub) and f (Sph) are the characteristic function of the cuboid and sphere respectively given by

f (Sph) (x) = 1 -H (|x| -r w ) , ( 35
)
f (Cub) (x) = 3 i=1 1 -H x.e i - l w i 2 , ( 36
)
where H (x) denotes the Heaviside step function. The representation of weighting function along e 1 and e 3 is shown in Figure 3.

Elongations and simple shears are prescribed to the molecular box to identify the atomistic elastic tensor This post-processing has shown that the stiffness tensor components with Voigt notation verify that ∀x, 4a,b show the evolution of (C w ) 11 , (C w ) 12 , (C w ) 13 , (C w ) 33 and (C w ) 66 with the distance to the graphene z. Substantial increase of one order of magnitude of the stiffness components in graphene plane ((C w ) 11 , (C w ) 12 and (C w ) 66 ) is shown for z < 2.5Å due to graphene sheet (see Figure 4a).

(C w ) 11 (x) ≈ (C w ) 22 (x), (C w ) 31 (x) ≈ (C w ) 32 (x) and ∀I ≥ 4, (C w ) I1 (x) ≈ (C w ) I2 (x) ≈ 0; and ∀I < 6, (C w ) I6 (x) ≈ 0. The Figure
In the interphase area, 2.5Å ≤ z ≤ 10.5Å, the stiffness components (C w ) 11 , (C w ) 12 and (C w ) 33 are 40% higher than in the bulk part z ≥ 10.5Å (see Figure 4b), while (C w ) 13 and (C w ) 66 seem to be the same as the polymer bulk. Thus, the thickness of the interphase region is identified as 10.5 Å.

Due to the slip at the polymer/graphene interface (see Figure 5b), the local strain, w (x) and the local stress increment, σ w (x), in the polymer remains zero for all shears in the graphene plane. It is therefore impossible to identify with this method the stiffness tensor components (C w ) I4 and (C w ) I5 . For convenient reasons, we assume that (C w ) 44 (x) = (C w ) 55 (x) = (C w ) 66 (x) and C I4 (x) = C I5 (x) = 0 for all x in the polymer region.

The stiffness increase in the interphase region appears to be related to the increase in atomic density, ρ w (x) which results in an increase in the van der Waals interaction density. Figure 4d shows the local density of graphene, polymer and composite systems as a function of the z distance from the graphene sheet. The density of the graphene polymer composite with sandwich structure is dominated by the density of graphene in the midplane and by the density of polymer for the rest part. We can see a dense packing of polymer mass close to graphene layer called interphase zone, which is due to the adsorption of the chains on faces of graphene sheets. In addition to the density effect, the stiffness increase in the interphase along the graphene sheet plan ((C w ) 11 , (C w ) 22 and (C w ) 12 ) is also due to the orientation of the covalent bonds inside of polymer macromolecules into the plane (e 1 , e 2 ). This local conformation change is evidenced by the nematic order parameter s w (see Figure 4c), which is defined by

s w (x) = 1 ρ b (x) B (p) 3 (n αβ .e 3 ) 2 -1 2 b w (x; r α (t) , r β (t)) (38)
where B (p) is the set of covalent bonds in the polymer; e 3 is the unit normal of graphene; n αβ = (r αr β )/r αβ is the unit vector which defined the covalent bond direction between two CH 2 atom group α and β along the polymer chain; b w (x; r α (t) , r β (t)) is the bond function defined Eq. ( 27) and

ρ b (x) = B (p) b w (x; r α (t) , r β (t)
) is a normalization constant which corresponds to the local bond density at x. We can see on the Figure 4c a decrease of the nematic order parameter s w in the interphase. Generally, the nematic order s w is used to describe the orientational order of a nematic liquid crystal. For a completely random and isotropic sample, s w = 0, whereas for a perfectly aligned sample along e 3 axis, s = 1. In addition, s = -1/2 denotes that all the bonds are in the plan perpendicular to e 3 . Therefore, it indicates that in the interfacial zone close to the graphene sheets, the polymer chains turn to be paralleled to the graphene layer. Moreover, the ALIAS method makes it possible to define a displacement field u w . Figure 5 represents the components of displacement fields, u w along the normal of graphene sheet for an elongation of n = 0.3% normal to the graphene sheet, ¯ = n e 3 ⊗ e 3 (see Figure 5a) and for an simple shear of γ = 3% in a direction e 1 on a graphene plane, ¯ = γ 2 (e 3 ⊗ e 1 + e 1 ⊗ e 3 ) (see Figure 5b).

As shown in Figure 5a, the component u w • e 3 of displacement field evolves linear in polymer bulk (z < 10.5 and z > 10.5), which means that the strain is homogeneous there (≈ 0.32%). In the central part corresponding to the interphase, discontinuities are observed. Note that the average slope is lower there (≈ 0.23%), which is in agreement with the fact that the interphase is stiffer.

Figure 5b shows sharp discontinuities of the displacement fields component u w • e 1 on both sides of graphene sheet. Moreover, the displacement field remains constant in the polymer region. These two facts indicate that the deformation applied to the atomic box is accommodated by slips at the interfaces between graphene and polymer. For the identification of interface parameters, we write the surface free energy φ s as the sum of three contributions as shown in Figure 6: one for the graphene and two for the interfaces polymer/graphene.

Identification of imperfect interface elastic parameters

Therefore, the three layers description of the imperfect interface is used. We assume that the elastic parameters of the imperfect interface are constant along the graphene sheet for the sake of simplicity and consistency with the previous section.

Following the classification in [START_REF] Chatzigeorgiou | Generalized interfacial energy and size effects in composites[END_REF], the graphene layer can be considered as an elastic interface, denoted by s2, i.e. it is kinematically coherent ( u | s2 = 0) but kinetically non-coherent ( t | s2 0). The free energy of s2 is defined by

φ s2 = 1 2 s2 : C s2 : s2 + τ s2 : s2 + φ s2 0 , ( 39
)
where C s2 is the fourth order surface stiffness tensor of graphene, τ s2 is the internal surface stress and s2

is the surface strain tensor into the graphane sheet.

The graphene/polymer interfacial regions can be assumed as cohesive interfaces, denoted s1 and s3, i.e.

it is kinematically non-coherent ( u | s1 0 and u

| s3 0) but kinetically coherent ( t | s1 = t | s3 = 0).
The cohesive laws for interface s1 and s3 is

σ| s1 .e 3 = K s1 u | s1 , σ| s3 .e 3 = K s3 u | s3 , ( 40
)
where K s1 and K s3 are the second order cohesive stiffness matrix of the interface s1 and s3. The free energy of s1 and s3 are defined by

φ s1 = 1 2 u | s1 .K s1 . u | s1 + φ s1 0 , φ s3 = 1 2 u | s3 .K s3 . u | s3 + φ s3 0 , ( 41
)
where

φ s 0 = φ s1 0 + φ s2 0 + φ s3 0 .
The kinematic compatibility across the graphene sheet s2 enforces that the jump of displacement field in the equivalent imperfect interface is

u = u| + s3 -u| - s1 = u| + s3 -u g + u g -u| - s1 = u | s1 + u | s3 . ( 42
)
In equation [START_REF] Gu | Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfaces[END_REF], "+" and "-" correspond to the upper and lower surface of the interface with respect to the normal direction e 3 . The results of atomistic simulations suggest that the displacement jumps across interface s1 and s3 are equal. For the sake of simplicity, we assume that u

| s1 = u | s3 = u /2, which implies that the displacement field, u g (u g = u| + s2 = u| - s2 = u| + s1 = u| - s3
), and the strain fields, s2 , on the graphene sheet are given by

u g = u| + s3 + u| - s1 2 , s2 (x) = 1 2 ∇ s u g (x) + (∇ s u g ) T (x) . ( 43
)
Keeping in mind the cohesive law for the equivalent imperfect interface as

K s u = t s = ( σ| s1 + σ| s3 ) .e 3 2 , ( 44
)
introducing the cohesive laws of the interface s1 and s3, and using the equality of jump displacement across s1 and s3, we obtain

K s = 1 4 K s1 + K s3 = 1 2 K s1 = 1 2 K s3 . ( 45
)
The two last equalities come from the mirror symmetry of the system with respect to the graphene plane which implies that

K s1 = K s3 .
The additivity of specific free energy and the equation [START_REF] Chatzigeorgiou | Generalized interfacial energy and size effects in composites[END_REF] gives

φ s = φ s1 + φ s2 + φ s3 , ( 46
)
φ s = u | s1 .K s . u | s1 + 1 2 s : C s2 : s + u | s3 K s . u | s3 + φ s 0 . ( 47
)
After establishing the new expression of specific free energy of the imperfect interface φ s , the identification of stiffness tensor K s is detailed. At atomic scale we impose elementary rigid body translation of graphene sheet, d = d 1 e 1 + d 2 e 2 + d 3 e 3 (see e.g. Fig. 7a for translation along e 1 ). The polymer atoms displacements are kept fixed so that the variation of potential energy is only due to the displacement discontinuity at the interface and u | s1 = -u | s3 = d. Assuming that the temperature is 0 K in our simulations, the variation of free energy for a translation d of graphene is given by the variation of potential energy

∆U = φ -φ (0) = Ω φ -φ 0 dΩ + I φ s -φ s 0 dS (48) = S ( u | s1 .K s . u | s1 + u | s3 K s . u | s3 ) (49) = 2S (d • K s • d) ( 50
)
where S is the graphene sheet surface. The last equation shows that the potential energy is a quadratic function of d defined by the six independant components of the symetric second order stiffness tensor K s .

For the identification, we compute a set of 37 values of the variation of potential energy ∆U to estimate the stiffness tensor K s by the least squares method (see supplementary information for details). We assume that the stiffness tensor K s is diagonal (K s ) 23 = (K s ) 31 = (K s ) 12 = 0 MPa.nm -1 , for the sake of simplicity.

The cohesive part of the general imperfect interface is very anistropic because the stiffness associated to the opening mode perpendicular to the graphene (K s ) 33 = 52960 MPa.nm -1 is five order of magnitude than the stiffness associated to the slips along graphene plane (K s ) 11 = (K s ) 22 = 8.50 MPa.nm -1 .

This identification method allows to identify the Critical Resolved Shear Stress (CRSS) associated to the slip in the graphene plane by computation of the local shear stress. Fig. 7b shows the evolution of the The surface elasticity tensor is identified by the method used for the bulk. Indeed, the fourth-order stiffness tensor C w is identified on the graphene sheet (z = 0). Then, C s , is deduced from (see e.g. [START_REF] Gu | Coordinate-free derivation and weak formulation of a general imperfect interface model for thermal conduction in composites[END_REF]).

C s = h C w - (C w : N) ⊗ (C w : N) N : C w : N , ( 51
)
where N = e 3 ⊗e 3 and h = 4Å is the graphene sheet thickness. The h value is chosen as the full width at half maximum of the curve of graphene density ρ w vs. z (see Figure 4d). In this work, C s is given

C s 11 = C s 22 = 280
GPa.nm, C s 12 = 90 GPa.nm, C s 66 = 95 GPa.nm and C s ij = 0 GPa.nm for other components.

Discussion

Summary of the identification of elastic parameters

The results of the previous section suggest that polymer/graphene nanocomposite can be modeled by a five layer composite: an imperfect interface for graphene (I), two interphase layers (Ω 2 and Ω 3 ) with a thickness, t I = 10.5 Å and two layers for polymer bulk (Ω 1 and Ω 4 ) see Figure 8b. The mean values of the elastic constants in Voigt notation for each layer are given in Table 2.

The elastic behavior of the polymer bulk is isotropic and can be defined by the Lame's coefficients ), with a huge mechanical contrast between the out of graphene plane stiffness and the in plane stiffness (K s 33 /K s 11 ≈ 10 5 ).

It is important to note that the results depend on the weight function, w (x) and therefore on the shape and size of our measuring tools. Nevertheless, the mean strain computed with the atomistic displacement, u w (x) is equal to the the effective strain of the atomistic simulation box ¯ , and that the effective properties are the same for the atomic and continuum models (cf. following section and Table 2). Another area for improvement is taking the temperature into account. Indeed, all the results are obtained at zero temperature which explains why the Young's modulus of polymer bulk, E b ≈ 6.4 GPa, has the same order of magnitude as that of an amorphous polymer below its secondary transition temperature [START_REF] Halary | Polymer materials: macroscopic properties and molecular interpretations, chapter Chapter 4 Secondary relaxations in amorphous polymers[END_REF]. We estimate that the ratio between Young's modulus of interphase, E I and bulk remains qualitatively the same E I /E b ≈ 1.5 for temperatures below the glass transition temperature of the bulk and for much higher temperatures. We issue a large reserve for a temperature range around the glass transition temperature of the bulk because it is possible that the glass transition temperature of the interphase is different [START_REF] Fryer | Dependence of the glass transition temperature of polymer films on interfacial energy and thickness[END_REF][START_REF] Bansal | Quantitative equivalence between polymer nanocomposites and thin polymer films[END_REF], which would generate a change of one or two orders of magnitude of the ratio E I /E b .

Ω 1 Ω 2 Ω 3 Ω 4 e 1

Effective properties of sandwich structure

We use the framework derived by Chatzigeorgiou et al [START_REF] Chatzigeorgiou | Generalized interfacial energy and size effects in composites[END_REF] to define the macro strain tensor (or the mean strain tensor of the sandwich RVE of the Fig. 8 a):

= 1 |Ω| Ω dV + 1 |Ω| I 1 2 (n I ⊗ u + u ⊗ n I ) dS, ( 52
)
where n I is the unit normal vector to the imperfect interface at a given point. The macro stress σ (or the mean stress tensor of the sandwich RVE) is defined by

σ = 1 |Ω| Ω σ dV + 1 |Ω| I σ s dS ( 53
)
The effective stiffness tensor C C of the equivalent continuum model, which linearly connects the macro stress and strain (σ = C C : ), is computed by numerical homogenization with periodic boundary condi-

Polymer bulk Polymer interphase

C b =       
13.6 9.0 9.0 0 0 0 9.0 13.6 9.0 0 0 0 9.0 9.0 13.6 tions satisfying the extended Hill-Mandel condition [START_REF] Chatzigeorgiou | Generalized interfacial energy and size effects in composites[END_REF] (see Supplementary Information for details). We use the finite element method to discretize the solution space with linear tetrahedrons for bulk part and linear triangles for interfaces. The effective behavior of the sandwich structure, presented in Table 2, is transversally isotropic because the microstructure and the elastic properties of phases are isotropic along the graphene plane. In addition, the effective stiffness depends on the RVE size L z . It is possible to obtain an analytic expression of the effective stiffness coefficient (C C ) ij in the case of RVE with a sandwich structure.

0 0 0 0 0 0 2.3 0 0 0 0 0 0 2.3 0 0 0 0 0 0 2.3        GPa C I =        17
The theoretical study of this size effect will be the subject of a future work, but here we preferred to study the influence of the interface and the interphase on the RVE anisotropy.

To validate the ALIAS methodolgy, we estimate directly by atomic simulations the effective stiffness tensor, C A , following the procedure introduced for polymer based materials by Theodorou and Suter [START_REF] Theodorou | Atomistic modeling of mechanical properties of polymeric glasses[END_REF] (see Supplementary Information for details). We believe that our approach is validated because we find the same effective module To show this anisotropy, we plot in Figure 9a the Young's modulus E(θ) normalized by Young's modulus of polymer E b as function of the θ angle (see Figure 8) between the tensile direction and the vector e 1 into the plane (e 1 , e 3 ). We use the effective stiffness tensor C C to compute E(θ), the demonstration of the used formula is provided in Supplementary Informaiton. We notice that for tensile directions that are perpendicular to the graphene (θ = 90 • ) and aligned with graphene (θ = 0 • ) the Young's modulus is greater than that of the polymeric matrix. It is even five times larger in the graphene plane due to the extraordinary stiffness of graphene. On the other hand, for a wide angular range (6 • ≤ θ ≤ 86 • ), the Young's modulus is much lower than that of the polymer, going up to more than one order of magnitude for angular range

C C = C A ,
(22 • ≤ θ ≤ 68 •
). This is due to the superlubricity of graphene (see e.g. [START_REF] Chen | Superlubricity of carbon nanostructures[END_REF], for review) which generates very low shear stiffness and CRSS.

To compare the influence of the interphase and the interface on the elastic properties, we plot in Figure 9a the Young's modulus normalized as a function of the θ angle for one model "without interphase" (C b = C i ) and another one "without interface", i.e. without displacement discontinuities at graphene interfaces

( u = 0, K s ii → ∞).
It is clear that the interphase has no significant effects on the Young modulus because the curves of the models with and without interphase are superposed; with the exception of an increase in question question the Young's modulus of 10% for a tensile direction perpendicular to graphene (θ = 90 • ), which is not visible on the graph. Conversely, the interface acts as a notable softener effect because the model without interface has a Young's modulus always higher than that of the polymer.

L z = 10 nm 2t i = 2 nm L x = 20 nm 18 nm L x = 4 nm 2t i = 2 nm ∆ = 0.2 nm
We can conclude that one of the ways to increase the Young's modulus of nanocomposite without microstructure change is to increase the interfacial cohesion between the polymer and the graphene by chemical route such as grafting. Nevertheless, this involves a compromise for the material designer between stiffness and electrical conductivity. In fact, the exceptional electrical conductivity of graphene is altered when the perfection of the 2D crystal is broken. The chemical interface modification will also modify the CRSS, which would tend to increase the yield stress and decrease the toughness. Indeed, the toughness material is able to dissipate energy by anelastic mechanisms such as interfacial sliding before propagating a critical crack. Therefore, the CRSS increasing, induced by the chemical grafting of graphene, should decrease the number of active slip sites in a volume element and thus reduce the toughness. We fall back into one of the classic conflicts between strength and toughness [START_REF] Ritchie | The conflicts between strength and toughness[END_REF].

Another way to increase the Young's modulus is to change the geometry of the system by either adding the GNP extremities or introducing a graphene wrinkle. Figure 10 shows model microstructures, one with the GNP extremities and the other with wrinkled graphene. Figure 9b shows the effect of these geometric changes on the Young's modulus anisotropy. Taking into account the GNP extremities increases the Young's modulus by one order of magnitude for angular range 20 • ≤ θ ≤ 70 • . However the same simulation with "without interface", i.e. without displacement discontinuity ( u = 0), shows an increase in Young's modulus going up to double for θ = 45 • . The other geometric effect which increases the Young's modulus is the GNP wrinkling. In this case, the effect on Young's modulus is similar to that of taking into account the GNP extremity (see Figure 9b) with higher moduli in directions (θ = 0 • and 90 • ). Note that the curvature effects is taken into account in the imperfect interface model which is a generalization of the Young-Laplace equation [see, e.g. 105]. In this simulation, we negled also the bending stiffness of graphene which could be considered extremely low [START_REF] Wei | Bending rigidity and gaussian bending stiffness of single-layered graphene[END_REF][START_REF] Sajadi | Size-and temperature-dependent bending rigidity of graphene using modal analysis[END_REF]. The wrinkling of graphene also induces a decrease in the GNP conductivity due to the presence of midgap state [START_REF] Xu | Scanning tunneling microscopy characterization of the electrical properties of wrinkles in exfoliated graphene monolayers[END_REF]. Here, we also find the conflicts between stiffness and conductivity. The influence of these geometric changes on the yield stress and toughness of the GNP polymer nanocomposite is not obvious. In our opinion, it would merit special study because of deformation mechanism such as cavitation [START_REF] Li | Atomistic simulations on multilayer graphene reinforced epoxy composites[END_REF] or crazing [START_REF] Lu | Multiscale study of influence of interfacial decohesion on piezoresistivity of graphene/polymer nanocomposites[END_REF]. The Table 3 summarizes all the conclusions resulting from this discussion, giving the supposed influence of changes in interface properties and geometry on the effective properties of GNP polymer nanocomposites.

Conclusion

In this paper, we have tackled the problem of identification of elastic continuum model by atomistic simulations for graphene polymer nanocomposite. The ALIAS method was developed to estimate the local stiffness tensor on all points of polymer graphene nanocomposite with a sandwich structure model. Results

suggest that the graphene can be modeled at continuum scale by a general imperfect interface with zero thickness. The cohesive part of the general imperfect interface is very anistropic because the stiffness associated to the slip along graphene plane is five orders of magnitude lower than the stiffness associated to the opening mode perpendicular to the graphene. Moreover, the identification procedure revealed the existence of interphase on either side of the graphene with a thickness of 1 nm, which is one and a half times stiffer than the polymer bulk matrix.

The identified continuum model has been used to study the effective elastic properties of nanocomposites with sandwich microstructure by numerical homogenization. This study at continuum scale showed a softening effect due to the very low stiffness of slips along graphene plane is preponderant in relation to the interphase stiffening. Finally, the continuum model suggests that the wrinkling of graphene increases the stiffness of nanocomposites, as well as the increasing of interfacial cohesion by chemical route such as grafting of graphene.

Aknowledgements

X. Lu thanks the financial support of China Scholarship Council(CSC) for the Ph.D. thesis. The molecular simulations were carried out using the LAMMPS molecular dynamics software [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF]. cas of coarse-grain PE and graphene, these potentials are written as

θ αβγ r α r β r γ φ αβγδ r α r β r γ r δ a) b)
U (s) ({r α }) = Bond k (s) αβ 2 r αβ -r 0 αβ 2 , ( 2
)
U (b) ({r α }) = Bending k (b) αβγ 2 cos (θ αβγ ) -cos θ 0 αβγ 2 , (3)
U (tor) ({r α }) = Dihedral k (tor) βγ 2 1 -3 cos (φ αβγδ ) + 4 cos 3 (φ αβγδ ) , (4)
where k are parameters which depends on the atom types. The bending angle θ αβγ and the dihedral angle φ αβγδ are difened on Figure 1 and they can be computed using interatomic distances(see, e.g. [2] p. 468 for expression of dihedral angle cosine):

cos (θ αβγ ) = r βα .r βγ r αβ r βγ = r 2 αβ + r 2 βγ -r 2 αγ 2r αβ r βγ , (5) cos
(φ αβγδ ) = (r αβ × r βγ ) • (r βγ × r γδ ) |r αβ × r βγ | |r βγ × r γδ | (6) = r 2 αβ + r 2 βγ -r 2 αγ r 2 βγ + r 2 γδ -r 2 βδ -2r 2 βγ r 2 βγ -r 2 αγ -r 2 βδ + r 2 αδ 4 r 2 αβ r 2 βγ -1 4 r 2 αβ + r 2 βγ -r 2 αγ 2 r 2 γδ r 2 βγ -1 4 r 2 γδ + r 2 βγ -r 2 βδ 2 . ( 7
)
The total potential energy of the system is represented by the sum of all these contributions:

U ({r α } , ¯ ) = U (vdw) ({r α }) + U (s) ({r α }) + U (b) ({r α }) + U (tor) ({r α }) . ( 8
)

Fitting details of the cohesive stiffness of interface K s

The variation of potential energy ∆U of the atomistic box for a translation d of graphene is quadratic form given by

∆U = 2S (d -d 0 ) • K s • (d -d 0 ) ( 9
)
where S is the graphene sheet surface and d 0 a small shift vector to correct the numerical error due to the numerical minimization of the initial configuration of atomistic box. The cohesive stiffness of interface K s and the shift vector d 0 matrix tensor K s is defined by the least squares method, for a given set of K = 37 translations d (k) , such as

{K s , d 0 } = argmin K s ∈M diag 3 ,d0∈R 3 K k=1 ∆U (k) -2S d (k) -d 0 • K s • d (k) -d 0 2
(10) where M diag 3 is the set of diagonal matrix in basis (e 1 , e 2 , e 3 ) and ∆U (k) is the variation of potential energy for the translation d (k) . The 37 values of translation vector are chosen along six directions of translation {e 1 , e 2 , e 3 , e 2 + e 3 , e 3 + e 1 , e 1 + e 2 }. The Figure 2 shows the atomistic value of the variation of potential energy ∆U and the quadratic fit for a translation into graphene plane along e 1 and another one perpendicular to graphene plane e 3 ). The cohesive part of the general imperfect interface is very anistropic because the stiffness associated to the opening mode perpendicular to the graphene K s 33 = 52960 MPa.nm -1 is five order of magnitude than the stiffness associated to the slips along graphene plane K s 11 = K s 22 = 8.50 MPa.nm -1 .

Computation of effective stiffness tensor by atomistic simulations

The method presented in the previous section is a transposition of the procedure introduced by Theodorou and Suter [3] to estimate the stiffness tensor of polymer based materials. The starting point of this procedure is to write the variation of the total free energy ∆φ as function of the effective infinitesimal strain tensor ¯ , such as

∆φ = |Ω| 2 ¯ : C A : ¯ (11)
Assuming that the temperature is 0 K in our simulations, the variation of free energy for a prescibed effectice strain ¯ is given by the variation of potential energy

∆φ = ∆U = U (¯ ) -U (0) (12)
The effective stiffness C A is defined by the least squares method, for a given set of K = 43 effective strain ¯ (k) , such as

C A = argmin C A ∈Orth K k=1 ∆U (k) - |Ω| 2 ¯ (k) : C A : ¯ (k) 2 ( 13
)
where Orth is the set of fourth order elastic tensor for orthorhombic material in basis (e 1 , e 2 , e 3 ) and ∆U (k) is the variation of potential energy for the prescribe effective strain 3 ¯ (k) . Note that this method is not limited to orthorhombic material, this simplifying assumption that guided by laminate geometry of atomistic box allowed us to reduce the coefficient number of 2 (e 3 ⊗ e 3 + e 1 ⊗ e 1 ) . For each direction α (j) , we prescribe six effective strains, such as

¯ (k) = λ i α (j) , k = i + 7 (j -1) (14)
with λ i = {-0.015; -0.01; -0.005; 0.005; 0.01; 0.015}

The effective stiffness tensor computes by atomistic simulations is

C A =         43

Computation of effective stiffness tensor by numerical homogenization

In this section, the definition and computation of effective properties in the context of linear elasticity with imperfect interface are presented. The practical calculation of the effective elastic tensor with 3D Finite Element Method (FEM) is detailed. An overview of computational homogenization with FEM can be found in [4]. The Representative Volume Element (RVE) for periodic laminate polymer/graphene nanocomposte is composed by a five layers: an imperfect interface for graphene (I), two interphase layers (Ω 2 and Ω 3 ) with a thickness, t I = 10.5 Å and two layers for polymer bulk (Ω 1 and Ω 4 ) see Figure 3.

Ω 1 Ω 2 Ω 3 Ω 4
The continuum homogenization problem being linear, the effective stiffness of polymer/graphene nanocomposite C C given the macroscopic constitutive relationship betwen macro stress σ and macro strain σ = C C : [START_REF] Brown | A molecular dynamics study of a model nanoparticle embedded in a polymer matrix[END_REF] is obtained by computing the macro stress for six elementary localization problems on RVE with the following prescribed macro strains :

(1) = e 1 ⊗ e 1 ;

(4) = 2 (e 2 ⊗ e 3 + e 3 ⊗ e 2 ) ;

(2) = e 2 ⊗ e 2 ;

(5) = 2 (e 3 ⊗ e 1 + e 1 ⊗ e 3 ) ;

(3) = e 3 ⊗ e 3 ;

(6) = 2 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) .

Strong form

The localization problem based assuming that the RVE is subjected to the prediodic boundary condition is given as follows:

Given a macroscopic strain , find the displacement field u (x) in Ω such that:

∇ • σ = 0 ∀x ∈ Ω, (18)
∇ s • σ s + t = 0 ∀x ∈ I, (19)
σ = C (x) : (x) ∀x ∈ Ω, (20)
σ s = C s (x) : s (x) ∀x ∈ I, (21)
t s = 1 2 σ| + I + σ| - I • n I = K s • u ∀x ∈ I, (22)
with the kinematic

(x) = 1 2 ∇u (x) + ∇u T (x) , ∀x ∈ Ω, (23) s
(x) = 1 2 ∇ s u s (x) + (∇ s u s (x)) T ∀x ∈ I, (24)
u s (x) = 1 2 u| + I + u| - I ∀x ∈ I, (25)
u = u| + I -u| - I ∀x ∈ I. (26)
and verifying the periodic boundary conditions

u (x) = • x + ũ (x) , ∀x ∈ ∂Ω ( 27
)
where the fluctuation ũ (x) is periodic on Ω and discoutinous over the imperfect interface I. In other words, it takes the same values at two homologous points on opposite faces of a parallelepipedic domain Ω. The traction vectors σ • n and σ s • n are antiperiodic.

Weak form

Given a macroscopic strain , find ũ ∈ H p 1 (Ω) for all δu ∈ H p 1 (Ω), such as :

Ω σ [u] : [δu] dV + I σ s [u] : s [δu] + t s [ u ] • δu dS = 0 (28)
where u = • x + ũ and H p 1 (Ω) is the set {w (x) periodic on Ω and continous on Ω\I, such as Ω w 2 dV < ∞ , Ω ∇w : ∇w dV < ∞ }.

The equation (28) becomes after using the elastic constitutive laws

Ω [ũ] : C (x) : [δu] dV + I s [ũ] : C s (x) : s [δu] + ũ • K s • δu dS (29) = - Ω : C (x) : [δu] dV - I s : C s (x) : s [δu] dS
where s = P • • P is the projection of the macro strain onto the tangent plane of the interface I and P = In I ⊗ n I .

Finite element discretization

For the discretization of this discontinous weak form it suffices to divide the domain into two subdomains, Ω -and Ω + , separated by the interface I and write the weak as follow :

Given a macroscopic strain , find ũ ∈ H p 1 (Ω) for all δu ∈ H p 1 (Ω), such as :

m ++ (ũ, δu) + m --(ũ, δu) + m +-(ũ, δu) + m -+ (ũ, δu) = f + (δu) + f -(δu) (30)
where

m ++ (ũ, δu) = Ω + [ũ] : C (x) : [δu] dV + I + 1 4 s ũ| + I : C s (x) : s δu| + I dS + I + ũ| + I • K s • δu| + I dS, m --(ũ, δu) = Ω - [ũ] : C (x) : [δu] dV + I - 1
• K s • δu| + I dS, f + (δu) = - Ω + : C (x) : [δu] dV - I + 1 2 s : C s (x) : s δu| + I dS, f -(δu) = - Ω - : C (x) : [δu] dV - I - 1 2 s : C s (x) : s δu| - I dS.
The bilinear form m ++ (ũ, δu) and the linear form f + (δu) (reps. m --(ũ, δu) and f -(δu)) are discretized by Lagrangian finite-element method with linear tetrahedron P1-elements of subdomain Ω + (resp. Ω -) and computed by Gauss quadrature of order 1 and order 2 for the term I + ũ|

+ I • K s • δu| + I dS (resp. I -ũ| - I • K s • δu| - I dS).
To compute the coupling term m -+ (ũ, δu) and m -+ (ũ, δu), we impose the same discretization of interface I for the domains Ω + ans Ω -.

The using 3D meshes are plotted on number of degree of freedom (dof) is doubled. Note that the peridodic boundary conditions and the continuous translational symmetry of all RVE geometry along e 2 (along e 1 and e 2 for flat graphene RVE) imply that the effective stiffness C C is independent of the lenght L y for all studied RVE (and independent of L x for flat graphene RVE). We chose L y = 4Å with 4 rows of elements following the direction e 2 to reduce the computational cost. The choice of L x is arbitrary for all RVE while the L z = 10 nm is chosen to have the same size as the atomic simulation box. The shape of wrinkled graphene is sinusoidal given by the parametric equations :

x = u u ∈ - L x 2 ; L x 2 y = v v ∈ - L y 2 ; L y 2 z = ∆ 2 cos 2π u L x
where ∆ = 2Å and L x = 4 nm. This choice gives the same order of magnitude graphene curvature that this one observed by dynamics molecular simulations at 300K [5].

Case of finite size graphene. In this case, we add a Neumman boundary condition at the interface border ∂I, such as

σ s • m = 0 ( 31
)
where m is the normal vector to the line ∂I and perpendicular to the normal unit vector of the interface n I . Note that adding the boundary condition doesn't change the expresion of weak form (Eq. 28).

Case of wrinkled graphene.

where C s ij and K s ii are the corresponding components for the flat graphene case; and where n I [x] is unit normal to the imperfect interface at point x and (t 1 [x] , t 2 [x]) is the orthonormal basis of tangent plane to the imperfect interface at point x. The vector t 2 [x] = e 2 ∀x ∈ I because of the RVE is invariant by translation along e 2 .

Post-processing

The macro stress σ defined by is evaluated by Gauss quadrature of order 1.

The prescribed macro strain are chosen the Voigt's tensors basis so the effective stiffness tensor is given in Voigt notation by : C C = σ (1) ; σ (2) ; σ (3) ; σ (4) ; σ (5) ; σ (6) (35

)
where σ (i) is the macro stress in Voigt notation obtain for prescribed the macro strain (i) . The effective stiffness and effective compliance for studied cases are given in Table 1. For macroscopic tensile loading along direction e t , the macro stress tensor is :

σ = σ t e t ⊗ e t , ( 36
)
and the the macro strain tensor is given by = S : σ [START_REF] Marcadon | Confrontation between molecular dynamics and micromechanical approaches to investigate particle size effects on the mechanical behaviour of polymer nanocomposites[END_REF] where the fourth order tensor S = C -1

CM is the effective compliance tensor. We define the effective Young's modulus associated to the tensile loading along direction e t as E (e t ) = σ t t [START_REF] Le | Stochastic continuum modeling of random interphases from atomistic simulations. Application to a polymer nanocomposite[END_REF] where t is the normal strain associated to the tensile loading along direction e t defined by

t = : e t ⊗ e t . ( 39
)
In case of tensile loading along direction e t = cos (θ) e 1 + sin (θ) e 3 , the macro stress tensor is σ = σ t cos 2 (θ) e 1 ⊗ e 1 + sin 2 (θ) e 3 ⊗ e 3 + cos (θ) sin (θ) (e 3 ⊗ e 1 + e 1 ⊗ e 3 ) . [START_REF] Bövik | On the modelling of thin interface layers in elastic and acoustic scattering problems[END_REF] Therefore, the macro stress and macro strain are in Voigt notation associated to the basis (e 1 , e 2 , e 3 ).

where S ij are effective compliance tensor components for transversally isotropic material. The macro strain can be rewritten as Finally, the effective Young's modulus E(θ) as function of the angle, θ, between the tensile direction and the vector e 1 into the plane (e 1 , e 3 ) is

=
E (θ) = 1
S

  are based on the following conventions. Scalars are in italics lowercase (a), vectors are boldface in lowercase (a), second-order tensors are bold-face in uppercase or bold-face in greek symbol (A, α), and fourth-order tensors are blackboard bold in uppercase (A). The dyadic product of two vectors a and b is a second-order tensor D = a ⊗ b with (D) ij = (a) i (b) j . The scalar product of two vectors a and b is denoted a • b = (a) i (b) i . The scalar product of two second-order tensors A and B is denoted

Figure 1 :

 1 Figure 1: Sketch of domains Ω with the imperfect interface I and the unit normal n I .

+ I and u| -I

 u| are the values of displacement in the upper and lower surfaces of interface I respectively.

Figure 2 :

 2 Figure 2: Graphene polymer nanocomposite with sandwich structure under periodical boundary condition, single layer graphene is located in the middle of the box.At the atomic scale, the graphene polymer nanocomposite is modeled by a set of N classical interacticting point particules with sandwich structure in a rectangular box under periodic boundary condition (see Figure2). Two types of particles are considered: the carbon in the graphene sheet and the -CH 2 -atom group

Figure 3 :

 3 Figure 3: Profiles for weighting function a) along e 1 and b) along e 3 .

CFigure 4 :

 4 Figure 4: (a) and (b) Evolution of stiffness tensor components with Voigt notation as a function of the normal distance, z, from the graphene layer; (c) Nematic order parameter, sw vs z; (d) Evolution of density ρw of polymer, graphenes and nano composite as a function of the normal distance, z.

Figure 5 :

 5 Figure 5: Evolution of displacement components as a function of the normal distance, z, from the graphene layer. (a) uw • e 3 for a elongation of n = 0.3% normal to the graphene sheet, ¯ = ne3 ⊗ e 3 . (b) uw • e 1 for a simple shear of γ = 3% in a direction e 1 on a graphene plane, ¯ = γ 2 (e 3 ⊗ e 1 + e 1 ⊗ e 3 ).

IFigure 6 :

 6 Figure 6: Scheme of the three layers of the imperfect interface I.

Figure 7

 7 Figure 7: a) Scheme of graphene sheet translation along e 1 . b) Evolution of shear stress in a graphene plane for a translation along e 1 .

λ b = 9 . 0

 90 GPa and µ b = 2.3 GPa, such as C b 11 = C b 22 = C b 33 = λ b + 2µ b ; C b 12 = C b 23 = C b 31 = λ b and C b 44 = C b 55 = C b 66 = µ b . The elastic behavior of the polymer interphase is transverse isotropic, it can be defined by five coefficients, such as C I 11 = C I 22 = 17 GPa; C I 12 = 11.4 GPa, C b 23 = C b 31 = 9.1 GPa; C I 44 = C I 55 = 2.3 GPa and C I 66 = (C I 11 -C I 12 )/2. The surface elasticity of imperfect interface is isotropic, it is defined by surfacic Lame's coefficients λ b = 95 GPa.nm and µ b = 90 GPa.nm, such as C s 11 = C s 22 = λ s + 2µ s ; C b 12 = λ s and C b 66 = µ s . The cohesive stiffness of imperfect interface is isotropic into the graphene plane (K s 11 = K s 22 and K s 12 = K s 23 = K s 31 = 0 MPa.nm -1

Figure 8 :

 8 Figure 8: a) Atomistic configuration; b) Equivalent continuum model with imperfect interface and interphases.

  with the exception of the components (C A ) 44 and (C A ) 55 which cannot be obtained by atomistic simulations because of the sliding at the polymer/graphene interface (cf. the end of section 3

Figure 9 :

 9 Figure 9: Young's modulus normalized by that of the bulk polymer E b as function of the angle, θ, between the tensile direction and the vector e 1 into the plane (e 1 , e 3 ) (see Supplementary Information for compuation details of E(θ)). a) Comparison between the full continuum model (with interface and interphase) and a model without interphase and another one without the cohesive part of imperfect interface ( u = 0) b) Illustration of the effect finite size of graphene and wrinkle of graphene on the stiffness anisotropy.

Figure 10 :

 10 Figure 10: RVE scheme for finite element simulations: a) with finite size graphene; b) with the wrinkled graphene. Note that the 3D RVE model is obtained by an extrusion of these 2D sketches. The associated meshes are given in supplementary information.

Figure 1 :

 1 Figure 1: Definition of: a) the bending angle θ αβγ and b) the dihedral angle φ αβγδ .

Figure 2 :

 2 Figure 2: a) Scheme of graphene sheet translation along e 1 . b) Evolution of potential energy as function of d for translation of graphene along e 1 . c) Scheme of graphene sheet translation along e 3 . d) Evolution of potential energy as function of d 3 for translation of graphene along e 3 .

1 √ 2 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) , 1 √ 2 (e 1 ⊗ e 1 + e 2 ⊗ e 2 ), 1 √ 2 (e 2 ⊗ e 2 + e 3 ⊗

 121221211212223 C A from 21 to 9 ((C A ) 11 , (C A ) 22 , (C A ) 33 , (C A ) 44 , (C A ) 55 , (C A ) 66 , (C A ) 12 , (C A ) 23 , (C A ) 31 ). Moreover, the elastic constant (C A ) 44 and (C A ) 55 cannot be obtained by atomistic simulations because of the sliding at the polymer/graphene interface, so we impose (C A ) 44 = (C A ) 55 = 0. The 43 values of effective strain ¯ (k) are chossen along seven directions α (j) in space of symetric second order tensor {e 1 ⊗ e 1 , e 2 ⊗e 2 , e 3 ⊗e 3 , e 3 ), 1 √

Figure 3 :

 3 Figure 3: Representative Volume Element for periodic laminate polymer/graphene nanocomposte (Lz = 10 nm).

Figure 4 .

 4 The element size is chosen when the values of the components of the stiffness tensors C C changed by minus 0.1% when the

Figure 4 :

 4 Figure 4: 3D meshes with linear tetrahedron P1-element of RVE whit : a) flat graphene, b) wrinkled graphene, c) finite size graphene.

- 1 Flat10 5 .

 15 graphene without interphase (C I = C b ) Expression of Young's modulus as function of angle θ

S

  11 cos 2 (θ) + S 31 sin 2 (θ) S 12 cos 2 (θ) + S 23 sin 2 (θ) S 31 sin 2 (θ) + S 33 cos 2 (θ) 0 S 55 cos (θ) sin (θ

Table 1 :

 1 Symbols tableL x , L y , L zSide length of the RVE t

I Thickness of interphase n I

Unit normal vector to the interface I φ and φ s Bulk/interfacial free energy density u (x), u s (x)

Displacement field and surfacic displacement field in interface I u Displacement jump across the interface I t s , t

Average traction vector and traction vector jump across the interface I (x), s (x)

Infinitesimal strain tensor and infinitesimal surface strain tensor in interface σ (x), σ s (x)

Cauchy stress tensor in the bulk and surface Cauchy stress tensor in interface C b , C I

The fourth-order symmetric stiffness tensor for bulk polymer and interphase C s

The fourth-order symmetric interfacial stiffness tensors K s

Second-order symmetric cohesive interfacial stiffness tensor ¯ , σ

Effective infinitesimal strain tensor and effective Cauchy stress tensor C C , C A Effective stiffness tensor compute by continuum model and by atomistic model m α Masse of αth atom r α , r (0) α Current and initial position of αth atom r αβ , n αβ Inter-atomic distance and covalent bond direction between αth atom and βth atom v α Velocity of αth atom f αβ , f α = β f αβ Force of βth atom on the αth atom and force on αth atom U ({r α } , ¯ ) Potential energy of atomistic system for prescribe macro strain ¯ w (x), b w (x, r α , r β ) Weight function and bond function Murdoch-Hardy procedure ρ w (x)

Table 2 :

 2 Identified elastic constants of the continuum model in Voigt notation.

Table 3 :

 3 Influence on the effective properties of graphene modifications

	Stiffness Yield stress Toughness Electrical resistivity
	(E)	(σ y )	(K Ic )	(ρ e )
	Chemical modified
	graphene
	Geometrical modification	open	open
	(wrinkling, border)

  In this case, surface elastic parameter are express into the local basisB I [x] = (t 1 [x] , t 2 [x] , n I [x]) C s (x) =C s 11 t 1 [x] ⊗ t 1 [x] ⊗ t 1 [x] ⊗ t 1 [x](32)+ C s 22 t 2 [x] ⊗ t 2 [x] ⊗ t 2 [x] ⊗ t 2 [x] + C s 12 (t 1 [x] ⊗ t 1 [x] ⊗ t 2 [x] ⊗ t 2 [x] + t 2 [x] ⊗ t 2 [x] ⊗ t 1 [x] ⊗ t 1 [x]) + 2C s 66 (t 1 [x] ⊗ t 2 [x] + t 2 [x] ⊗ t 1 [x]) ⊗ (t 1 [x] ⊗ t 2 [x] + t 2 [x] ⊗ t 1 [x]) , K s (x) =K s 11 t 1 [x] ⊗ t 1 [x] + K s 22 t 2 [x] ⊗ t 2 [x] + K s 33 n

I [x] ⊗ n I [x] ,

Table 1 :

 1 Effective elastic constants in Voigt notation.

	Effective stiffness						Effective compliance
	Flat graphene
		 43.3 18.7 8.77	0	0	0				30.2	-10.6 -12.2	0	0	0
	CC =	      18.7 43.3 8.77 8.77 8.77 14.1 0 0 0 0 0 0	0 0 0.079 0	0 0 0 0.079	0 0 0 0	     	GPa S =	-10.6       -12.2 -12.30.2	-12.2	0	0	0
		0	0	0	0	0	12.3

  σ t S 11 cos 2 (θ) + S 31 sin 2 (θ) e 1 ⊗ e 1 + σ t S 12 cos 2 (θ) + S 23 sin 2 (θ) e 2 ⊗ e 2 + σ t S 31 sin 2 (θ) + S 33 cos 2 (θ) e 3 ⊗ e 3 + σ t S 55 cos (θ) sin (θ) 2 (e 3 ⊗ e 1 + e 1 ⊗ e 3 ) , and we deduce the normal strain t t = σ t S 11 cos 2 (θ) + S 31 sin 2 (θ) cos 2 (θ) + σ t S 31 sin 2 (θ) + S 33 cos 2 (θ) sin 2 (θ)

	+ σ t	S 55 2	(cos (θ) sin (θ))	2 .

  11 cos 4 (θ) + S 33 sin 4 (θ) + 2S 31 + S55

	2	sin 2 (θ) cos 2 (θ)	.	(42)

Preprint submitted to ElsevierMay 30, 2021

Supplementary Information for : Identification of elastic properties of interphase and interface in graphene-polymer nanocomposites by atomistic simulations

Abstract

This document provides supplementary information regarding:

1. Expressions of DREIDING potential in term of interatomic distance 2. Fitting details of the cohesive stiffness of interface K s 3. Computation of effective stiffness tensor by atomistic simulations 4. Computation of effective stiffness tensor by numerical homogenization 5. Expression of Young's modulus as function of angle θ

Expressions of DREIDING potential in term of interatomic distance

In the DREIDING potential [START_REF] Van Lier | Ab initio study of the elastic properties of single-walled carbon nanotubes and graphene[END_REF], the interactions of all the atoms are described by the valence (or bonded) interactions and the nonbonded interactions. For convienent reasons, we consider only the van der Waals interactions in the nonbonded part of DREIDING potential, which is modeled by the Lennard-Jones 12-6 type expression

where αβ and σ αβ are parameters which depends on the atom types. For computational reasons, we use a cut off radius r c = 14.0Å. We assume that the valence interactions consist of bond stretch (U (s) , two-body), bond-angle bend (U (b) , three-body), dihedral angle torsion (U (tor) , four-body). In the


J H Osório

F Amrani

F Delahaye

A Dhaybi

K Vasko

G Tessier

F Giovanardi

L Vincetti

B Debord

F Gérôme

F Benabid
email: f.benabid@xlim.fr

Hollow-core fibers with ultralow loss in the ultraviolet range and sub-thermodynamic equilibrium surfaceroughness

We report on hollow-core fibers displaying rms core roughness under the surface capillary waves thermodynamic limit and record loss in the short-wavelength range (50.0 dB/km at 290 nm, 9.7 dB/km at 369 nm, 5.0 dB/km at 480 nm, 0.9 dB/km at 558 nm, 1.8 dB/km at 719 nm).

Among the novel frontiers in photonics are ultraviolet (UV) photonics and next-generation optical fibers. The former holds a prominent position within technological and scientific areas, as it includes activities that entail or require the generation and control of light at wavelengths between 400 nm and 100 nm. Such activities encompass materials synthesis, micro-and nano-processing, precision spectroscopy, and light transport. The latter saw a huge progress, exemplified by the latest achievements in inhibited-coupling (IC) hollow-core photonic crystal fibers (HCPCFs) technology. Still, although the recent results on these fibers' performances are remarkable, further advances should rely on the improvement of the HCPCFs core surface quality as it imposes a fundamental hindrance to loss decreasing due to scattering processes. Differently from the IC HCPCFs guiding in the infrared (IR), whose loss remains limited by the fiber design (i.e., confinement loss, CL), the current state-of-the-art in HCPCFs with guidance at short wavelengths (λ < 800-1000 nm) is constrained by the surface scattering loss (SSL) [START_REF] Debord | Ultralow transmission loss in inhibited-coupling guiding hollow fibers[END_REF]. Thus, although the cladding design optimization efforts have entailed significant lessening of the loss for fibers guiding in the IR, reducing it in the visible and UV ranges lingers as a more challenging task owing to SSL limitation. SSL is set by the roughness of the core surface, which stems from frozen-in thermal surface capillary waves (SCW) that are forged during the fiber draw. SSL takes the form of Eq. ( 1), where F is the core mode and core contour overlap, 𝜆 is the wavelength, and 𝜆 0 is a calibration constant, and scales with the square of surface roughness root-mean-square (rms) height, i.e., 𝛼 𝑆𝑆𝐿 ∝ ℎ 𝑟𝑚𝑠 2 [START_REF] Roberts | Ultimate low loss of hollow-core photonic crystal fibres[END_REF]. The multiplying factor 𝜂 in Eq. (1) relates to the surface quality and is proportional to ℎ 𝑟𝑚𝑠 2 . The thermodynamic equilibrium surface roughness (TESR), ℎ 𝑇𝐸𝑆𝑅 , amounts to √𝑘 𝐵 𝑇 𝐺 𝛾 ⁄ ~ 0.4 nm for silica (kB: Boltzmann constant, TG: glass transition temperature; γ: interfacial tension).

𝛼 𝑆𝑆𝐿 = 𝜂 × 𝐹 × ( 𝜆 𝜆 0 ) -3 (1)
Thus, diminishing the SSL levels must rely on manipulating the surface roughness. In this context, recent results have shown that the roughness along the fiber drawing direction can be lowered by controlling the drawing stress [START_REF] Bresson | Anisotropic superattenuation of capillary waves on driven glass interfaces[END_REF]. Furthermore, SSL reduction enables the advent of optical fiber with unprecedentedly low transmission-loss and solarization-resistance, which, in turn, gives UV-photonics a strong leverage to further industrial and scientific applications.

Here, we report on the fabrication of HCPCFs with improved core surface quality and record low loss in the short-wavelength range. To accomplish it, we revisited the usual HCPCF fabrication methods and obtained a rms roughness reduction from ~0.4 nm in fibers drawn with standard fabrication methods to 0.15 nm in fibers fabricated by using an innovative technique. The new fibers have, therefore, sub-TESR levels. The reduction in the core roughness allowed attaining fibers with ultralow loss: 50.0 dB/km at 290 nm, 9.7 dB/km at 369 nm, 5.0 dB/km at 480 nm, 0.9 dB/km at 558 nm, 1.8 dB/km at 719 nm. We envisage that our investigation will confer a new panorama for fibers working in short-wavelengths and a new path for attaining ultra-smooth surfaces for next-generation applications in photonics. Fig. 1 shows the microstructures of the fibers we report herein (F#1 and F#2) and a zoom in the lowest loss regions. The fibers are tubular amorphous lattice cladding HCPCF formed by 8 non-touching tubes. Table 1 lists the core diameter, Dcore, the thickness and diameter of the cladding tubes t and Dtubes, and the spacing between the lattice tubes, g, of the fibers. F#1 exhibits minimum loss values of 50.0 dB/km at 290 nm, 9.7 dB/km at 369 nm, and 5.0 dB/km at 480 nm. The minimum attenuation values for F#2 are 0.9 dB/km at 558 nm and 1.8 dB/km at 719 nm. Fig. 2 displays the measured loss of F#1 and F#2, the silica Rayleigh scattering limit (SRSL) trend, and a typical output near field profile. Noteworthily, the values reported herein are new record low loss figures for fibers working in the short-wavelength range [START_REF] Chafer | 1 km hollowcore fiber with loss at the silica Rayleigh limit in the green spectral region[END_REF][START_REF] Gao | Conquering the Rayleigh scattering limit of silica glass fiber at visible wavelengths with a hollow-core fiber approach[END_REF][START_REF] Sakr | Hollow core optical fibres with comparable attenuations to silica fibres between 600 nm and 1100 nm[END_REF][START_REF] Gao | Nodeless hollow-core fiber for the visible spectral range[END_REF][START_REF] Yu | Single-mode solarization-free hollow-core fiber for ultraviolet pulse delivery[END_REF] and lie below the Rayleigh scattering limit of bulk silica (in red). The remarkable ultralow loss reported above could be attained thanks to an amelioration of the core surface quality, which was accomplished by controlling the shear stress acting onto the cladding microstructure during the draw [START_REF] Roberts | Ultimate low loss of hollow-core photonic crystal fibres[END_REF]. To study the effectiveness of the new fabrication approach, we assessed, by using a picometer resolution profilometer [START_REF] Brun | Picometer-scale surface roughness measurements inside hollow glass fibres[END_REF], the core roughness in two groups of tubular amorphous lattice cladding HCPCFs, G#1 and G#2. G#1 fibers have been drawn by using the standard procedure and G#2 fibers have been drawn by using our new methods. Fig. 3a exhibits typical surface profiles along the HCPCF axis of G#1 and G#2 fibers. Fig. 3a apprises that the peak-to-peak roughness of G#1 fibers has values around 1.5 nm while, for G#2 fibers, values around 0.5 nm. Fig. 3b displays the rms roughness of G#1 and G#2 fibers with different cladding tube thicknesses estimated from several surface roughness measurements. The rms roughness of G#1 fibers fluctuates around 0.4 nm. G#1 fibers have, thus, TESR. In turn, G#2 fibers have rms roughness around 0.15 nm, which is considerably lower than the TESR limit set by the SCW scenario in amorphous silica (Fig. 3b, blue line). Fig. 3c shows the typical loss for G#1 and G#2 fibers with similar tubes thickness. The distinction between the loss trends is patent. The loss in the reported G#1 fiber (Dcore = 41 μm) presents a considerable loss increase for λ < 750 nm due to SSL influence. In the reported G#2 fiber (Dcore = 27 μm), a decreasing loss trend is observed even for λ < 750 nm. Fig. 3c also shows the CL, SSL, and total loss (TL = CL + SSL) trends calculated by considering the scaling laws in [START_REF] Vincetti | Empirical formulas for calculating loss in hollow core tube lattice fibers[END_REF]. Here, we used η = 2133, corresponding to a fiber surface with TESR (in blue), and η = 300, corresponding to a surface-roughness rms that is 2.7 lower than the TESR rms (in green). The η values have been estimated by adjusting its value to the loss of fibers displaying TESR and sub-TESR [START_REF] Roberts | Ultimate low loss of hollow-core photonic crystal fibres[END_REF][START_REF] Bresson | Anisotropic superattenuation of capillary waves on driven glass interfaces[END_REF][START_REF] Chen | Multi-kilometer long, longitudinally uniform hollow core photonic bandgap fibers for broadband low latency data transmission[END_REF]. It is seen that the TL turning point happens at shorter wavelengths for smaller SSL. The difference between the loss behaviors of G#1 and G#2 fibers is, hence, a signature of SSL reduction in the new fibers.

In summary, we developed HCPCFs exhibiting ultralow loss at short wavelengths by improving the fiber cores' surface quality and by attaining sub-TESR levels. We envision that our research will provide a new scenery for the next HCPCFs for the visible and ultraviolet ranges and a new avenue for attaining nextgeneration photonic devices with ultralow surface roughness levels.

Fig. 1 .

 1 Fig. 1. F#1 and F#2 cross-sections and measured loss. Minimum loss values are 50.0 dB/km at 290 nm, 9.7 dB/km at 369 nm, 5.0 dB/km at 480 nm, 0.9 dB/km at 558 nm, 1.8 dB/km at 719 nm.

Fig. 2 .

 2 Fig. 2. F#1 and F#2 measured loss together with the silica Rayleigh scattering limit (SRSL) trend and a typical output near field profile.

Fig. 3 .

 3 Fig. 3. (a) Typical core surface roughness profiles (Δh: roughness height variation) and (b) measured rms roughness. (c) Typical loss spectra for fibers in G#1 and G#2; Confinement loss (CL), surface scattering loss (SSL) and total loss (TL) trends.

Table 1 .

 1 F#1 and F#2 geometrical parameters (Dcore: core diameter; t: thickness of the cladding tubes; Dtubes: diameter of the cladding tubes; g: spacing between the lattice tubes).

		F#1	F#2
	Dcore (μm)	27	42
	t (μm)	0.6	0.9
	g (μm)	2.1-4.7	3.6-5.2
	Dtubes (μm)	11

Acknowledgements

This research was funded through PIA program (grant 4F) and la region de la Nouvelle Aquitaine.


est destinée au dépôt

Learning of the scientific approach at university : the case of research situations from problems of discrete mathematics

Mickael da Ronch, Michèle Gandit, Sylvain Gravier

To cite this version:

Mickael da Ronch, Michèle Gandit, Sylvain Gravier. Learning of the scientific approach at university : the case of research situations from problems of discrete mathematics. 14th International Congress on Mathematical Education, Jul 2021, Shanghai, China. , 2021, Proceedings of the ICME-14 Topic Study Group 15 : Teaching and learning of discrete mathematics. hal-03285073


W Gomes

P Maurice

N Settembre

J Theurel

L Wioland

J J Atain-Kouadio

L Claudon

B Chenuel

H Hani

A Kimmoun

B Levy

J M Sellal

J Paysant

S Malikov

S Ivaldi

Biomechanical effects of using a passive back support exoskeleton during prone-positioning maneuver: A pilot study

Keywords: Exoskeleton, Ergonomics, Workload Assessment, Digital Human Simulation, Healthcare Sector

Introduction

Musculoskeletal disorders associated with patients care and manipulation are frequent among healthcare workers, due to force exertion and sustained awkward postures (EU Commission 2011[START_REF] Davis | Prevalence of musculoskeletal disorders for nurses in hospitals, long-term care facilities, and home health care: A comprehensive review[END_REF]. During the COVID-19 pandemic, intensive care units (ICUs) staffs had to perform prone-positioning (PP) of a large number of intubated patients suffering from Acute Respiratory Distress Syndrome (ARDS) many times a day. PP maneuver consists in repositioning a patient from a supine position (on the back) to a prone position (on the front side), which improves oxygenation and ventilatory mechanisms in ARDS patients [START_REF] Guérin | Prone positioning in severe acute respiratory distress syndrome[END_REF]. PP is a delicate procedure (risk of extubation) which requires the medical staff to remain bent forward for several minutes (Fig. 1). This posture result in back loading and in the long-term potential injuries in the low-back. Similar postures are common in industrial tasks, and recently occupational exoskeletons have started to be commercialized to help reduce physical load on workers (DeLooze et al. 2016). Exoskeletons therefore appear as a promising solution to support healthcare workers, though few studies focused on this sector (e.g. [START_REF] Miura | Successful Use of the hybrid assistive limb for care support to reduce lumbar load in a simulated patient transfer[END_REF][START_REF] Iishi | Development of a new type of lightweight power assist suit for transfer work[END_REF]. But while studies overall agree that back-support exoskeletons, and specifically passive devices, have the potential to reduce lumbar muscle activity in tasks involving trunk bending, the exact efficacy and usability of such devices appear largely device-and task-dependent [START_REF] Theurel | Occupational exoskeletons: Overview of their benefits and limitations in preventing work-related musculoskeletal disorders[END_REF]. A preliminary study of our team [START_REF] Settembre | The use of exoskeletons to help with prone positionning in the intensive care unit during COVID-19[END_REF] identified the Laevo passive exoskeleton for lumbar support (Laevo, the Netherlands) as the exoskeleton that best matched the PP maneuver and ICU constraints, among 4 exoskeletons that were available at the time of the test during the COVID-19 first wave. The present study presents a pilot investigation of the biomechanical effects of using the Laevo during PP maneuver. Two healthy volunteers with PP experience performed PP maneuvers both with and without wearing the Laevo, on a patient simulator (manikin) at the Hospital Simulation Center of the University of Lorraine (Fig. 1). Each maneuver lasted about 4 min. In a first experiment, the participants were equipped with the Xsens MVN Link inertial motion capture system (capture rate 240 Hz) to record whole-body kinematics. In a subsequent experiment, the participants were equipped with a Delsys Trigno EKG sensor to measure heart rate, and 12 Delsys Trigno sEMG sensors (sampling rate 4370 sa/sec) placed according to the Seniam recommendations bilaterally on the erector spinae longissimus (ESL), erector spinae iliocostalis (ESI), trapezius ascendens (TA), and biceps femoris long head (BF), and laterally on the right rectus abdominis at T10 level (RA), rectus femoris (RF), gluteus maximus (GM) and tibialis anterior (TAL).

Methods

Experimental protocol

Data analysis

In order to analyze the kinematics and dynamics of the maneuver, we replayed the recorded motions in a physicsengine based simulation using a 43 degrees of freedom digital human model (DHM). The recorded motions were retargeted on the DHM directly in the Cartesian space by using a quadratic programming controller (QP) typically used for humanoid robots [START_REF] Maurice | Assessing and improving human movements using sensitivity analysis and digital human simulation[END_REF]. The QP control directly estimates the joint torques associated with the motion. This method has the advantage to blend the traditional inverse kinematics and inverse dynamics steps into one faster step, since the computation runs in close to real-time. Since we did not have a direct measurement of the support torque provided by the Laevo, we used the empirical calibration curve published by [START_REF] Koopman | Effects of a passive exoskeleton on the mechanical loading of the low back in static holding tasks[END_REF] to estimate the lumbar support torque depending on the trunk flexion angle. We used the resulting L5/S1 flexion/extension estimated joint torque as a complement to the physiological measurements obtained with the sEMG.

Results and discussion

Kinematic results: Using the Laevo did not significantly modify the kinematics of the motion during PP. Specifically, the profile of the lumbar flexion angle remained overall similar (Fig. 2). This result was confirmed by the subjective feedback of the participants who reported that they did not feel the exoskeleton was affecting the way they performed the task. Dynamic results: The median value of the estimated L5/S1 flexion/extension torque decreased by about 12% when using the Laevo (Fig. 2). This result agrees with the reduced perceived lumbar effort reported by participants, and is aligned with the 15% reduction observed by [START_REF] Koopman | Effects of a passive exoskeleton on the mechanical loading of the low back in static holding tasks[END_REF]) in a similar static forward bending task. Yet those results should be considered cautiously, since only the assistive torque of the exoskeleton was modeled, and not the full effect of the load transfer. Physiological results: We did not observe any major change in heart rate when using the Laevo. This might be due to the short duration of the experimental task (4 min). The activity of the back extensor muscles (except ESI) and hip extensor muscles overall decreased when using the Laevo (Table 1). The other recorded muscles remained mostly unaffected by the exoskeleton. Those results are consistent with the perceived reduced effort reported by the participants.

Conclusions

Our pilot study, though limited to a small number of participants, suggests that using the Laevo could help reduce the musculoskeletal load on the low back during PP maneuver, without causing any significant negative side-effect nor modifying the practice. Following this study, 2 Laevos have been successfully deployed in the ICU of CHRU Nancy and have been used there for the past 6 months. Recording quantitative measurements during this field-testing is not possible. However subjective feedback from users is collected through questionnaires and interviews, and will be analyzed in a future study. during PP maneuver with and without using the Laevo.

Figure 1 :

 1 Figure 1: Medical staff equipped with the Laevo exoskeleton performing a PP maneuver on a manikin.

Figure 2 :

 2 Figure 2: Time-series (top) and distribution across time (bottom) of the L5/S1 flexion angle (left) and torque (right)during PP maneuver with and without using the Laevo.

Table 1 :

 1 Relative difference (%) in RMS value of the EMG signal with the Laevo vs. without it, for both participants and both PP positions (patient's head (H) and side (S)).

		TA	ESL	ESI	GM	RF	BF RA TAL
	P1 H -35	-12	-1	-30	0	-16	-1	-2
	P2 H	+2	+5	0	+5	NA -20	0	+6
	P1 S -15	-9	-1	-12	-6	-53	0	-10
	P2 S	-5	-11	-2	+8	NA -34	0	-20


On the Schwarz Symmetry Principle in a Model Case

Joël Merker, Francine Meylan

To cite this version:

Joël Merker, Francine Meylan. On the Schwarz Symmetry Principle in a Model Case. Proceedings of the American Mathematical Society, 1999, 127 (4), pp.1097-1102. hal-03285202