Eduard Marenić

Adnan Ibrahimbegovic

Homogenized elastic properties of graphene for moderate deformations

Keywords: graphene, elastic properties, molecular mechanics, homogenization

This paper presents a simple procedure of obtaining a substitute, homogenized mechanical response of single layer graphene sheet. The procedure is based on the judicious combination of molecular mechanics simulation results and homogenization method. Moreover, a series of virtual experiments are performed on the representative graphene lattice. Following these results, the constitutive model development is based on the well established continuum mechanics framework, that is, the non-linear membrane theory which includes the hyperelastic model in terms of principal stretches. A proof-of-concept and performance is shown on a simple model problem where the hyperelastic strain energy density function is chosen in polynomial form.

Introduction

Nanoscale systems and processes based on graphene [START_REF] Novoselov | Two-dimensional atomic crystals[END_REF], Singh et al., 2011] are becoming more viable for engineering applications, however, our ability to model their performance remains limited. The main challenge is that parts of graphene-based devices modelled with discreet models (such as molecular mechanics (MM) and molecular dynamics (MD) ) typically contain extremely large number of particles, even though the actual physical dimension may be quite small. For instance, a simple square shaped Single Layer Graphene Sheet (SLGS) with the side length of approximately 500 nm has already nearly one million of carbon atoms. Following explicitly the trajectory for large number of degrees of freedom easily becomes intractable. Thus we reach for a substitute, continuum model which models the average mechanical behaviour of atomic system.

Averaged continuum properties of graphene in the context of infinitesimal deformation is the subject of research for nearly 10 past years (e.g. see [START_REF] Caillerie | Discrete homogenization in graphene sheet modeling[END_REF], Reddy et al., 2006, Huang et al., 2006, Zhao et al., 2009, Georgantzinos et al., 2010]). However, there is a large discrepancy in the results obtained by means of the different simulation methods and experimental studies (e.g. see [START_REF] Reddy | Equilibrium configuration and continuum elastic properties of finite sized graphene[END_REF],Huang et al., 2006]). We reported also in [START_REF] Marenić | Homogenized elastic properties of graphene for small deformations[END_REF] an overview of the mechanisms causing the result scattering which can be summarized as: formulation differences (which concerns MM, MD, ab initio methods, continuum mechanics); choice of interatomic potential driving the atomic system (Tersoff-Brenner [Tersoff, 1986, Brenner, 1990] usually, also Morse [Morse, 1929], AMBER and second generation REBO [START_REF] Brenner | A second-generation reactive empirical bond order (rebo) potential energy expression for hydrocarbons[END_REF]); uncertainty of the thickness (yet known as Yakobson's paradox) [START_REF] Huang | Thickness of graphene and single-wall carbon nanotubes[END_REF]; and also size effect, relaxation (minimisation of the energy due to coordination), chirality, and edge passivation [START_REF] Lu | Atomistic simulation and continuum modeling of graphene nanoribbons under uniaxial tension[END_REF]. In [START_REF] Marenić | Homogenized elastic properties of graphene for small deformations[END_REF], we also pointed out the influence of the boundary conditions on the (linear-elastic) stiffness of SLGS, which is related to the homogenized stiffness bounds as introduced in [Huet, 1990] and latter discussed in [START_REF] Markovic | On micro-macro interface conditions for micro scale based FEM for inelastic behavior of heterogeneous materials[END_REF].

We tend here to give a combined MM and continuum approach for the study of the in plane mechanical behaviour of SLGS under large deformations. In other words we present equivalent continuum modelling of large deformations of graphene, which goes beyond what linear theory can handle. We note here that continuum models of nano-structures are extremely simplified, however they bring insight as well as quantitative information on the relevant physical phenomena . Thus, these models seam optimal for the current and potential future applications of graphene materials, such as, e.g., reinforcement agents to strengthen composites or structural parts. Recently, graphene has attracted both academic and industrial interest in graphene-based polymer nanocomposites because it can produce a dramatic improvement in properties at very low filler content. Thus it can be used either in graphene-based nano-electro-mechanical systems (NEMS) devices [START_REF] Geim | The rise of graphene[END_REF], or graphene-conducting polymer nanocomposites (see [START_REF] Hu | Graphene-polymer nanocomposites for structural and functional applications[END_REF], Stankovich et al., 2006, Gmez et al., 2011]. Micro-mechanical models are conveniently used to predict macroscopic properties of conventional reinforced polymer composites, which behave as heterogeneous material (Mori and Tanaka, 1973) on the scale somewhat higher then the nano length scale. However, when the multi-scale model has a goal to predict properties of nano-reinforced composites, a detailed knowledge of the response on the scale of the nano-structure is required (see e.g. an example of the multiscale approach for the nanoclay reinforced polymer [Bdoui andCauvin, 2012, Gelineau et al., 2015]. Similarly for the graphene-based nanocomposites, and mentioned application, the mechanical response of graphene under different loading programs, boundary conditions and large strain regime is crucial, and should still be better understood. Experimental investigation (see e.g. [START_REF] Lee | Measurement of the elastic and intrinsic strength of monolayer graphene[END_REF]) and first-principles calculations show that graphene undergo very large deformations, with highly nonlinear behavior and still remain elastic, with stable bonds and intact bond topology. This fact allows us to adapt and use a finite strain elastic model for the large range of phenomena and motivates us to benefit from somewhat simplified nonlinear membrane theory. The main difficulty in constructing the solutions for the problems of the deformable body subjected to large deformations concerns the choice of the suitable formulation like Lagrangian, Eulerian, updated Lagrangian, etc. Choosing a particular formulation implies the corresponding choice of the reference frame and the properly invariant measure of large strain [START_REF] Ibrahimbegovic | Nonlinear Solid Mechanics[END_REF], Wriggers, 2008, Belytschko et al., 2000]. There are generally two ways of dealing with this problem. First considers the theoretical development placed in the framework of differential manifolds [START_REF] Abraham | Foundations of Mechanics[END_REF] which enables the switch between configuration by means of the metric tensor. However, the conceptual clarity of the theoretical formulation of the finite deformation elasticity set on a differential manifold does not imply simplicity with respect to numerical implementation, mainly due to locking. In [START_REF] Arroyo | An atomistic-based finite deformation membrane for single layer crystalline films[END_REF] a general methodology to develop hyper-elastic membrane models for single-atomic-layered films is presented with an extension of the Cauchy-Born (CB) rule [Ericksen, 1984, Zanzotto, 1996, Ericksen, 2008] based on the exponential map. The exponential map is added in the formulation as an extension of CB rule to account for the curvature of the lattice vectors. Thus, they are using curvilinear coordinates and giving general formulation capable to treat even buckling problems of carbon nano tubes.

In this paper we are performing virtual experiments (MM simulation) on the representative sample of the graphene lattice and thus we are not making any assumption regarding the kinematics of the atoms (CB hypothesis does not have to be fulfilled). These virtual experiments are used to fit energy potential. Similar procedure is performed in e.g. [START_REF] Reddy | Equilibrium configuration and continuum elastic properties of finite sized graphene[END_REF], Lu and Huang, 2009, Cadelano et al., 2009, Lu et al., 2011] where only uniaxial tests are performed giving an elastic potential with respect to nominal strain. In [START_REF] Cadelano | Nonlinear elasticity of monolayer graphene[END_REF] continuum elasticity theory and tight-binding atomistic simulations are combined to determine the constitutive nonlinear stress-strain relation, and the corresponding nonlinear elastic moduli for graphene. We show here biaxial tests and, moreover, we seek to adopt the nonlinear membrane theory which includes, as a special case, the hyperelastic model in terms of principal stretches. The latter was often used to characterise rubberlike materials, see [START_REF] Ibrahimbegović | A consistent finite element formulation of nonlinear membrane shell theory with particular reference to elastic rubberlike material[END_REF], with finite element formulation of elastic membrane shells with co-axial energy-conjugate pairs of stress and strain measures. The main advantage of the formulation in [START_REF] Ibrahimbegović | A consistent finite element formulation of nonlinear membrane shell theory with particular reference to elastic rubberlike material[END_REF] pertains to the finite element analysis of elastic membranes. The problem of the finite deformation of elastic membranes can be solved relaying on the Cartesian structure without going to a more general setting of manifolds (note that this is not possible for the elastoplastic membranes [Ibrahimbegović, 1994]). As a first step, we present the modelling of the in plane large deformation of SLGS , however the proposed theory can be extended to membranes whose reference and current configurations can be arbitrary space-curved surfaces. Thus, the development of homogenized constitutive model is based on the well established continuum mechanics framework. The main novelty concerns the specific application to graphene with the finite element implementation being straightforward when relying upon previous works on large deformation model for rubber-like materials. The series of virtual tests can be costly, but needs to be performed only once. Moreover, the developed model is fully capable to reproduce the linear elastic behaviour in small strain regime as well as the nonlinear that occurs in large strain regime.

The outline of the paper is as follows. We will address the problem of analysis of finite deformation of general, space-curved, elastic membranes in the next section. In the section 2 we present the virtual experiments performed on the SLGS and the development of the substitute continuum model in terms of principal stretches. In the section 3 we elaborate on the performance of the developed model, and give a closing remarks and future research in the section 5.

2 Continuum elastic membrane in finite deformations

Generalities and motivation

We will first briefly revisit continuum model problem in large displacements and solution strategy. Let the position of each point be denoted with X and x = ϕ(X) in reference (Ω) and current (Ω ϕ ) configuration, respectively, where ϕ(•) denotes the motion as a point transformation. For each point X we define the displacement vector u(X) = x -X. When the SLGS is submitted to large deformations the difference between the initial configuration at the beginning of the load program and the deformed configuration, can no longer be ignored as for the case of small deformations, characterized by small strain tensor, . There is a large variety (theoretically infinite) of possible choices for stress and strain tensors available for the continuum large strain problem formulation as presented in most of the textbooks covering the subject (see e.g. [START_REF] Ibrahimbegovic | Nonlinear Solid Mechanics[END_REF],Wriggers, 2008]). For the continuum large strain problem the constitutive model formulation depends on the choice of the particular work-conjugate pair. Usually first Piola-Kirchhoff stress and deformation gradient (P , F ) or second Piola-Kirchhoff stress and Green-Lagrange strain (S, E) are chosen to express internal work. We now construct the weak form of the continuum boundary value problem in Ω for the case of large displacements. This immediately introduces a solution strategy by weakening the way of satisfying the equilibrium, i.e. it is satisfied only in average sense. Therefore, we assume that Dirichlet boundary conditions u = ū are prescribed on the part Γ u of the boundary Γ. The nanostructure system treated as surrogate continuum is subjected to tractions t on the part Γ σ of the boundary and to a volume forces b in Ω. We choose a virtual displacement field v as infinitesimal and kinematically admissible with respect to Dirichlet boundary conditions, thus each component v i takes a zero value on the Γ u i.e.

V 0 := {v i : Ω → R | [v i ] Γu i = 0}.
We also suppose that the virtual displacement is supperposed on the deformed configuration and par ametrized by the coordinates in the deformed configuration (Ω ϕ ). For the real displacement vector field u the components u i are defined within

V := {u i : Ω → R | [u i ] Γu i = ūi }.
The weak form of equilibrium at large displacements in material description (Ω) states

0 = G(ϕ; v) := Ω Γ • SdΩ - Ω v • bdΩ - Γσ v • tdΓ, ( 1
)
where Γ is the virtual Green-Lagrange strain given as the directional derivative of the Green-Lagrange strain measure.

The constitutive model formulation depends on the choice of the particular work-conjugate pair, however the unique form of the constitutive relation can be written for the hyperelastic material model in terms of the strain energy potential, W (•). Energy potential is an unique choice, since all possibilities of the stress and strain measures are only different material representations of the same work. The generalised approach for establishing the well-posed form of the strain energy is given in terms of the polyconvexity conditions. The role of the polyconvexity conditions is to ensure that the large strain remain accompanied by large stress. Moreover, the polyconvexity conditions impose that the strain energy remains a convex function which can be written as

W (αa 1 + (1 -α)a 2 ) ≤ αW (a 1 ) + (1 -α)W (a 2 ), (2)
where a 1 and a 2 are functions representing certain intrinsic measure of deformation and 0 < α < 1. Besides the polyconvexity conditions applicable only to hyperelastic materials, there is a number of invariance restrictions on the general elastic response, which any constitutive model ought to respect. In order to describe the elastic response that satisfies the invariance requirements, the strain energy potential is usually expressed as a function of principal invariants

W (i 1C , i 2C , i 3C ), (3)
where the principal invariants of the right (C = F T F ) and left (B = F F T ) Cauchy-Green deformation tensors are given as

i 1C := tr[F T F ] ≡ tr[F F T ] =: i 1B i 2C := 1 2 (tr[C]) 2 -tr[C 2 ] ≡ 1 2 (tr[B]) 2 -tr[B 2 ] =: i 2B (4) i 3C := det[F T F ] = (det[F ]) 2 ≡ det[F F T ] =: i 3B .

Constitutive law in terms of prinipal stretches for large deformation of graphene

An elegant alternative to (3) is the strain energy potential defined in terms of the principal stretches λ i , i = 1, 2, 3. These values correspond to the principal values of the stretch tensors, right U or left V . This derives from the standard eigenvalue problem which can be written either in material (a) or in spatial (b) description

(a) (U -λ i I)n i = 0, (b) (V -λ i I ϕ )m i = 0. (5)
Note that the computed principal (eigen) values λ i from the equation ( 5) remain the same in both descriptions, but the corresponding eigenvectors n i and m i change. By solving these eigenvalue problems, we can obtain spectral decomposition of the deformation gradient, rotation tensor

(a) F = 3 i=1 λ i m i ⊗ n i , (b) R = 3 i=1 m i ⊗ n i , (6)
and both stretch tensors

(a) U = 3 i=1 λ i n i ⊗ n i , (b) V = 3 i=1 λ i m i ⊗ m i . (7)
Note that these results hold for the principal vectors that form the ortho-normal principal frames, i.e., n i • In j = δ ij , where δ ij is the Kronecker delta. We further discuss the spectral decomposition of the Cauchy-Green tensors, related to the choice of the class of constitutive equations. Considering (6a) and it's transpose

F T = 3 i=1 λ i n i ⊗ m i , (8)
the spectral decomposition for both Cauchy-Green tensors is given as

(a) C = 3 i=1 λ 2 i n i ⊗ n i , (b) B = 3 i=1 λ 2 i m i ⊗ m i . (9)
With these results in hand we can easily express the principal invariants from (4) in terms of the principal stretches

i 1C := λ 2 1 + λ 2 2 + λ 2 3 , i 2C := λ 2 1 λ 2 2 + λ 2 2 λ 2 3 + λ 2 3 λ 2 1 , i 2C := λ 2 1 λ 2 2 λ 2 3 . (10)
Thus, any isotropic hyperelastic response that satisfies material invariance restriction can be expressed in terms of strain energy potential as a function of principal stretches. Thus, starting from the strain energy potential written as a function of principle invariants (3), we can express using (10) strain energy as a function of principal stretches

W (i 1C , i 2C , i 3C ) → W (λ 1 , λ 2 , λ 3 ). ( 11
)
For the elastic material behaviour the unit principal stretch corresponds to the case of no deformation which is accompanied by the zero value of strain energy:

λ i → 1 ⇒ W (λ i ) → 0. ( 12
)
Formulating the strain energy potential as in ( 11) makes it simple to check the polyconvexity conditions described above. These conditions enforce that large stresses should accompany large values of strains which is written in terms of principal stretches as

W (λ i ) → ∞ if {λ 1 , λ 2 , λ 3 } → ∞ (in tension), W (λ i ) → ∞ if {λ 1 , λ 2 , λ 3 } → 0 + (in compression). (13)
The last result states that polyconvexity conditions require the strain energy convexity with respect to each principal stretch. Non-convexity is responsible for many phenomena like the development of dislocations and phase changes, however we are not dealing with these issues in this work.

We turn now to 2D case formulation that describes the SLGS sheet. Thus, putting conveniently graphene sheet in the plane x 1 x 2 , we are neglecting the out of plane stretch λ 3 , following the usual hypothesis for the membrane theory. Since the sheet is made of a single atomic layer and the out-of-plane strains and stresses might be difficult to interpret, we restrict the considerations to the in-plane 2D continuum mechanics (this assumption is made by majority of authors dealing with graphene, see e.g. [Volokh, 2012]). Considering the mentioned assumption the strain energy density (SED) in (11) becomes W (λ 1 , λ 2 ). We further present the procedure to calculate second Piola-Kirchhoff stress (S) and tangent elasticity (C) tensors from the strain potential written in terms of the principal stretches. This computation is still performed in the conventional manner as

(a) S = ∂W ∂E = 2 ∂W ∂C , (b) C = ∂S ∂E = 2 ∂S ∂C . (14)
However, the computation of the stress and tangent elasticity tensors from the material model given by W (λ 1 , λ 2 ) is not performed directly, but rather using a simple chain rule. Thus, an important role is played by the auxiliary result pertaining to derivatives of the principal values λ i . This result can be obtained by applying the Gâteaux derivative formalism to the corresponding eigenvalue problem leading to

∂λ i ∂C = 1 2λ i n i ⊗ n i . ( 15
)
With this result in hand, we can calculate the second Piola-Kirchhoff stress tensor from the SED potential written in terms of principal stretches

S = 2 i=1 1 λ i ∂W (λ 1 , λ 2 ) ∂λ i n i ⊗ n i . ( 16
)
Comparing equations ( 16) and (9a) one can notice that the second Piola-Kirchhoff stress tensor is coaxial with the right Cauchy-Green tensor. Existence of the the coaxial energy-conjugate pair enables the solution of the problem of finite deformation membrane relaying on the Cartesian structure. The second Piola-Kirchhoff stress tensor can be further decomposed as

S = 2 i=1 s i n i ⊗ n i , (17)
where the n i represent the same eigenvectors as in first expression in (9), and the term s i denoting the principal stresses can be written as

s i = 1 λ i ∂W (λ 1 , λ 2 ) ∂λ i . ( 18
)
Next, we turn to the calculation of the elastic tangent modulus. This is done in the same manner, i.e., by performing a next step of directional derivative computation, which gives

C = 2 2 i=1 ∂s i ∂C n i ⊗ n i + 2 2 i=1 s i ∂ ∂C (n i ⊗ n i ). ( 19
)
The first and the second terms on the right hand side in ( 19) correspond to material (C mat ) and geometric (C geo ) part of the tangent elasticity tensor, respectively. Using the auxiliary result in ( 15) we obtain the expression for the material part of the tangent elasticity tensor:

C mat = 2 i=1 2 j=1 1 λ j ∂s i ∂λ j D ij [n i ⊗ n i ][n j ⊗ n j ]. (20)
Note that the material part of the tangent elasticity tensor is usually given in terms of its reduced form D ij in principal axes, see [START_REF] Ibrahimbegović | A consistent finite element formulation of nonlinear membrane shell theory with particular reference to elastic rubberlike material[END_REF]. The derivation of explicit form of the geometric part of the tangent elasticity tensor starts from the spectral decomposition of the right Cauchy-Green strain tensor (9) and considers a systematic usage of the auxiliary result in (15). Due to brevity we omitted this derivation and we give the final expression of the elastic tangent modulus in tensor notation

C := 2 i=1 2 j=1 D ij [n i ⊗ n i ][n j ⊗ n j ] + 2 s 1 -s 2 λ 2 1 -λ 2 2 [I -(n 1 ⊗ n 1 ) ⊗ (n 1 ⊗ n 1 ) -(n 2 ⊗ n 2 ) ⊗ (n 2 ⊗ n 2 )] ,
(21) where I = 1 2 (δ ik δ jl + δ il δ jk ). The constitutive law in terms of prinipal stretches for large deformation of membrane (used here for the SLGS), namely the expression ( 18) and ( 21), are used directly in the finite element solution procedure. The details about the 2D elastic membrane finite element can be found in most of the books dealing with nonlinear solid mechanics, and will not be discussed herein. For the completeness, we summarize the main steps needed in the finite element approximation, which for the presented constitutive law in terms of prinipal stretches boils down simply to the matrix form of the results obtained for principal stresses and elasticity, see Appendix A.

3 Development of constitutive law from virtual experiments on molecular model of graphene

Molecular model

In Section 2 a continuum model problem in large displacements is presented. Likewise, we introduce here briefly the molecular model problem. We consider a reference domain Ω a in R 3 which is occupied by N atoms placed within graphene nanostructure. A major feature of the structure of graphene is the hexagon pattern that repeats itself periodically in space. As a result of the periodicity, each atom is bonded to three neighbouring atoms. Such structure is mainly due to the process of sp 2 hybridization during which one s-orbital and two p-orbitals combine to form three hybrid sp 2 -orbitals at 120 • to each other within a plane [Atkins andDe Paula, 2006, Ruoff et al., 2003]. This covalent bond, often referred to as the σ-bond, is a strong chemical bond and plays an important role in the impressive mechanical properties of graphene, while the out-of-plane bond (the π-bond) that is relatively weak contributes to the interaction between the layers of graphene. Let X i and x i denote, respectively, the position vectors in the reference and the current configurations of atom i, where i = 1, . . . , N . The corresponding displacement vector of atom i is given by d i = x i -X i . The boundary conditions are imposed atom-wise in a quasi-static manner, such that either the displacement, di , or the external point force, fi , is given. The total energy (E tot ) stored in the atomic structure is given by

E tot = U (x 1 , . . . , x N ) - N i fi • d i , (22)
where U denotes the energy stored in the atomic bonds, as presented in sequel, and the second term on the right-hand side represents the external energy. The state of equilibrium of the atomistic system requires the variation of the total energy to be equal to zero

δE tot = N i ∂U ∂x i -fi • δx i = 0, (23)
where δx i represents the kinematically admissible virtual motion. Linearizing (23) and writing the result in matrix notation leads to

K (k) ∆d (k) = F (k) , (24)
where ∆d (k) is the displacement increment corresponding to the k-th load increment, whereas K (k) and F (k) are global stiffness and the residual vector, respectively. The latter can explicitly be defined as:

K ij = ∂ 2 U ∂x i ∂x j , F i = ∂U ∂x i -fi . ( 25
)
Unlike conventional FE method for continuum mechanics, we derive and assemble the stiffness and residual matrices by looping over all atoms. Due to the non-linear nature of the interatomic potentials and geometrically nonlinear kinematics, an incremental-iterative solver is used. For each load increment, several Newton iterations are performed, until convergence criteria are met. At each iteration (k) the atomic positions are updated as x

(k+1) i = x (k) i + ∆d (k) . ( 26
)
The initial iteration (k) = 0 starts at the initial configuration of the atomic system, with the position vector, x

i = X i .

Looking at the equation ( 25) it is obvious that the heart of molecular model lays in the internal energy governed by the interatomic potential. For the molecular simulation of a graphene sheet, a Morse potential [Morse, 1929] is used, more precisely we use it's modified version following [START_REF] Belytschko | Atomistic simulations of nanotube fracture[END_REF]. The modified Morse potential is a sum of pair and angular parts given as

U = bonds V M p (r) + angles V M θ (θ). (27)
The first term is a function of the chosen atom distance r to its first neighbour, whereas the second depends upon angle θ between particular atom bonds (see Fig. 1 b) in sequel). For the modified Morse potential, the bond energy terms are given as

V M p (r) = D e (1 -e -β(r-r 0 ) ) 2 -1 , (28)
V M θ (θ) = 1 2 k θ (θ -θ 0 ) 2 [1 + k sext (θ -θ 0 ) 4 ], (29)
where the constants of the potential according to [START_REF] Belytschko | Atomistic simulations of nanotube fracture[END_REF] are: D e = 6.03105 × 10 -19 Nm, β = 2.625×10 10 m -1 , k θ = 0.9×10 -18 Nm rad -2 , k sext = 0.754 rad -4 , the initial value of the bond length r 0 = 1.39 × 10 -10 m and the bond angle θ 0 = 2π/3 rad. The numerical implementation of MM model based on modified Morse potential, i.e., the forming of the residual force and tangent stiffness matrices, and the assembly procedure are given in more detail in [START_REF] Marenić | Homogenized elastic properties of graphene for small deformations[END_REF]. A general procedure for carbon structures is shown in [Wackerfuss, 2009], where the molecular model is plugged within the formalism of the finite element method.

Virtual experiments on representative lattice

We show subsequently the biaxial tensile tests performed on the graphene lattice sample using the presented molecular model implemented in the inhouse MATLAB code. The scheme of the representative lattice element with the symmetry boundary conditions is illustrated in Fig. 1 A zoom on the bulk atom i is shown with the bond length and angle given in the undeformed (r 0 , θ 0 ) and deformed (r ij , θ ikl ) configurations.

the boundary of the graphene sheet is composed of lines L 1 to L 4 . Atoms which are on the lower and left lines L 1 and L 3 of the sample are pinned with u 2 = 0 and u 1 = 0, respectively. The boundary atoms which belong to the upper and right lines, L 2 and L 4 , have a given displacement u 2 = ū1 and u 2 = ū2 , respectively, to produce the stretch λ. Overlay plot of the undeformed and deformed lattice sample is shown in Fig. 1 b) for the case where λ1 = 1 and λ2 = 1.15 (deformed shape is given with scale factor, and really small lattice sample is shown due to visibility). Potential parameters related to the initial, undeformed geometry of the lattice (r 0 , θ 0 ) are shown in the zoom on the bulk atom i. Given deformation ( λ1 = 1, λ2 = 1.15) maps the atom i and his neighborhood to deformed configuration (r 0 → r ij , θ 0 → θ ikl ), as schematically shown on the right of Fig. 1 b).

Equivalent continuum

In order to construct equivalent continuum potential W fit (λ 1 , λ 2 ) we determine the equilibrium potential energy of atomistic system. Note that the boundary atoms have two neighboring atoms instead of three for the bulk ones (on the right side of Fig. 1 b) bulk atom i together with it's neighbors j, k, l is shown). Thus, the equilibrium energy density of the finite size lattice sample depends on the ratio of boundary and bulk atoms. The smaller the ratio the less is the number of the boundary with respect to bulk atoms, and the energy density converges to the infinite lattice (where all atoms have full neighborhood). This convergence tendency is visible, e.g., on Fig. 10 b) in [START_REF] Marenić | Homogenized elastic properties of graphene for small deformations[END_REF]. Our approach which considers the development of the surrogate continuum model based on the virtual test on representative graphene samples works well for both really small and/or narrow lattice (e.g. graphene nano-ribbons [Xu, 2009]) as well as for the larger ones. We compute the equilibrium potential energy of atomistic system for the series of loading cases. These loading cases are designed to form the uniform grid in the space of λ 1 , λ 2 in the range

λ 1 = {1, . . . , λ}, λ 2 = {1, . . . , λ}, (30)
resulting with the cloud of points Ŵ (λ 1 , λ 2 ), shown as dots in Figure 2. Note that in the above equation the given values of stretch λ ≥ 1, which corresponds only to in-plane tension1 . The energy distribution for the series of loading cases is further used to perform a least squares, polynomial surface fitting (see Figure 2) with SED potential given as

W fit (λ 1 , λ 2 ) = n i n j a ij λ i 1 λ j 2 , ( 31
)
where i and j are the degree in λ 1 and in λ 2 , respectively. The total degree of the polynomial is the maximum of i and j. Note that the total degree of the polynomial cannot exceed the maximum of i and j. Hence if i = j = 5, for instance, the coefficients a ij = 0 if i + j > 5. As mentioned above for the initial iteration we have x

i = X i , that is the r = r 0 and θ = θ 0 , and considering our choice of the modified Morse potential and equations ( 28) and ( 29), the value of SED turns to be zero for the unit stretch. No other constraints were introduced in the model fit.

Using fitted, polynomial potential (31), and following the equations ( 16) to ( 18) we can calculate the stress components. We present in Figure 3 surface plots (closed form, polynomials obtained by energy fit) of the nonzero stress components. Note that stress components are conveniently written in terms of the stress vector, as shown in appendix A. In addition, using spectral decomposition (17) one computes directly principal stresses using (18). In our consideration, this further simplifies stress vector to two components. Knowing principal stresses and directions, finite element implementation further relays to the internal force vector, and tangent elasticity matrix computation. Stress vector (in terms of principal stresses) is directly plugged into membrane finite element equations, for details see appendix, equations (A.5) and (A.6). From the plots in Figure 3 we can note that the maximum stress does not correspond to the maximum biaxial

1 1.05 1.1 1 1.05 1.1 -0.1 0 0.1 0.2 0.3 0.4 λ 1 W , eV/ Å2 λ 2 Figure 2:
The surface plot of the strain energy density polynomial surface fit W fit (λ 1 , λ 2 ). The fit is obtained by means of series of biaxial tests performed by molecular mechanics simulation ( Ŵ (λ 1 , λ 2 )) represented with circular markers on the grid 0 ≤ λi ≤ 1.15, i = 1, 2.

stretches, that is, to the case when λ i = λ. This phenomena is related to the given hexagonal lattice structure, chosen interatomic potential, as well as the finite deformation regime and will be further analyzed and explained in sequel.

Analogously to the stress, using fitted potential and following (20) we present the surface plots of the components of material part of the tangent elasticity tensor using its reduced form D ij , on the Figure 4. We note that D 12 = D 21 due to symmetry.

Analysis and verification of the constitutive model

We seek first to explain the effect illustrated in Figure 3, i.e., to answer why the stress is not maximum for the maximum biaxial stretch. Thus, we present a variation of the stress component S 11 with λ 1 , taking the prestretch in the perpendicular direction (λ 2 ) as a parameter, see left plot on the Figure 5. We can observe that for larger deformation (roughly λ 1 > 10%), value of the S 11 decreases for higher pre-stretch (λ 2 ). We can also note this effect by plotting S 11 vs. λ 2 with the parameter λ 1 . The latter is depicted in the right plot in the Figure 5, where we clearly see the decrease of S 11 with evolving stretch in direction 2, for higher pre-stretch in direction 1. Needless to say, this effect is fully captured in the developed homogenised, continuum model, i.e., the stress decrease is noticeable in terms of the tangent elastic modulus. Namely, the component D 12 of the tange nt elastic modulus (shown on the bottom plot in Figure 4) which governs the relation between S 11 and strain in direction 2, becomes negative for large deformation. An analogous effect can be seen for the stress component S 22 , thus these plots are not shown.

In order to further explain this effect of stress decrease, we turn to the study of the nanostructure of graphene. Moreover, we designed an equibiaxial, so-called half snail loading program depicted in Figure 6 a) and we trace the corresponding lattice deformation. Given loading program consider stretching first half of the loading time in direction 2 (up to λ 2 = λ 2,max ) while holding the lattice sample with lambda 1 = 1. Second half of the loading program considers holding the stretch λ 2 = λ 2,max while loading in direction 1. Note that the stretches we did in our virtual experiments did not go more than about 15% in each principal direction. This value is related to the chosen interatomic potential, more precisely to its pair part. Namely the derivative of the pair part (28) with respect to the pair bond separation r gives a pair force F p (r). This function has a peak value after which the force diminishes with further bond separation. There are two reasons why we are not taking this in our consideration. The first one pertains to the choice of interatomic potential which is not suitable for the bond fracture modeling. The second is related to the implemented Newton iterative algorithm which is not suitable when the effect similar to material softening occurs. During the loading program we follow the deformation in C-C bonds by selecting the bulk atom i and his neighbours j and k (see Figure 6 b)), omitting l due to symmetry. We note that the the choice of the atom we trace is free, given that it is far enough from the boundary. The bond separation (∆r) evolving with the loading (given as pseudo-time) is shown on Figure 6 c). We observe that for the first half of the load program with stretch in direction 2, both ∆r ij and ∆r ik are increasing. However, in the second part of the load program the bond separation ∆r ij passes a peak value. This effect is related to the interplay of the pair and angular parts of the used modified Morse potential, as well as to the to the geometric nonlinearity, i.e. the large rotations of the bonds. The latter is causing the influence of the angular term of the modified Morse potential (eq. ( 29)) to overcome the influence of the pair part (eq. ( 28)). Finally, this causes the global response of the nanostructure defined in terms of stress-stretch diagram to show the stress decrease in large deformation regime which is noticed in the developed continuum model. Next, we show the cross-sections from the Figure 4 similarly like it was done for the stresses in Figure 5. This gives the evolution of the components of reduced tangent elastic modulus, namely we show D 11 versus λ 1 with the pre-stretch λ 2 as a parameter. Likewise, D 22 versus λ 2 with the parameter λ 1 is shown in Figure 7. The thick lines with markers (in Figure 7) denote the evolution of D ii without pre-stretch, i.e. dashed line with circular markers represents the D 11 (λ 1 ) for λ 2 = 1. Analogously, full line with square markers shows D 22 (λ 2 ) for λ 1 = 1. Any material description of an elastic constitutive law for large deformations should reduce to Hooke's law for the case of small deformation. Thus the we proceed with the constitutive model verification. We perform this verification by comparison with the limiting case of small deformations considered in [START_REF] Marenić | Homogenized elastic properties of graphene for small deformations[END_REF]. The numerical values of the initial stiffness obtained from the fitted continuum model are as follows

D 11 | λ i =1 = 956.95GPa, D 22 | λ i =1 = 876.35GPa. (32)
We note that this corresponds to the results presented in [START_REF] Marenić | Homogenized elastic properties of graphene for small deformations[END_REF] considering the case E| 'V' . Moreover, plotting the evolution of the components of reduced tangent elastic modulus with the increasing (corresponding) pre-stretch we end up with a band which is depicted for the whole range λ i = 1 . . . λ in

Conclusion and perspective

We have shown the continuum formulation of the in plane behaviour of graphene based on the hyperelastic potential given in terms of the principal stretches. The developed model is based on the well established continuum mechanics framework, and fully capable to reproduce the linear elastic behaviour in small strain regime as well as the stress release caused by intrinsic geometric non-linearity of the interatomic bonds that occurs in large strain regime. Moreover, this formulation fits perfectly to the finite element implementation of the elastic membranes based on the Cartesian structure, i.e., co-axial energy-conjugate pairs of stress and strain.

The procedure performed here concerns virtual experiments performed on the graphene samples. Thus, this approach should work for lattices of other two-dimensional materials like e.g. boron nitride [Novoselov et al., 2005,Topsakal andCiraci, 2010] which may have more complex lattice. The latter yields at the bottom line more complicated deformation mechanism on the lattice level and may preclude the Cauchy-Born rule, as a common link between atomistic and continuum scales, to be valid. Thus we plan to confront our large strain surrogate continuum model based on the numerical homogenization procedure with the Cauchy-Born based approach. In addition, we presented here only the modeling of the in plane large deformation of single layer graphene sheet, however, the extension to axisymmetric (like the Carbon Nano Tube (CNT) ) or arbitrary curved membrane is possible (see e.g. [Ibrahimbegović, 1994]), and will be concerned in future research.

A Finite element implementation

The details about the 2D plane elastic membrane finite element can be found in most of the books dealing with nonlinear solid mechanics e.g. [START_REF] Ibrahimbegovic | Nonlinear Solid Mechanics[END_REF],Wriggers, 2008,Belytschko et al., 2000] and will not be discussed in detail. We rather illustrate the main steps needed in the finite element approximation concerning, i.e. we recast in matrix form the results obtained for the constitutive law in terms of principal stretches. First, we define the coordinate representation of the principal vectors in the two-dimensional setting under consideration as andb) the bulk atom i with its neighbourhood. The given load program causes the lattice deformation as presented in c). Due to symmetry, the bond separation ∆r il is equal as ∆r ij , and is thus omitted.

In the above equation the matrix T is created by ordering the tensor product of eigenvectors (A.1) in vector notation where a = 1, . . . , n nodes denotes the node, i = 1, 2 the degree of freedom (n dof ), index p = n dof (a -1) + i, and matrix B stands for the derivatives of the shape functions.

n 1 ⊗ n 1 → n 1 n T 1 =
It is also convenient to write tangent elasticity tensor, given in (19), that connects stress and strain through C = ∂S ∂E in a matrix form, using the transformation matrix T from (A.3), and D ij from (20)

C → C (3×3) = TD ij T T + s 1 -s 2 λ 2 1 -λ 2 2 gg T . (A.6)
In the above equation auxiliary term g T = [-sin 2α sin 2α cos 2α] is used to express the geometric part of the tangent elasticity tensor in more compact form. We note further, that the element tangent stiffness matrix K e also consists of a material and of a geometric part, and is for the 2D plane elastic membrane finite element given as where a, b = 1, . . . , n nodes denote the node, i, j = 1, 2 the degree of freedom, and indexes p, and q are p = n dof (a -1) + i and q = n dof (b -1) + j. where N a represents element shape function.

Figure 1 :

 1 Figure 1: Scheme of the lattice sample with symmetry BCs used for biaxial tensile tests a). The envelope of the sample is composed of lines L 1 . . . L 4 which coincides with boundary atoms. The overlay plot of the undeformed (red) and deformed (grey) shapes of the small graphene lattice sample for the case λ1 = 1, λ2 = 1.15 b).A zoom on the bulk atom i is shown with the bond length and angle given in the undeformed (r 0 , θ 0 ) and deformed (r ij , θ ikl ) configurations.

Figure 3 :

 3 Figure 3: Surface plot of the nonzero stress components vs. principal stretches.

Figure 4 :

 4 Figure 4: Surface plot of the components of reduced tangent elastic modulus. The values are given in GPa.

Figure 7 .

 7 Figure 7. Note that the increase of the pre-stretch causes the decrease of stiffness. Acknowledgements: This work was supported by Labex MS2T post-doctoral fellowship.

Figure 5 :

 5 Figure5: The stress-stretch plots showing the component S 11 versus: stretch λ 1 with parameter λ 2 (left plot), and stretch λ 2 with parameter λ 1 (right plot). The parameter is in the range λ i = 1, . . . λ, where the lower-most stress curve corresponds to the value of the paremeter λ i = 1, while the upper-most corresponds to λ i = λ.

  α denotes the angle between the first principal direction and axis x 1 . Using the component form of the (9) the value of α is choose to order the second Piola-Kirchhoff stress and Green-Lagrange strain tensor components in a vector as S → s T = [S 11 , S 22 , S 12 ], E → e T = [E 11 , E 22 , 2E 12 ], respectively (so that their inner product is preserved). The latter enables to recast the stress spectral decomposition from equation (17sin 2 α sin 2 α cos 2 α sin α cos α -sin α cos α

Figure 6 :

 6 Figure 6: Load program showing the given displacements (upper and right edges) with respect to the pseudo time a),and b) the bulk atom i with its neighbourhood. The given load program causes the lattice deformation as presented in c). Due to symmetry, the bond separation ∆r il is equal as ∆r ij , and is thus omitted.

  as the columns in T (analogously for the n 2 ⊗ n 2 ). Note that the last result for the stress vector (s) can further be directly used for the calculation of the internal force vector of the finite element f

Figure 7 :

 7 Figure 7: Evolution of the diagonal components of the reduced tangent elastic modulus with the associated stretches.The thick lines with markers denote the evolution of D ii without pre-stretch.

  The geometric part of the tangent stiffness appears only in large displacement problems, and it depends directly on current stress values through

						D11, param. λ2; D22, param. λ1
			1000
										D11, λ1
			900						D22, λ2
			800
			700
		Dii, GPa	400 500 600
			300
			200
			100
				0
				1		1.05	1.1			1.15
								λi
	G ab = S 11	∂N a ∂x	∂N b ∂x	+ S 22	∂N a ∂y	∂N b ∂y	+ S 12	∂N a ∂x	∂N b ∂y	+	∂N a ∂y	∂N b ∂x	,	(A.8)

The compressive stresses even the ones transmitted by the substrate causes out of plane buckling of the SLGS , see[Zhang and Arroyo,

2013] for the analysis of the interplay between localized folds and distributed wrinkling of graphene deposited on planar surfaces.


Eduard Mareni
email: eduard.marenic@utc.fr

Delphine Brancherie
email: delphine.brancherie@utc.fr

Marc Bonnet
email: mbonnet@ensta.fr

Multiscale asymptotic-based modeling of local material inhomogeneities

This research is a first step towards designing a numerical strategy capable of assessing the nocivity of a small defect in terms of its size and position in the structure with low computational cost, using only a mesh of the defect-free reference structure. The proposed strategy aims at taking into account the modification induced by the presence of a small inhomogeneity through displacement field correction using an asymptotic analysis. Such an approach would allow to assess the criticality of defects by introducing trial micro-defects with varying positions, sizes and mechanical properties.

Introduction

The role played by defects in the initiation and development of rupture is crucial and has to be taken into account in order to realistically describe the behavior till complete failure. The difficulties in that context revolve around (i) the fact that the defect length scale is much smaller than the structure length scale, and (ii) the random nature of their position and size. Even in a purely deterministic approach, taking those defects into consideration by standard models imposes to resort to geometrical discretisations at the defect scale, leading to very costly computations and hindering parametric studies in terms of defect location and characteristics.

Our current goal is to design an efficient two-scale numerical strategy which can accurately predict the perturbation in terms of displacement concentration caused by an inhomogeneity in elastic (background) material given by elasticity tensor (C). To make it computationally efficient, the analysis uses only a mesh for the defect-free structure, i.e. the mesh size does not depend on the (small) defect scale. The latter is instead taken into account by means of a multiscale asymptotic expansion (see e.g. [START_REF] Brancherie | Effect of surface defects on structure failure: a two-scale approach[END_REF][START_REF] Dambrine | A multiscale correction method for local singular perturbations of the boundary[END_REF]), in which the concept of elastic moment tensor (EMT) [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF][START_REF] Ammari | Polarization and Moment Tensors: with Applications to Inverse Problems and Effective Medium Theory[END_REF] plays an important role.

Problem definition

The main feature of the proposed approach is as briefly stated in the introduction to accurately predict the perturbation caused by an inhomogeneity in elastic (background) material, without explicitly modeling it. That is, we consider an elastic body occupying a smooth bounded domain Ω ⊂ R 3 characterized by elasticity tensor C, as depicted in Figure 1 a), and we seek the solution of the background problem, given below. Having this solution in hand, we proceed with the asymptotic expansion of the displacement perturbation. Adding these two elements together should correspond to what we call the transmission (perturbed) problem, without really solving it. The transmission problem for a small trial inhomogeneity involves a small elastic (C ) inhomogeneity located at z ∈ Ω, of characteristic linear size a, occupying the domain B a := z + aB, where B ⊂ R 3 is a smooth fixed domain centered at the origin and defines the inhomogeneity shape, see Figure 1 b). The elastic properties of the whole solid are defined as

C a := Cχ Ω\Ba + C χ Ba , (1)
where χ D is the characteristic function of a domain D. For the later reference we also introduce the material contrast, or elastic tensor perturbation, as ∆C = C -C.

Background solution. The background solution in terms of displacement field u arising in the reference solid Ω due to prescribed excitation (f , t, ū), is defined by the problem div(C :

ε[u]) + f = 0 in Ω, t[u] = t on Γ N , u = ū on Γ D , (2)
where Γ N and Γ D typically denote Neumann and Dirichlet parts of boundary, while ε[w] and t[w] are linearized strain tensor and traction vector associated with a given displacement w, which are defined as

(a) ε[u] = ∇ sym u, (b) t[u] = (C : ε[u]) • n, (3)
with n being the unit outward normal to Γ N .

Transmission solution

The displacement field u a arising in the solid containing the small inhomogeneity B a due to the same prescribed excitation (f , t, ū), solves the transmission problem with elastic properties distributed as in (1)

(a) div(C a : ε[u a ]) + f = 0 in Ω, (b) t[u a ] = t on Γ N , (c) u = ūa on Γ D , (4)
where the perfect-bonding transmission conditions are implicitly enforced via (4a) written in the distributional sense and the usual smoothness assumption u a ∈ H 1 (Ω).

Multiscale approach by asymptotic expansion. Let us first define the displacement perturbation

v a := u a -u, (5)
as the difference of the total and unperturbed displacement corresponding to problems (4) and (2), respectively. Two distinct asymptotic expansions of v a arise, namely a near-field and a far-field expansion. This distinction naturally involves two scales: the "scale a" corresponding to characteristic length a of the inhomogeneity, and the "scale 1" on the size of structure. This multiscale character of the asymptotic expansions is schematically depicted in the Figure 2. Furthermore, this description requests two variables [START_REF] Dambrine | A multiscale correction method for local singular perturbations of the boundary[END_REF][START_REF] Brancherie | Effect of surface defects on structure failure: a two-scale approach[END_REF]: • 'x' representing coordinates at the structure scale (see Fig. 2 a)), and

a) b)
• 'x := (x -z)/a' representing coordinates at the defect scale (see Fig. 2 b)).

The near-field expansion is given by

v a (x) = av B (x) + o(a), (6)
where v B is the solution in terms of displacement perturbation of the auxiliary problem of a perfectly bonded inhomogeneity B embedded in an infinite elastic medium Ω = R 3 and subjected to the constant remote stress C : ∇u(z). Such solutions v B are known analytically for simple inhomogeneity shapes (see [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]) in which case they correspond to the famous Eshelby inclusion problem [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion and related problems[END_REF]. Expansion ( 6) is valid for finite x, i.e. within a neighbourhood of B a whose linear size is O(a). More precisely, the solution (6) represents a zoom on the ihomogeneity and thus it exists on the "scale a" (see Fig. 2 a)).

The far-field expansion, on the other hand, is given by

v a (x) = -∇ 1 G(z, x) : A(B, C, ∆C) : ∇u(z)a 3 + o(a 3 ), x = z, (7)
and is valid at a finite (independent on a) distance from the inhomogeneity, that is, on the scale of the structure. In this work, as a first step, we address the computation of the far-field approximation [START_REF] Marenić | Asymptotic analysis for the multiscale modeling of defects in mechanical structures[END_REF] for given inclusion location z (thus purely deterministic approach) and varying evaluation point x. The key "ingredients" of the far-field expansion [START_REF] Marenić | Asymptotic analysis for the multiscale modeling of defects in mechanical structures[END_REF] in terms of displacement perturbation are: the gradient of the background solution on the inhomogeneity site ∇u(z), Green's tensor for the domain Ω denoted as G(z, x) (i.e. it's gradient with respect to the first argument ∇ 1 G(z, x)), and the elastic moment tensor (EMT) A(B, C, ∆C). We proceed with the brief introduction, definition and properties of these quantities.

Elastostatic Green's tensor. The elastostatic Green's tensor G(ξ, x) is defined as a solution of

div(C : ε[G(•, x)]) + δ(• -x)I = 0 in Ω, t[G(•, x)] = 0 on Γ N , G(•, x) = 0 on Γ D . (8)
Defined this way Green's tensor gathers the three linearly independent elastostatic displacement fields G k (•, x) resulting from unit point forces δ(• -x)e k applied at x ∈ Ω along direction k. We further introduce a decomposition of G as

G(•, x) = G ∞ (• -x) + G c (•, x), (9)
where G ∞ is the singular, infinite-space Green's tensor, such that

div(C : ε[G ∞ ]) + δI = 0 in R 3 , |G ∞ (ξ -x)| → 0 as |ξ -x| → ∞, (10)
and the complementary Green's tensor G c is bounded at ξ = x and represents a correction of G ∞ due to the finite size of Ω. More precisely, by a superposition argument, G c solves the elastostatic boundary-value problem (BVP) with regular boundary data and zero body force density:

div(C : ε[∇G c ]) = 0, ∇G c = -∇G ∞ on Γ D , t[∇G c ] = -t[∇G ∞ ] on Γ N , (11)
where the boundary data is related to the trace of ∇G ∞ on ∂Ω, treated columnwise as displacement fields, and the traction vectors associated with those displacements. Note that problem (11) is posed directly i terms of the gradient of complementary Green's tensor, mainly to avoid the loss of accuracy caused by numerical derivative (see [START_REF] Marenić | Asymptotic analysis for the multiscale modeling of defects in mechanical structures[END_REF] for details). Note also that problem (11) in fact involves 9 loading cases for 3D conditions (and 4 for 2D conditions).

Elastic moment tensor (EMT).

The EMT [START_REF] Ammari | Polarization and Moment Tensors: with Applications to Inverse Problems and Effective Medium Theory[END_REF] is the key mathematical concept in the asymptotic expansion of v a . It carries important microstructural information, namely the material contrast ∆C, the inhomogeneity shape B and its orientation (see e.g. Definition 2.3. in [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF] for details and properties of the EMT). We focus in this work on the EMT A associated with an ellipsoidal inhomogeneity (B, C + ∆C) embedded in a medium with elasticity tensor C, which is given by [3]

A = |B|C : (C + ∆C : S) -1 : ∆C, (12)
where S = S(B, C) denotes the (fourth-order) Eshelby tensor of the inclusion B [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]. We will be dealing in this paper with the special case of plane strain, corresponding to an ellipsoidal inclusion B infinitely elongated in the x 3 direction, assuming C and C to be both isotropic. In this case, the nonzero components of the Eshelby tensor S are given explicitly as

S 1111 = A(1 -m)(3 + γ + m), S 1122 = A(1 -m)(1 -γ -m), S 2222 = A(1 -m)(3 + γ -m), S 2211 = A(1 + m)(1 -γ + m), (13)
S 1212 = A(1 + m 2 + γ),
where

A = 1/[8(1 -ν)], γ = 2(1 -2ν)
, and m = (a 1 -a 2 )/(a 1 + a 2 ), with a 1 and a 2 denoting the semi-axes of the ellipse B in principal directions x 1 and x 2 , respectively (while a 3 → ∞).

Numerical example

The capabilities of the proposed strategy are illustrated on the academic example involving the biaxial load of the square-shaped membrane as shown on Fig. 3. This pre- liminary example, concerns a plane strain case involving a single, local, elliptic inhomogeneity, and will be used to further clarify the procedure and demonstrate the performance of the proposed strategy. The example of finite element meshes of unperturbed (M 0 ) and reference (M r ) models are shown on the Figure 4. Proposed approach considers first the FE computation of the background solution u defined by problem (2) using a coarse mesh, M 0 shown on Fig. 4a). In this computation the inhomogeneity is not taken into account. Following steps of the proposed strategy aim to correct u by adding far field correction [START_REF] Marenić | Asymptotic analysis for the multiscale modeling of defects in mechanical structures[END_REF]. Thus, for a given inhomogeneity position (z ∈ Ω) we first extract the gradient of the background displacement ∇u(z). Next, we turn to the computation of the Green's tensor considering a decomposition (9). The infinite part is in our case of isotropic background material a simple derivative of the textbook fundamental (Kelvin) solution, given in (4.57) in [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF]. The same infinite part serves as a BC for complementary part (11), thus a special attention needs to be given to the preparation and correct imposition of these BC's (see [START_REF] Marenić | Asymptotic analysis for the multiscale modeling of defects in mechanical structures[END_REF]). With the properly defined BC's, the computation of the complementary part ( 11) is performed numerically using finite element method. This step also considers a finite element computations without any (fine) meshing at the defect scale. Finally, for the given shape and material properties (B, C ) of the inhomogeneity we compute the EMT A(B, C, ∆C) using (12). Note, that for the case of isotropic background and inhomogeneity this can be solved analticaly, that is without any additional computation.

Having these results in hand, all the ingredients required for evaluation of the far-field asymptotic approximation [START_REF] Marenić | Asymptotic analysis for the multiscale modeling of defects in mechanical structures[END_REF] of v a are available, and this solution may be post-processed for evaluating the defect criticality. Obtained results are to be compared with reference solution performed on a really fine mesh M r given in th Fig. 4b). Reference model explicitly models the inhomogeneity, and is considered as a solution of the transmission problem (4), denoted as u num a (or in terms of perturbation v num a ). Note that the proposed approach avoids any computation on the fine mesh, both u and ∇G c are computed on M 0 (other steps pertain merely on the evaluation of analytical expressions). Approach based on asymptotic expansion takes, thus, a fraction (about 11%) of CPU time needed for the reference solution. We will discuss further the accuracy and convergence of the proposed approach.

Accuracy. The comparison of our asymptotic approximation v a and referential numerical solution v num a are given first along the characteristic crossection A -A, shown on the scheme in Fig. 3. The inhomogeneity is in this case taken simply as circular hole, and a perfect agreement is achieved, solution obtained by asymptotic expansion coincides with the referential solution on the evaluation points x 0 (see Figure 5 (left)). Moreover, to verify further proposed approach we proceed with computing the far-field asymptotic approximation (7) also on fine mesh M r , i.e., for each node x r of the fine mesh. As mentioned above, background solution and complementary part of the Green's tensor are computed on the coarse mesh M 0 , and we can using the FE interpolation ϕ write

v a (x r ) = -(∇G ∞ (x = x r ) + ϕ∇G c (x 0 )) : A : ∇u(z)a 3 , x r , x 0 = z. (14)
Let the difference of the two fields on the fine mesh M r be given as v e = v num a -v a . Then, the comparison can be given as contour plot in Figure 5 (right). This plot shows a good agreement over the whole domain, with the negligible error appearing only in the close vicinity of the inhomogeneity. Convergence. The proposed approach shows a good performance both in terms of efficiency and accuracy. We proceed here with the convergence test again by comparison of the asymptotic approximation with the reference results, but here for the varying characteristic inhomogeneity size a. The reference results for various inhomoeneity types and varying sizes a were made using really fine FE model (like on Figure 4 (right)). While each of this computations inherits the mentioned problems of the severe mesh refinement and cost, we can note in this convergence study the key advantage of the proposed approach based on asymptotic expansion [START_REF] Marenić | Asymptotic analysis for the multiscale modeling of defects in mechanical structures[END_REF] and the multiplicative decomposition therein.

The decomposition makes possible that the information about the (virtual) inhomogeneity is stored in EMT independently from Green's tensor and background solution given by ∇u. Thus, asymptotic approximation for varying inhomogeneity sizes, types (here shown only ellipsoidal) and material contrasts, boils down to simple recomputing of the EMT following (12). The computation of EMT for ellipsoidal inhomogeneity (12) mainly hinges on the evaluation of Eshelby tensor S, which for the isotropic background material covered within this work concerns merely analytic expressions. While already for a single run our approach was taking just a fraction of the computational cost of the fine FE model, for multiple runs there is practically no added cost! This way we can achieve a significant speed-up in predicting the perturbation caused by varying inhomogeneity sizes, types and material contrasts, which is of great practical interest.

For the convergence analysis we define discrepancy of reference and asymptotic solutions as

R(a) = ||v -v num a || a 2 = ||v e || a 2 = o(1). (15)
In order to find a O(a 2 ) remainder in (15), which is the outer asymptotic behavior of the solution, we evaluated R(a) using a L 2 norm

||v e || L 2 = Ω\D v e v T e dΩ, (16)
over a domain Ω\D, where D is a fixed neighbourhood around inhomogeneity. In this case we have choosen a disk with a radius R d centered at the inhomogeneity (as shown for one a size in Figure 4b), with R d kept fixed as a varies. The plot of R(a) versus the normalized characteristic size of the inhomohenity a/L is depicted in Figure 6 for the varying inhomogeneities. We show the results for soft (E /E = 0.5), hard (E /E = 1.5), and nearly rigid (E /E = 10) inhomogeneities of circular (a 1 /a 2 = 1) and ellipsoidal (a 1 /a 2 = 2) shape. A desired convergence trend (O(a 2 )) is clearly visible for all considered cases.

Conclusion and perspectives

A numerical strategy for predicting the perturbation caused by an inhomogeneity in an elastic solid is outlined. We focus on the evaluation of the far field displacement correction. Thanks to the use of an asymptotic expansion, we achieve the major advantage which pertains to the fact that the defect scale is not meshed. The computations involved are straightforward and include the computation of background solution, followed by the Green's and elastic moment tensors. A very good performance in terms of accuracy and convergence is presented on the academic example.

The next steps of this work include (i) matching the near and far-field asymptotic expansions to obtain uniform expansions, and (ii) applying the developed strategy to assessing the criticality of defects by considering virtual micro-defects and varying their positions, sizes and mechanical properties.

Figure 1 :

 1 Figure 1: Model scheme showing the geometric setting, that is, unperturbed a) and perturbed b) domains, Ω and Ω a , respectively. The inhomogeneity B a located at z is the shaded subdomain in (b).

Figure 2 :

 2 Figure 2: Scheme showing the far-field a) and near-field b) expansions. Far-field is represented with coordinates at the structure scale, while the near-field represents zoom on the ihomogenity represented by coordinates at the structure scale.

Figure 3 :

 3 Figure 3: Scheme of the geometry and boundary conditions of the biaxially loaded square membrane.

Figure 4 :

 4 Figure 4: Finite element mesh of unperturbed i.e. background a) and reference models b). The unperturbed mesh consists of 289, and reference of 2464 nodes.

Figure 5 :

 5 Figure 5: The comparison of the asymptotic approximation v a and reference solution v num a , along crossection A -A (left), and in terms of contour plot of perturbation error v e (right). Both are given for a membrane with a circular hole on site z of radius a = 2 mm.

1 /a 2 = 1 ,Figure 6 :

 1216 Figure 6: The plot of discrepancy versus the normalized inhomogeneity size. The soft and hard inclusions with elliptic and circular shape are considered.


Eduard Marenić
email: eduard.marenic@utc.fr

Adnan Ibrahimbegović
email: adnan.ibrahimbegovic@utc.fr

on Multi-scale Computational Methods for Solids and Fluids Multiscale elastic properties of graphene for moderate deformations

Keywords: graphene, virtual tests, elastic properties, molecular mechanics

A multiscale approach to simulate graphene combining molecular mechanics and continuum mechanics is presented. A series of virtual expermiments is performed first, to develop computationally cheaper, substitute continuum model of the discrete molecular model of pristine graphene. Continuum and molecular models are then coupled providing the means to simulate a defected graphene while keeping the computational cost.

Introduction

Nanoscale systems and processes based on graphene [START_REF] Singh | Graphene based materials: Past, present and future[END_REF] are becoming more viable for engineering applications, however, our ability to model their performance remains limited. The main challenge is that parts of graphene-based devices modelled with discreet models (such as molecular mechanics (MM) and molecular dynamics (MD) ) typically contain extremely large number of particles, even though the actual physical dimension may be quite small. We focus in this work upon the MM neglecting both the dynamic effects and the thermal effects, used for quasi-static loading applications with the assumption of the zero Kelvin temperature. The equilibrium configuration of graphene corresponds to a state of minimum energy of the particle system. It is assumed here that the initial configuration is at equilibrium. Moreover, for the MM model we use classical interatomic potential depending on the position of atomic nucleus, namely modified Morse potential (see [START_REF] Belytschko | Atomistic simulations of nanotube fracture[END_REF][START_REF] Marenić | Homogenized elastic properties of graphene for small deformations[END_REF]). Regardless of the above assumptions, following explicitly the trajectory for large number of degrees of freedom easily becomes intractable. Thus we reach for: 1. A substitute, continuum model which models the average mechanical behaviour of atomic system (in Sec. 2), and 2. A coupled multiscale (MS) model consisting of developed continuum and molecular model (in Sec. 3).

Continuum modeling of graphene

We tend here to give a combined MM and continuum approach for the study of the in plane mechanical behaviour of single layer graphene sheet under moderate deformations. We note that continuum model of graphene is extremely simplified, however, experimental investigation [START_REF] Lee | Measurement of the elastic properties and intrinsic strength of monolayer graphene[END_REF] show that graphene undergos large deformations, with highly nonlinear behavior and still remain elastic, with stable bonds and intact bond topology. Thus continuum modeling brings insight and quantitative information on the relevant physical phenomena, and seams optimal for the current and potential future applications. This further allows us to adapt and use a finite strain elastic model for the large range of phenomena and motivates us to benefit from somewhat simplified nonlinear membrane theory. We seek to adjust the nonlinear membrane theory which includes, as a special case, the hyperelastic model in terms of strain energy density (SED), W , as a function of principal stretches W (λ 1 , λ 2 ). This is an elegant alternative for the construction of the elastic constitutive response that satisfies the material indifference and isotropy restrictions (often used to characterise rubberlike materials [START_REF] Ibrahimbegović | A consistent finite element formulation of nonlinear membrane shell theory with particular reference to elastic rubberlike material[END_REF]). In order to construct SED potential W fit (λ 1 , λ 2 ), we determine the equilibrium potential energy of atomistic system for the series of biaxial loading cases performed on lattice sample, see Fig. 1 (left). The loading cases are designed to form the uniform grid in the space of λ 1 , λ 2 resulting with the cloud of points shown as dots in Fig. 1

(right).

These results are further used to perform a surface fit given as W fit (λ 1 , λ 2 ). With this result in hand, we

Zigzag (Z) c)
S := 2 ∑ i=1 1 λ i ∂W (λ 1 , λ 2 ) ∂λ i s i n i ⊗ n i , C := 2 2 ∑ i=1 ∂s i ∂C n i ⊗ n i + 2 2 ∑ i=1 s i ∂ ∂C (n i ⊗ n i ). (1)
Note that the matrix representations of these results can further be directly used for the calculation of the internal force vector and the element tangent stiffness matrix of the standard 2D large displacements elastic membrane finite element. Having the continuum description of the graphene behavior we turn to the coupling of molecular and continuum.

Multiscale atomistic-to-continuum coupling

In this paper coupling is based on the bridging domain (BD) method [START_REF] Belytschko | Coupling methods for continuum model with molecular model[END_REF]. The basic idea is to divide the computational domain Ω in three subdomains, atomistic Ω a , continuum Ω c and their overlap

Ω b = Ω a ∩ Ω c
, where the compatibility is enforced by Lagrange multiplier (LM) method. The total potential energy of the system can be written as

E tot,w = E a tot,w (d) + E c tot,w (u), (2)
where d and u are displacement vectors in the atomistic and continuum domains, respectively, and w in sequel denotes the energy weighting in Ω b . The weighted atomistic and continuum energies are defined as

E a tot,w = ∑ i ∑ j =i w a i j V p + ∑ j =k =i w a i V θ -∑ i∈Ω a w a i fi • d i , (3)
E c tot,w = Ω c w c (X)W dΩ c - Ω c w c (X)u • bdΩ c - Γ c σ w c (X)u • tdΓ c , (4)
where V p and V θ comes from interatomic potential (see [START_REF] Marenić | Homogenized elastic properties of graphene for small deformations[END_REF]), and t, b is related to continuum tractions and volume forces, while fi denotes point forces on atom i. In order to enforce the compatibility between the atomistic and continuum domains, the coupling term, in terms of energy, is added to total energy forming

W L := E tot,w +C = E tot,w + Ω b α 1 λ • (u -d b ) + α 2 ∇λ(∇u -∇d b )dΩ. ( 5
)
where the choice of the weighting parameters α 1 and α 2 determines the coupling by mixing the displacement and strain coupling terms, d b (X) is the interpolated atomistic displacement field in Ω b , and λ denotes LM field.

We give next the weak form of the coupling problem. Let V a and V a 0 be the sets of trial (d) and test (w) functions in Ω a , respectively. Analogously, for the continuum domain let the standard space of admissible solutions be V and the space of test functions V 0 . Furthermore, the space of LM is denoted as M = {λ, µ ∈ H 1 (Ω)}. We further proceed to the minimising of the functional in (5) which leads to the saddle point problem, which can be written in terms of its weak form:

Find (u, d, λ) ∈ V × V a × M such that G c w (u; v) + G a w (d; w) + (λ, v -w b i| i∈Ω b ) C = 0 ∀(v, w) ∈ V 0 × V a 0 , (µ, u -d b ) C = 0 ∀µ ∈ M , (6)
where w b denotes interpolated test atom displacement field, and the terms defining the weak form of equilibrium with the scaling in the overlap are as follows

G c w (u; v) := Ω c w c ∇ s v • σ(∇ s u)dΩ - Ω c w c v • bdΩ - Γ c σ w c v • tdΓ, (7)
G a w (d; w) := ∑ i∈Ω a ∑ j =i w a i j ∂V p ∂d i • w i + ∑ j =k =i w a i ∂V θ ∂d i • w i -∑ i∈Ω a w a i fi • w i . (8)
The numerical implementation of the given coupling formulation pertains to the choice of the standard finite element and LM field approximation with corresponding shape functions (for the details see [START_REF] Marenić | Atomistic-to-continuum modeling in solid mechanics[END_REF]). We give next one numerical example of graphene sheet with initial crack-like defect, where the defect is modeled simply by removing a line of bonds parallel with the vertical direction. The coupled model will be compared with the fully molecular, see both models in Fig. 2. Moreover, we present the energy error convergence, as a basis for adaptivity procedure, with the increase of the Ω a , for the mode-I type loading of the model from Fig. 2. Thus we define the global the global measure of the energy error as

e E = -ref Ω a \Ω b ref Ω , = 1 n a n a ∑ i i : C : i , (9)
where ref is the strain from the fully molecular model. The description of the models presented in the following results is given in Not surprisingly the convergence is achieved as the size of the molecular domain is increased see Figure 3, where the type of coupling and weighting influences the error (for details see [START_REF] Marenić | Atomistic-to-continuum modeling in solid mechanics[END_REF]). To conclude, we note that for the model denoted as 1 in the Table 1 n do f is reduced by 84%, while the corresponding solution yields negligible error (less then 0.25%) with respect to the fully atomistic model showing really good performance of the modeling strategy.

Figure 1 :

 1 Figure 1: Scheme of the lattice sample with symmetry BCs used for biaxial tensile tests. The envelope of the sample is composed of lines L i which coincides with boundary atoms (left). The surface fit W of SED obtained by series of biaxial tests performed by MM simulation.

Figure 2 :

 2 Figure 2: Graphene sheet with a hypothetical initial crack modeled using the fully molecular model (left) consisting of 10960 atoms and coupled model (right) with the size of atomistic domain 67.4 × 48.7 Å.

Figure 3 :

 3 Figure 3: Convergence of the energy error e E given for different molecular domain dimensions and the different couplings from (5) (L2 (α 1 = 1, α 2 = 0), H1 (α 1 = 1, α 2 = 1)).

Table 1 :

 1 Table1for the three different sizes of the Ω a , and referential model. Model parameters are: number of atoms n a , number of FE nodes n n , number of LM nodes n λ , number of degrees of freedom n do f = 2(n a + n n + n λ ). The data for the models 1, 2 and 3 used in convergence study. The size of the atomistic domain is defined by L 1 × L 3 and given in Å.

	id	ref	1	2	3
	n a	10960	1368	2080	2920
	n n	0	360	343	322
	n λ	0	38	46	54
	n do f	21920	3532	4938	6592
	L 1 163.7134 67.4100 77.0400 86.6700
	L 3 165.4100 48.6500 65.3300 82.0100