Huilong Ren

Xiaoying Zhuang

Nguyen-Thoi Trung

Timon Rabczuk
email: timon.rabczuk@tdtu.edu.vn

A nonlocal operator method for finite deformation higher-order gradient elasticity

Keywords: Nonlocal operator method, finite strain, second/third-gradient strain, invariant, variational principle

We present a general finite deformation higher-order gradient elasticity theory. The governing equations of the higher-order gradient solid along with boundary conditions of various orders are derived from a variational principle using integration by parts on the surface. The objectivity of the energy functional is achieved by carefully selecting the invariants under rigid-body transformation. The third-order gradient solid theory includes more than 10.000 material parameters. However, under certain simplifications, the material parameters can be greatly reduced; down to 3. With this simplified formulation, we develop a nonlocal operator method and apply it to several numerical examples. The numerical analysis shows that the high gradient solid theory exhibits a stiffer response compared to a 'conventional' hyperelastic solid. The numerical tests also demonstrate the capability of the nonlocal operator method in solving higher-order physical problems.

Introduction

Gradient theories have attracted increasing interest due to their capability of describing phenomena such as size effects, edge and skin effects as well as nonlocal effects in materials, which cannot be tackled by conventional continuum mechanics. Gradient elasticity theory introduces an internal length scale and higher-order gradients of the displacement field to account for size effects at the micro-or nano-scale. Gradient theories emerge from considerations of the microstructure in the material at micro-scale, where a mass point after homogenization is not the center of a micro-volume and the rotation of the micro-volume depends on the moment stress/couple stress as well as the Cauchy stress. The starting point of gradient elasticity theory can be traced back to Cosserat theory in 1909 [START_REF] Cosserat | Théorie des corps déformables[END_REF]. A variety of gradient elasticity theories have been proposed which include Mindlin solid theory [START_REF] Mindlin | On first strain-gradient theories in linear elasticity[END_REF][START_REF] R D Mindlin | Micro-structure in linear elasticity[END_REF], nonlocal elasticity [START_REF] A C Eringen | On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves[END_REF], couple stress theory [START_REF] Toupin | Elastic materials with couple-stresses[END_REF][START_REF] Toupin | Theories of elasticity with couple-stress[END_REF][START_REF] Hadjesfandiari | Boundary element formulation for plane problems in couple stress elasticity[END_REF], modified couple stress [START_REF] Yang | Couple stress based strain gradient theory for elasticity[END_REF][START_REF] Tsiatas | A new Kirchhoff plate model based on a modified couple stress theory[END_REF] and second-grade materials [START_REF] Polizzotto | A gradient elasticity theory for second-grade materials and higher order inertia[END_REF].

Gradient elasticity as a generalization of classical elasticity includes the contribution of strain gradients in the strain energy. Different from classical elasticity theory, such consideration enables gradient elasticity to model some interesting phenomena (such as size effect, the stress and strain effects on surface physics, nonlocal effect at micrometer/nanometer scale). Muller and Saul [START_REF] Müller | Elastic effects on surface physics[END_REF] reviewed the importance of surface and interface stress effects on thin films and nano-scaled structures, including the self-organization and elastic driven instabilities of nano-structures. Fischer etal. [START_REF] Fischer | On the role of surface energy and surface stress in phase-transforming nanoparticles[END_REF] studied the role of the surface energy and surface stress in phase-transforming nano-particles. Davydov etal. [START_REF] Davydov | On molecular statics and surface-enhanced continuum modeling of nano-structures[END_REF] showed that a continuum based on gradient elasticity with surface energy contributions can capture size effects that are observed in atomistic simulations. In Refs [START_REF] Artan | Rectangular rigid stamp on a nonlocal elastic half-plane[END_REF][START_REF] Zhou | Investigation of a griffith crack subject to anti-plane shear by using the non-local theory[END_REF][START_REF] Jc Reiher | Finite-element analysis of polyhedra under point and line forces in second-strain gradient elasticity[END_REF] , it is shown that gradient elasticity theory can circumvent stress singularities in local elasticity. Gradient elasticity is closely linked to flexoelectricity, where the strain gradient causes an electromechanical effect [START_REF] Pv Yudin | Fundamentals of flexoelectricity in solids[END_REF][START_REF] Km Hamdia | Sensitivity and uncertainty analysis for flexoelectric nanostructures[END_REF][START_REF] Xy Zhuang | Computational modeling of flexoelectricity-a review[END_REF]. Due to the regularity property of continuum mechanics, gradient elasticity has been applied to problems with strain localization [START_REF] Yang | Micromechanics based second gradient continuum theory for shear band modeling in cohesive granular materials following damage elasticity[END_REF][START_REF] Rodríguez-Ferran | A general framework for softening regularisation based on gradient elasticity[END_REF][START_REF] Placidi | A strain gradient variational approach to damage: a comparison with damage gradient models and numerical results[END_REF]. The micro-structure in continua also plays a crucial role in metamaterials [START_REF] Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF]. Besides strain gradients, velocity gradients enable a more realistic description of dispersive characteristics of the wave propagation in a nonhomogeneous medium such as polymer foams, porous materials, high-toughness ceramics [START_REF] R D Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Papargyri | Static analysis of gradient elastic bars, beams, plates and shells[END_REF] and carbon nanotubes [START_REF] Askes | Gradient elasticity and flexural wave dispersion in carbon nanotubes[END_REF].

Second-gradient elasticity, taking into account the Hessian of the strain tensor, can be viewed as a generalization of gradient elasticity. This theory has been firstly proposed by Mindlin in 1965 [26] in order to account for cohesive force and surface-tension in solids.

Beside the second-gradient term in the displacement field, Polizzotto [START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part i: Constitutive equations and quasi-static behavior[END_REF][START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part ii: Dynamic behavior[END_REF] studied the static/dynamic behavior of linear second-gradient elasticity with second velocity gradient inertia. Askes etal. [START_REF] Askes | A classification of higher-order strain-gradient models-linear analysis[END_REF][START_REF] Askes | Gradient elasticity and flexural wave dispersion in carbon nanotubes[END_REF] showed higher order inertia models are able to realistically describe wave dispersion phenomena in a nonhomogeneous medium. Javili etal. [START_REF] Javili | Geometrically nonlinear higher-gradient elasticity with energetic boundaries[END_REF] derived the governing equations and boundary conditions for third-gradient elasticity with geometrical nonlinearities from variational principles. In the derivation, the bulk and boundary (surface and curve) energies are considered as independent energy forms and three balance laws are established in their respective domains. Reiher etal. [START_REF] Reiher | Finite third-order gradient elasticity and thermoelasticity[END_REF][START_REF] Reiher | Finite third-order gradient elastoplasticity and thermoplasticity[END_REF] developed a finite third-strain gradient elasticity/elastoplasticity theory.

The higher order continuity in gradient elasticity theory imposes challenges on many numerical methods. In order to satisfy the C 1 or C 2 continuity, a variety of numerical methods have been developed, see for instance the mixed finite element method proposed in [START_REF] Shu | Finite elements for materials with strain gradient effects[END_REF][START_REF] Amanatidou | Mixed finite element formulations of strain-gradient elasticity problems[END_REF][START_REF] Askes | Implicit gradient elasticity[END_REF][START_REF] Jc Reiher | Finite-element analysis of polyhedra under point and line forces in second-strain gradient elasticity[END_REF], boundary element method [START_REF] Atroshchenko | Fundamental solutions and dual boundary element methods for fracture in plane cosserat elasticity[END_REF], meshless methods [START_REF] Askes | Numerical modeling of size effects with gradient elasticity-formulation, meshless discretization and examples[END_REF][START_REF] Tang | Analysis of materials with strain-gradient effects: a meshless local petrov-galerkin (mlpg) approach, with nodal displacements only[END_REF], isogeometric analysis (IGA) formulations [START_REF] Balobanov | Isogeometric analysis of gradient-elastic 1d and 2d problems[END_REF][START_REF] Khakalo | Isogeometric analysis of higher-order gradient elasticity by user elements of a commercial finite element software[END_REF][START_REF] Niiranen | Isogeometric analysis for sixth-order boundary value problems of gradient-elastic Kirchhoff plates[END_REF][START_REF] Cl Thanh | The size-dependent thermal bending and buckling analyses of composite laminate microplate based on new modified couple stress theory and isogeometric analysis[END_REF][START_REF] Phung-Van | Porosity-dependent nonlinear transient responses of functionally graded nanoplates using isogeometric analysis[END_REF][START_REF] H X Nguyen | A refined quasi-3d isogeometric analysis for functionally graded microplates based on the modified couple stress theory[END_REF] and nonlocal operator method (NOM) [START_REF] Hl Ren | Nonlocal operator method with numerical integration for gradient solid[END_REF][START_REF] Hl Ren | A higher order nonlocal operator method for solving partial differential equations[END_REF][START_REF] Hl Ren | A nonlocal operator method for solving partial differential equations[END_REF]. NOM is proposed as a generalization of dual-horizon peridynamics [START_REF] Hl Ren | Dual-horizon peridynamics[END_REF][START_REF] Hl Ren | Dual-horizon peridynamics: A stable solution to varying horizons[END_REF][START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF][START_REF] Silling | Peridynamic states and constitutive modeling[END_REF]. It uses an integral form (i.e.nonlocal operators) to replace the local partial differential operators of various orders. The nonlocal operators can be viewed as an alternative to the partial derivatives of shape functions in FEM. Combined with a variational principle or weighted residual method, NOM obtains the residual vector and tangent stiffness matrix in the same way as in FEM. There are three versions of NOM, first-order particle-based NOM [START_REF] Hl Ren | A nonlocal operator method for solving partial differential equations[END_REF][START_REF] Rabczuk | A nonlocal operator method for partial differential equations with application to electromagnetic waveguide problem[END_REF], higher order particle-based NOM [START_REF] Hl Ren | A higher order nonlocal operator method for solving partial differential equations[END_REF] and higher order NOM based on numerical integration [START_REF] Hl Ren | Nonlocal operator method with numerical integration for gradient solid[END_REF]. The particle-based version can be viewed as a special case of NOM with numerical integration when nodal integration is employed. NOM has been applied to the solution of the Poisson equations in high dimensional space, von-Karman thin plate equations, fracture problems based on phase fields [START_REF] Hl Ren | A higher order nonlocal operator method for solving partial differential equations[END_REF], waveguide problem in electromagnetic field [START_REF] Rabczuk | A nonlocal operator method for partial differential equations with application to electromagnetic waveguide problem[END_REF] and gradient solid problems [START_REF] Hl Ren | Nonlocal operator method with numerical integration for gradient solid[END_REF]. NOM is suitable for problems requiring higher order continuity though its application to higher order gradient elasticity has not been explored yet.

Current higher order gradient solid theory is limited to third order. Javili etal. [START_REF] Javili | Geometrically nonlinear higher-gradient elasticity with energetic boundaries[END_REF] implemented the variational derivation of the third-gradient elasticity without considering the specific forms of energy. The derivation follows the setting of first Piola-Kirchhoff stress.

The concrete form of third-strain gradient theory with finite deformation is proposed in [START_REF] Reiher | Finite third-order gradient elasticity and thermoelasticity[END_REF],

where the objectivity of the energy form is emphasized. However, an associated implementation is missing to our best knowledge. In this paper, we propose a different strain energy density with objectivity. The energy form is based on the second Piola-Kirchhoff stress and is invariant under rigid body transformation. The number of gradient order is extended to 5 in 2D and 3 in 3D. For the first time, the geometrical nonlinear fifth-order gradient elasticity in 2D and third-order gradient elasticity in 3D are studied by numerical experiments based on nonlocal operator method.

The content of the paper is outlined as follows. The general strain energy density for large deformations is proposed for the fourth-order gradient elasticity in section 2. In Section 3, we derive the governing equations and the associated boundary conditions for the thirdorder gradient elasticity by using variational principles and exploiting integration by parts on surfaces. In section 4, the framework of the particle-based nonlocal operator method is briefly summarized and its implementation for solving higher order gradient solids presented.

In section 5, several representative numerical tests, including a point displacement load, point force load and the influence of the length scale in linear/nonlinear gradient elasticity, are presented to study the physical response of higher order gradient elasticity. Finally, we conclude our manuscript in section 6.

Higher order gradient solid with finite deformation

Let us denote the material coordinates (in the initial configuration Ω) by X, the spatial coordinates (in the current configuration Ω t ) by x and the displacement field by u := x-X.

The deformation gradient F , right Cauchy Green tensor C and Green-Lagrange strain tensor E are written as

F = ∂x ∂X = ∇x = ∇u + I (1)
C = F T • F , E = 1 2 (C -I) ( 2
)
where I is the identity matrix and ∇ = ∂ ∂X . The principle of frame indifference requires the quantity remain invariant under rigid body transformation x = Q(t)x+c(t) where c(t) is the rigid translation vector, Q(t) the (orthogonal) rotation matrix satisfying

QQ T = Q T Q = I.
Vectors and second-order tensors are objective if they are related by the rotation tensor as

u = Q • u (3)
T = Q • T • Q T (4)
The deformation gradient under rigid body transformation is related by

F = ∂x ∂X = ∂x ∂x • ∂x ∂X = Q • F (5)
Quantities are invariant if they remain unchanged by the rigid body transformation. Apparently, the right Cauchy tensor is invariant as

C = F T F = F T Q T QF = F T F = C.
Let F ,i := ∂F ∂X i denote the partial derivative of F with respect to X i . The derivative of F and Q can be written as

F ,i = Q ,i F + QF ,i (6)
Q T ,i Q + Q T Q ,i = 0 (7)
The gradient of C can be derived as

C ,i = (F T F ) ,i = F T ,i F + F T F ,i = (F T Q T ,i + F T ,i Q T )QF + F T Q T (Q ,i F + QF ,i ) = F T Q T ,i QF + F T ,i Q T QF + F T Q T Q ,i F + F T Q T QF ,i = -F T Q T Q ,i F + F T Q T Q ,i F + F T F ,i + F T ,i F = (F T F ) ,i = C ,i . (8)
Therefore ∇C is invariant. Note that ∇F is neither objective nor invariant since

F ,i = (QF ) ,i = Q ,i F + QF ,i = F ,i (9)
∇F :∇F = i F ,i : F ,i = i F ,i : F ,i . (10)
Hence ∇F cannot be used directly to define the energy density. However, the invariant property makes ∇C a good choice. Different orders of strain gradient can be written as

H = ∂E ∂X = ∇E = 1 2 (F T ∇F + ∇F T F ) (11)
for the first gradient strain tensor,

G = ∂H ∂X = ∂ 2 E ∂X 2 = ∇ 2 E (12)
for the second-gradient strain tensor, and

L = ∂G ∂X = ∂ 3 E ∂X 3 = ∇ 3 E (13)
for the third-gradient strain tensor and with ∇ n = ∇ ⊗ ... ⊗ ∇ n times

. Most 'computational' contributions focus on applications with small strain gradients. In this case, strain gradients of different orders are decoupled. Note that the expression in L contains derivatives ranging from ∇u to ∇ 4 u, which is due to geometric nonlinearities arising from finite deformation. The exact forms of L, G can be derived with Mathematica. It can be verified that C, E, H, G, L are invariants under rigid body translations and rotations. Thus, these quantities can be used to define the stresses and generalized stresses

S ij = D ijkl E kl with D ijkl = ∂S ij ∂E kl (14a) R ijk = E ijklmn H lmn with E ijklmn = ∂R ijk ∂H lmn (14b) Q ijkl = F ijklmnst G mnst with F ijklmnst = ∂Q ijkl ∂G mnst (14c)
P ijklm = G ijklmαβγηδ L αβγηδ with G ijklmαβγηδ = ∂P ijklm ∂L αβγηδ (14d)
Eq.14a has the exact form of the Saint Venant-Kirchhoff model. R ijk , Q ijkl and P ijklm are defined as the generalization of the Saint Venant-Kirchhoff model.

The strain energy density in the initial configuration can be assumed as

φ = 1 2 (S ij E ij + R ijk H ijk + Q ijkl G ijkl + P ijklm L ijklm ) (15) = 1 2 (E ij D ijkl E kl + H ijk E ijklmn H lmn + G ijkl F ijklmnst G mnst + L ijklm G ijklmαβγηδ L αβγηδ ) ( 16
)
where D is a 4th-order tensor, E a 6th-order tensor, F an 8th-order tensor and G a 10th-order tensor. The second Piola-Kirchhoff stress is work conjugate to the Green-Lagrange strain [START_REF] Bonet | Nonlinear continuum mechanics for finite element analysis[END_REF][START_REF] Korelc | Automation ofFinite Element Methods[END_REF]. Therefore, the generalized second Piola-Kirchhoff stresses define the strain energy density in the initial configuration and are objective under any rigid body transformations.

The strain energy density, given by Eq.16, is among the most simplified quadratic energy functionals. This energy functional is required to be positive, which can be satisfied provided that the material tensors D, E, F and G are positive definite. For an n dimensional space, a kth-order tensor has n k entries, for example G has 2 10 = 1 024 elements in 2D and 3 10 = 59 049 elements in 3D. However, when symmetry conditions are exploited, the number of elements can be greatly reduced. We discussed the symmetry of a 6th-order tensor in Appendix A and in Voigt notation in Appendix B. For a third-order gradient solid, there are thousands of material parameters to be determined, which are complicated to resolve experimentally. For simplicity, we introduce only three material length scales.

φ(E, ∇E, ∇ 2 E, ∇ 3 E) = 1 2 S : E + l 2 1 ∇S:∇E + l 4 2 ∇ 2 S • (4) ∇ 2 E + l 6 3 ∇ 3 S • (5) ∇ 3 E (17)
where • (k) is the generalization of the inner product, for example,

• (1) = •, • (2) =:, • (3)
= :, and the stresses can be written as

S = D : E (18a) R = ∇S (18b) Q = ∇ 2 S (18c)
P = ∇ 3 S. (18d)
The total internal strain energy in the domain can be expressed as

F int = Ω φ(E, ∇E, ∇ 2 E, ∇ 3 E) (19)
We recall that Ω is the initial configuration. Here we used Ω {•} := Ω {•} dV . When small deformations are assumed, the Green-Lagrange strain and second Piola-Kirchhoff stress degenerate to the linear strain and Cauchy stress tensor:

E → ε = 1 2 (∇u) T + ∇u (20)
S → σ = D : ε (21)
Then Eq.17 can be written as

φ(ε, ∇ε, ∇ 2 ε, ∇ 3 ε) = 1 2 (σ : ε + l 2 1 ∇σ:∇ε + l 4 2 ∇ 2 σ • (4) ∇ 2 ε + l 6 3 ∇ 3 σ • (5) ∇ 3 ε) (22)
Based on the simplified higher gradient elasticity in Eq.17, the general nth-gradient elasticity can be written as

φ(E, ∇E, ∇ 2 E, ..., ∇ n E) = 1 2 S : E + n k=1 (l k ) 2k ∇ k S • (k+2) ∇ k E ( 23
)
where l k is the internal length scale of kth-order. We employ E n elasticity (n = 0, 1, 2, 3, 4, 5)

to abbreviate the nth-gradient hyperelasticity theory wherein ∇ n E is used to define the energy density functional. Accordingly, the conventional hyperelasticity is denoted by E 0 elasticity, gradient elasticity is abbreviated by E 1 elasticity and the second-gradient elasticity by E 2 elasticity. The highest order of the energy form of E n -elasticity is (n + 1), while that of the strong form is 2(n + 1). The highest order of gradient elasticity we implemented is the E 5 elasticity in 2D. For E 5 elasticity, the governing equations are a set of 12th-order nonlinear partial differential equations (PDEs).

3. Governing equations of second-gradient solid

Integration by parts on close surface

Before delving into the variational derivation of second-gradient solids, we briefly review the integration by parts in domains and on surfaces using the following abbreviations

∂Ω i ∂ 2 Ω ij Ω n m γ
Ω {•} := Ω {•} dV and ∂Ω {•} := ∂Ω {•} dΓ.
The integration by parts for tensor fields in Ω is

Ω S : ∇u = ∂Ω n • S • u - Ω ∇ • S • u ( 24
)
where S is the second-order tensor field and u is the vector field.

According to Refs [START_REF] Gurtin | A continuum theory of elastic material surfaces[END_REF][START_REF] Sp Shen | A theory of flexoelectricity with surface effect for elastic dielectrics[END_REF], the integration by parts under the assumption of smooth surfaces can be expressed as

∂Ω S : ∇u = ∂Ω S : ∇ n u + ∂Ω S : ∇ t u = ∂Ω S : ∇ n u + (gn • S -∇ t • S) • u (25)
where g = ∇ t • n is the Gauss curvature, ∇ t and ∇ n are the tangential and normal gradient

∇ t = (1 -nn T ) • ∇, ∇ n = nn T • ∇ (26)
In the above derivation, ∇u needs to be divided into the tangential and normal parts with respect to the surface, i.e.∇u = ∇ t u + ∇ n u. For the case of piecewise smooth surfaces shown in Fig. 1, the boundary term should be considered explicitly,

∂Ω S : ∇u = ∂Ω S : ∇ n u + (gn • S -∇ t • S) • u + ∂ 2 Ω m • S • u ( 27
)
where m is the outward normal direction of ∂Ω in the tangent plane defined by n. If both S and ∇u are defined based on the tangent space of the surface, i.e. n • S = 0 and ∇ n u = 0, the integration by parts on surface is the same as Eq.24.

Variational derivation of second-gradient solid

In this paper, we only consider the higher order bulk energy. For the boundary (surface and curve) energies, the reader is referred to [START_REF] Javili | Geometrically nonlinear higher-gradient elasticity with energetic boundaries[END_REF] for more details. The second-gradient solid for linear elasticity with second velocity gradient inertia can be found in [START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part i: Constitutive equations and quasi-static behavior[END_REF][START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part ii: Dynamic behavior[END_REF]. Let φ := φ(u, ∇u, ∇ 2 u, ∇ 3 u) denote the internal energy density of a second-gradient solid (E 2 elasticity). The variation of the internal energy in Ω is then given by

δF int = Ω ∂φ ∂u b •δu + ∂φ ∂∇u S 1 : ∇δu + ∂φ ∂∇ 2 u S 2 :∇ 2 δu + ∂φ ∂∇ 3 u S 3 :: ∇ 3 δu = Ω b • δu p1 + S 1 : ∇δu p2 + S 2 :∇ 2 δu p3 + S 3 :: ∇ 3 δu p4 ( 28
)
where b can be viewed as the body force density, S i denotes the work conjugate to ∇ i u, (i = 1, 2, 3) and S i the general stress defined in the initial configuration:

S 3 = ∂φ ∂∇ 3 u = ∂φ ∂∇ 2 E :: ∂∇ 2 E ∂∇ 3 u (29a) S 2 = ∂φ ∂∇ 2 u = ∂φ ∂∇E : ∂∇E ∂∇ 2 u + ∂φ ∂∇ 2 E :: ∂∇ 2 E ∂∇ 2 u (29b) S 1 = ∂φ ∂∇u = ∂φ ∂E : ∂E ∂∇u + ∂φ ∂∇E :∂∇E ∂∇u + ∂φ ∂∇ 2 E :: ∂∇ 2 E ∂∇u (29c)
where E, ∇E, ∇ 2 E have the forms

E = 1 2 ∇u T + ∇u + ∇u T ∇u (30a) ∇E = 1 2 ∇(∇u) T + ∇ 2 u + ∇(∇u) T ∇u + ∇u T ∇ 2 u (30b) ∇ 2 E = 1 2 ∇ 2 (∇u) T + ∇ 3 u + ∇ 2 (∇u) T ∇u+ ∇(∇u) T ∇ 2 u + ∇(∇u) T ∇ 2 u + (∇u) T ∇ 3 u (30c)
Obviously, higher order strain gradients contain a low order gradient due to the existence of a geometric nonlinearity. The explicit form of S 1 , S 2 , S 3 can be derived by mathematical software such as Mathematica. Since the order of the derivatives is reduced by one each time when applying integration by parts, we integrate the terms from high order to low order so that the derived low order term with the existing low order term can be handled together.

In other words, the term with the third-order gradient is firstly integrated by parts, then the second-order gradient term and the accumulated gradient terms at last.

The p4 part in Eq.28 via Eq.24 can be written as

Ω S 3 :: ∇ 3 δu = ∂Ω n • S 3 :∇ 2 δu + Ω -∇ • S 3 :∇ 2 δu = ∂Ω n • S 3 :∇ t (∇δu) p7 + ∂Ω n • S 3 :∇ n (∇δu) p8 + Ω -∇ • S 3 :∇ 2 δu p6 (31)
The p7 part in Eq.31 is obtained via integration by parts on surface

∂Ω n • S 3 :∇ t (∇δu) = ∂Ω gn • (n • S 3 ) -∇ t • (n • S 3 ) V 1 : ∇δu = ∂Ω V 1 : ∇δu p9 (32)
In above derivation, we abbreviate the long expression as V 1 . Based on Eq.25, the p9 part reads

∂Ω V 1 : ∇δu = ∂Ω V 1 : ∇ n δu + (gn • V 1 -∇ t • V 1 ) • δu (33)
Hence, p4 has the form

Ω S 3 :: ∇ 3 δu = ∂Ω V 1 : ∇δu + ∂Ω n • S 3 :∇ n (∇δu) p8 + Ω -∇ • S 3 :∇ 2 δu p6 (34)
The p6 term has the same form as part p3 in Eq.28 and it can be estimated that a new surface term similar to p9 part will arise. The summation of p 3 , p 6 , p 9 can be written as

Ω S 2 :∇ 2 δu -∇ • S 3 :∇ 2 δu + ∂Ω V 1 : ∇δu = Ω (S 2 -∇ • S 3 ):∇ 2 δu
by Eq.24

+ ∂Ω V 1 : ∇δu = ∂Ω n • (S 2 -∇ • S 3 ) : ∇δu + ∂Ω V 1 : ∇δu by Eq.25 - Ω ∇ • (S 2 -∇ • S 3 ) : ∇δu = ∂Ω (gn • V 1 -∇ t • V 1 + gn • (n • (S 2 -∇ • S 3 )) -∇ t • (n • (S 2 -∇ • S 3 ))) V 2
•δu

+ ∂Ω (V 1 + n • (S 2 -∇ • S 3 )) : ∇ n δu + Ω -∇ • (S 2 -∇ • S 3 ) : ∇δu p10 = ∂Ω V 2 • δu + ∂Ω (V 1 + n • (S 2 -∇ • S 3 )) : ∇ n δu + Ω -∇ • (S 2 -∇ • S 3 ) : ∇δu (35a)
In the above derivation, the long expression in the second line is abbreviated as V 2 . Both integration by parts in the domain and integration by parts on the surface are used in the derivation.

The gradient order of p10 is identical to that in the p2 part. For simplicity, the integration by parts of p10 + p2 is

Ω (S 1 -∇ • (S 2 -∇ • S 3 )) S 4 : ∇δu = ∂Ω n • S 4 • δu - Ω ∇ • S 4 • δu (36)
By using integration by parts in Eq.24 and Eq.25 several times, the variation of the internal energy can be finally written as

δF int = Ω (b -∇ • S 4 ) • δu + ∂Ω (V 1 + n • (S 2 -∇ • S 3 )) : ∇ n δu + ∂Ω V 2 • δu + ∂Ω n • S 4 • δu + ∂Ω n • S 3 :∇ n (∇δu) (37)
=

Ω (b -∇ • S 4 ) • δu + ∂Ω (V 2 + n • S 4 ) • δu + ∂Ω (V 1 + n • (S 2 -∇ • S 3 )) : ∇ n δu + ∂Ω n • S 3 :∇ n (∇δu) (38)
where

S 4 = S 1 -∇ • (S 2 -∇ • S 3 ) (39a) V 1 = gn • (n • S 3 ) -∇ t • (n • S 3 ) (39b) V 2 = gn • V 1 -∇ t • V 1 + gn • (n • (S 2 -∇ • S 3 )) -∇ t • (n • (S 2 -∇ • S 3 )) (39c)
In Eq.38, we used ∂Ω to denote the boundaries, which should be tailored based on the actual boundary conditions. The variation of the internal energy yields the work-conjugate pairs on the boundaries. The expression for V 2 describes the contribution from the curvature related terms (e.g. curvature and curvature gradient on surface) as well as the generalized stresses of different orders. The contribution from the surface curvature indicates that the gradient effect or nonlocal effect of the solid may be significant at the 'micro'-scale, where the surface-to-volume ratio is much larger compared to the macro scale and the surface curvature for small-scale objects is huge. This relation is consistent with the fact that the strength of a material at micro-scale is much larger than that at macro-scale.

For the gradient elasticity, the boundary conditions may contain essential boundary conditions such as translation, the gradient of the translation and force boundary conditions like stress and couple stress. The gradient of the translation is similar to the prescribed rotation on the boundary in plate/shell theory, while the couple stress is the work conjugate to the gradient of the translation.

Based on the boundary work-conjugate pairs in the variation of the internal energy, the external energy can be constructed as

F ext = ∂Ω 0 D P • (u -ū) + ∂Ω 0 N P • u + ∂Ω 1 D Q : (∇ n u -∇ n u) + ∂Ω 1 N Q : ∇ n u + ∂Ω 2 D R:(∇ n ∇u -∇ n ∇u) + ∂Ω 2 N R:∇ n ∇u (40)
where

P = (V 2 + n • S 4 ), Q = (V 1 + n • (S 2 -∇ • S 3 )), R = n • S 3 , ∂Ω i D , ∂Ω i N , (i = 0, 1,
2) refer to the Dirichlet and Neumann boundary conditions for u of different partial derivative orders; ∂Ω 0 D designates the constraints of the displacement, ∂Ω 1 D denotes the constraints of the displacement gradient (e.g. fixed rotation state) and ∂Ω 1 D describes the displacement second-gradient; P is the traction load, Q refers to the couple stress load and R is the higher order couple stress load.

Gradient elasticity deals not only with the gradient of the deformation but also the gradient of inertia terms [START_REF] R D Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Askes | Gradient elasticity and flexural wave dispersion in carbon nanotubes[END_REF][START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part i: Constitutive equations and quasi-static behavior[END_REF][START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part ii: Dynamic behavior[END_REF]. The kinetic energy with velocity gradient can be written as [START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part ii: Dynamic behavior[END_REF]

K = t 1 t 0 Ω 1 2 ρ u • u + 1 2 ρl 2 d ∇ u : ∇ u ( 41
)
The variation of the kinetic energy can be written as

δK = t 1 t 0 Ω ρ u • δ u + ρl 2 d ∇ u : ∇δ u = t 1 t 0 Ω -ρ ü • δu -ρl 2 d ∇ ü : ∇δu = t 1 t 0 Ω -ρ ü • δu + ρl 2 d ∇ • ∇ ü • δu (42)
For any δu, δ∇ n u, δ∇ n ∇u, the Hamilton principle δK -

t 1 t 0 δF int + t 1 t 0 δF ext = 0 (43) leads to ρ ü -l 2 d ∇ • ∇ ü = -b + ∇ • (S 1 -∇ • (S 2 -∇ • S 3 )) in Ω (44a) u = ū on ∂Ω 0 D (44b) P = P on ∂Ω 0 N (44c) ∇ n u = ∇ n u on ∂Ω 1 D (44d) Q = Q on ∂Ω 1 N (44e) ∇ n ∇u = ∇ n ∇u on ∂Ω 2 D (44f) R = R on ∂Ω 2 N (44g)
The derivation of energies based on variational principle leads very naturally to the governing equations and various boundary conditions. The maximal order of derivatives in Eq.44a and V 2 is 6, 5, respectively. Similarly, the variation of the strain energy density of strain gradient elasticity (E 1 elasticity) can be derived. By setting S 3 = 0 in Eq.44 , we obtain the governing equations and boundary conditions of gradient elasticity

ρ ü -l 2 d ∇ • ∇ ü = -b + ∇ • (S 1 -∇ • S 2 ) in Ω (45a) u = ū on ∂Ω 0 D ( 45b
)
P = P on ∂Ω 0 N (45c) ∇ n u = ∇ n u on ∂Ω 1 D (45d) Q = Q on ∂Ω 1 N (45e)
where

Q = n • S 2 , P = (gn • (n • S 2 ) -∇ t • (n • S 2 ) + n • (S 1 -∇ • S 2 )), ∂Ω i D , ∂Ω i N , (i = 0,

Numerical implementation

Review of Nonlocal Operator Method

NOM uses the integral form to replace the partial differential derivatives of different orders. We adopted a Total Lagrangian description of motion for the higher order gradient elasticity NOM. Consider a domain as shown in Fig. 2(a), let X i be spatial coordinates in the domain Ω; r := X j -X i is a spatial vector ranging from X i to X j ; v i := v(X i , t) and v j := v(X j , t) are the field values for X i and X j , respectively; v ij := v j -v i is the relative field vector for spatial vector r. Support S i of point X i is the neighborhood of point X i . A point X j in support S i forms the spatial vector r(= X j -X i ). The support in the NOM can be a spherical domain, a cube, semi-spherical domain and so on.

Dual-support is defined as a union of points whose supports include X, denoted by

S i = {X j |X i ∈ S j }. ( 46
) x i x j � v(X i ) v(X j ) o r =X j -X i S i (a) S x x 1 x 2 x 3 x 5 x x 4 (b)
S X = {X 1 , X 2 , X 4 }, S X = {X 1 , X 2 , X 3 }.
Point X j forms the dual-vector r (= X i -X j = -r) in S i . On the other hand, r is the spatial vector formed in S j . One example to illustrate the support and dual-support is shown in Fig. 2(b).

The first-order nonlocal operator method uses the basic nonlocal operators to replace the local operator in calculus such as the gradient, divergence and curl operators. The functional formulated by the local differential operator can be used to construct the residual or tangent stiffness matrix by replacing the local operator with the corresponding nonlocal operator.

The nonlocal gradient of a vector field v for point X i in support S i is defined as

∇v i := S i w(r)v ij ⊗ r dV j • S i w(r)r ⊗ r dV j -1 . ( 47
)
The nonlocal gradient operator and its variation in discrete form are given by

∇v i = j∈S i w(r)v ij ⊗ r∆V j • j∈S i w(r)r ⊗ r∆V j -1 , (48)
∇δv i = j∈S i w(r)δv ij ⊗ r∆V j • j∈S i w(r)r ⊗ r∆V j -1 . ( 49
)
The operator energy functional for a vector field at point x i is

F hg i = p hg S i w(r)( ∇v i • r -v ij ) • ( ∇v i • r -v ij ) dV j ( 50
)
where p hg is a penalty coefficient. The residual and tangent stiffness matrix of F hg i can be obtained with ease, see [START_REF] Hl Ren | A nonlocal operator method for solving partial differential equations[END_REF] for more details. For problems that require higher order continuity, the higher order NOM is needed. According to Ref [START_REF] Hl Ren | A higher order nonlocal operator method for solving partial differential equations[END_REF], a scalar field u j at a point j ∈ S i can be obtained by a Taylor series expansion at u i in d dimensions with maximal derivative order not higher than n,

u j = u i + (n 1 ,...,n d )∈α n d r n 1 1 ...r n d d n 1 !...n d ! u i,n 1 ...n d + O(r |α|+1 ) (51)
with

r = (r 1 , ..., r d ) = (X j1 -X i1 , ..., X jd -X id ) (52a) u i,n 1 ...n d = ∂ n 1 +...+n d u i ∂X n 1 i1 ...∂X n d id (52b) |α| = max (n 1 + ... + n d ) (52c)
α n d being the list of multi-indexes, given by

α n d = {(n 1 , ..., n d )|1 ≤ d i=1 n i ≤ n, n i ∈ N 0 , 1 ≤ i ≤ d} (53)
and N 0 = {0, 1, 2, 3, ...}. The number of multi-indices in α n d is (n + d)!/(n!d!) -1 and all elements in α n d of Eq.53 can be obtained easily by Mathematica [START_REF] Hl Ren | A higher order nonlocal operator method for solving partial differential equations[END_REF]. For any multi-index (n 1 , ..., n d ) ∈ α n d , the partial derivative and the polynomial are

u i,n 1 ...n d , r n 1 1 ...r n d d n 1 !...n d ! , ∀(n 1 , ..., n d ) ∈ α n d . (54)
When the length scale of support S i at u i is taken into account, the Taylor series expansion in Eq.51 can be written as

u j = u i + (n 1 ,...,n d )∈α n d r n 1 1 ...r n d d h n 1 +...+n d i h n 1 +...+n d i n 1 !...n d ! u i,n 1 ...n d + O(r n+1 ) = u i + (n 1 ,...,n d )∈α n d r n 1 1 ...r n d d h n 1 +...+n d i u h i,n 1 ...n d + O(r n+1 ) ( 55
)
where h i is the characteristic length of S i , and

u h i,n 1 ...n d = h n 1 +...+n d i n 1 !...n d ! u i,n 1 ...n d (56)
Let p h j , ∂ h α u i and ∂ α u i denote the list of the polynomials, scaled partial derivatives, partial derivatives, respectively, based on multi-index notation α n d in Eq.53,

p h j = ( r d h , ..., r n 1 1 ...r n d d h n 1 +...+n d , ..., r n 1 h n ) T (57a) ∂ h α u i = (u h i,0...1 , ..., u h i,n 1 ...n d , ..., u h i,n...0 ) T (57b)
∂ α u i = (u i,0...1 , ..., u i,n 1 ...n d , ..., u i,n...0 ) T . (57c)
∂ h α u i and ∂ α u i are related by

∂ α u i = H -1 i ∂ h α u i , with H i = diag h i , ..., h n 1 +...+n d i n 1 !...n d ! , ..., h n i n! (58)
where diag[a 1 , ..., a n ] denotes a diagonal matrix whose diagonal entries starting in the upper left corner are a 1 , ..., a n . Therefore, the Taylor series expansion with u i being moved to the left side of the equation can be written as

u ij = (∂ h α u i ) T p h j , ∀j ∈ S i (59)
where u ij = u j -u i . Integrating u ij with weighted coefficient w(r)(p h j ) T in support S i , we obtain

S i w(r)u ij (p h j ) T dV j = (∂ h α u i ) T S i w(r)p h j ⊗ (p h j ) T dV j = (∂ α u i ) T H i S i w(r)p h j ⊗ (p h j ) T dV j ( 60
)
where w(r) is the weight function. Thus, the nonlocal operator ∂α u i can be obtained as

∂α u i := H -1 i S i w(r)p h j ⊗ (p h j ) T dV j -1 S i w(r)u ij p h j dV j = K i • S i w(r)p h j u ij dV j (61)
with

K i := H -1 i S i w(r)p h j ⊗ (p h j ) T dV j -1 . ( 62
)
The variation of ∂α u i is

∂α δu i := K i • S i w(r)p h j (δu j -δu i ) dV j (63)
In the continuous form, the number of dimensions of ∂δu i is infinite and a discretization is required. After discretization of the domain by particles, the whole domain is represented by

Ω = N i=1 ∆V i ( 64
)
where i is the global index of volume ∆V i and N is the number of particles in Ω. Particles in S i are represented by

S i = {j 1 , ..., j k , ..., j n i } (65)
where j 1 , ..., j k , ..., j n i are the global indices of neighboring particle i and n i is the number of neighbors of i in S i . The discrete form of Eq.61 and its variation are

∂α u i = K i • j∈S i u ij w(r j )p h j ∆V j = K i p h wi ∆u i (66a) ∂α δu i = K i • j∈S i δu ij w(r j )p h j ∆V j = K i p h wi δ∆u i (66b)
with

K i = H -1 i j∈S i w(r)p h j ⊗ (p h j ) T ∆V j -1 , (67a)
p h wi = w(r j 1 )p h j 1 ∆V j 1 , ..., w(r jn i )p h jn i ∆V jn i (67b)

∆u i = (u ij 1 , ..., u ij k , ..., u ijn i ) T (67c)
The nonlocal operator provides all partial derivatives with maximal order for a single index up to n. The set of derivatives in PDEs of real application is a subset of the nonlocal operator. Together with the weak formulation (weighted residual method or variational principles (i.e. [START_REF] Hl Ren | A nonlocal operator method for solving partial differential equations[END_REF])), Eq.66a can be employed to solve many linear (nonlinear) PDEs. Eq.66a can be written more concisely as

∂α u i = K i p h wi ∆u i = B αi u i (68)
with B αi being the operator matrix for point i based on multi-index α n d

B αi =   -(1, • • • , 1) np K i p h wi K i p h wi   ( 69
)
u i = (u i , u j 1 , u j 2 , • • • , u jn i ) T (70)
where (1,

• • • , 1) np K i p h wi is the column sum of K i p h wi , n p is the length of α n d .
The operator matrix obtains all partial derivatives of maximal order less than |α| + 1 by the nodal values in the support. For real applications, one can select the specific rows in the operator matrix based on the partial derivatives contained in the specific PDEs. For example, if the order of derivatives in the given PDEs are ∂u i = (u i,Y , u i,Y Y , u i,XY ) T ⊂ ∂α u i , one can select the blue lines in Eq.71 to form the actually operator matrix given by B i in Eq.72.

              u i,Y u i,X u i,Y Y u i,XY u i,XX . . .               ∂αui =               b 11 b 12 • • • b 1(n+1) b 21 b 22 • • • b 2(n+1) b 31 b 32 • • • b 3(n+1) b 41 b 42 • • • b 4(n+1) . . . . . . . . . . . . b m1 b m2 • • • b m(n+1)               B αi         u i u j 1 . . . u jn         u i (71)      u i,Y u i,Y Y u i,XY      ∂u i =      b 11 b 12 • • • b 1(n+1) b 31 b 32 • • • b 3(n+1) b 41 b 42 • • • b 4(n+1)      B i         u i u j 1 . . . u jn         u i (72)
For a given maximal differential order and number of space dimensions, NOM offers the derivatives of all orders in discrete form 'automatically'. These nonlocal derivatives are similar to the derivatives of the shape functions in IGA. When the selected derivatives are inserted into the equivalent functional of the physical problem in discrete form, the residual and tangent stiffness matrix of the functional can be derived.

Besides considering the functional for the physical problem, the functional for the nonlocal operators should be considered explicitly. The energy functional for all nonlocal operators is defined as [START_REF] Hl Ren | A higher order nonlocal operator method for solving partial differential equations[END_REF] F i (u) =

j∈S i w(r) u ij -(p h j ) T ∂h α u i 2 ∆V j (73)
Based on Eq.66a, F i (u) can be simplified as

F i (u) = j∈S i w(r)u 2 ij ∆V j -∆u T i (p h wi ) T j∈S i w(r)p h j (p h j ) T ∆V j -1 p h wi ∆u i =∆u T i diag w(r j 1 )∆V j 1 , ..., w(r jn i )∆V jn i -(p h wi ) T j∈S i w(r)p h j (p h j ) T ∆V j -1 p h wi ∆u i =∆u T i M i ∆u i (74)
with

M i = diag w(r j 1 )∆V j 1 , ..., w(r jn i )∆V jn i -(p h wi ) T j∈S i w(r)p h j (p h j ) T ∆V j -1 p h wi (75)
Apparently, M i is a symmetric matrix. The expression of F i (u) is quadratic, and its Hessian matrix can be extracted as

K hg i = p hg m i   v i -v T i -v i M i   (76)
where v i (j) = n i k=1 M i (j, k) is the sum of the row of M i ; the first row (column) denotes the entries for point i, while the neighbors start from the second row (column), p hg is a penalty coefficient and m i the normalization coefficient given by m i = j∈S i w(r)r • r∆V j .

The reader is referred to [START_REF] Hl Ren | A higher order nonlocal operator method for solving partial differential equations[END_REF] for more details of the NOM.

Newton-Raphson method

The governing equations and boundary conditions in Eq.44 are quite complicated. The highest continuity in Q is C 4 and the gradient and Hessian matrix of the functional on boundary ∂Ω 1 D are cumbersome. Note that NOM does not satisfy the Kronecker-Delta property, and the order of NOM should be at least C 5 in order to satisfy the Dirichlet boundary conditions on ∂Ω 0 D , where the continuity order in P is C 5 . Therefore, we employ the penalty method to enforce both Dirichlet boundary conditions and the normal Dirichlet boundary conditions. The equivalent energy functional of second-gradient elasticity then becomes

F = Ω φ(u, ∇u, ∇ 2 u, ∇ 3 u) dV + ∂Ω 0 D α 1 (u -ū) • (u -ū) dS - ∂Ω 0 N P • u dS + ∂Ω 1 D α 2 (∇ n u -∇ n u) : (∇ n u -∇ n u) dS - ∂Ω 1 N Q : ∇ n u dS (77)
where α 1 , α 2 are penalty parameters. One advantage of the penalty method is that the highest order of partial derivatives is 4 for third-gradient elasticity, while the formulation based on the modified variational principle requires C 7 continuity. We neglect the terms on ∂Ω 2 D and ∂Ω 2 N for simplicity. After discretization, the discrete form of the functional in Eq.77 becomes

F = ∆V i ∈Ω φ(u i , ∇u i , ∇ 2 u i , ∇ 3 u i )∆V i + ∆S i ∈∂Ω 0 D α 1 (u i -ūi ) • (u i -ūi )∆S i - ∆S i ∈∂Ω 0 N Pi • u i ∆S i + ∆S i ∈∂Ω 1 D α 2 (∇ n u i -∇ n u i ) : (∇ n u i -∇ n u i )∆S i - ∆S i ∈∂Ω 1 N Qi : ∇ n u i ∆S i (78)
The differential derivatives in φ of E 3 elasticity with unknowns u = (u, v) T in material coordinates X = (X, Y ) are

∂U 2d = u ,Y , v ,Y , u ,X , v ,X , u ,Y Y , v ,Y Y , u ,XY , v ,XY , u ,XX , v ,XX , u ,Y Y Y , v ,Y Y Y , u ,XY Y , v ,XY Y , u ,XXY , v ,XXY , u ,XXX , v ,XXX , u ,Y Y Y Y , v ,Y Y Y Y , u ,XY Y Y , v ,XY Y Y , u ,XXY Y , v ,XXY Y , u ,XXXY , v ,XXXY , u ,XXXX , v ,XXXX T (79)
The differential derivatives in F of E 3 elasticity with unknowns u = (u, v, w) T in material coordinates X = (X, Y, Z) are

∂U 3d = u ,Z , v ,Z ,
In NOM, the differential derivatives can be written as ∂U s = B s U , s ∈ {2d, 3d}, where B s is the operator matrix constructed with steps similar to that in Eq.71 and U is the vector form of all unknowns in one support. In other words, for a given equivalent energy functional, the independent derivatives of various orders can be extracted, which form a subset of the list of nonlocal derivatives provided by NOM. The operator matrix B s is formed by selecting a row with the same index of derivative in ∂α u i . The residual vector and tangent stiffness matrix for one particle can be obtained as

R i = ∂φ i ∂U = ∂(∂U ) ∂U ∂φ i ∂(∂U ) = B T ∂φ i ∂(∂U ) (81)
K i = ∂R i ∂U = B T ∂ 2 φ i ∂(∂U ) 2 B (82)
The explicit forms of ∂φ ∂(∂U ) and ∂ 2 φ ∂(∂U ) 2 can be obtained by softwares such as Mathematica [57] allowing symbolic operations. For simplicity, we omit these lengthy expressions in the paper.

However, the code will be made available. One can see that the construction of the residual vector and tangent stiffness matrix for each particle is a series of matrix multiplications.

The global tangent stiffness matrix for the functional in domain Ω can be expressed as

R Ω = ∆V i ∈Ω R i ∆V i , K Ω = ∆V i ∈Ω K i ∆V i (83)
The assembly of R i , K i is based on the global indices of all unknowns in one support. The global tangent stiffness matrices (e.g. N can be applied directly on the particles. The moment boundary condition on ∂Ω 1 N can be enforced by calculating the residual

K ∂Ω 0 D , K ∂Ω 1 D ) and residuals (e.g. R ∂Ω 0 D , R ∂Ω 1 D ) for functionals on boundaries ∂Ω 0 D , ∂Ω
R ∂Ω 1 N = ∆S i ∈∂Ω 1 N ∂φ N ∂(∂U N ) ∂(∂U N ) ∂U ∆S i = ∆S i ∈∂Ω 1 N B T N ∂φ N ∂(∂U ) ∆S i ( 84
)
where φ N = Q : ∇ n u and ∂U N = ∇ ⊗ u = B N U ; B N are constructed by selecting the first 2 rows or 3 rows of B in 2D or 3D, respectively; ∂φ N ∂(∂U N ) can be obtained by Mathematica [START_REF] Wolfram | The MATHEMATICA® book, version 4[END_REF]. Then, the global tangent stiffness matrix and residual are

R = R Ω + R ∂Ω 0 D + R ∂Ω 1 D + R ∂Ω 1 N (85) K = K Ω + K ∂Ω 0 D + K ∂Ω 1 D ( 86
)
With the global residual and tangent stiffness available, a standard Newton Raphson method can be used to find the solution.

Numerical examples

In this section, we present several representative numerical examples to study the property of the E n elasticity theory. The setup of the 2D/3D examples and the associated boundary conditions are outlined in Fig. 3. The domain is discretized with a Cartesian grid.

The particles in the void domain are removed to form holes. Each particle has the same number of neighbouring particles in the support, and the support size is selected as the distance between the furthest neighbour particle and the master particle in the support.

So the central particles have a smaller support size compared to the particles close to the boundary. The number of neighbours in support is selected as 1/2(n 2 + 33n + 32) in 2D,

where n is the order of gradient elasticity. The material parameters and length scales will be given in the subsections.

Convergence of strain energy in E 3 elasticity

The first example tests the strain energy distribution of E 3 elasticity for different dis- of the strain energy of different orders. The engineering strain is approximately 0.0288 as depicted in Fig. 4(a). The distribution of total strain energy density on each particle is given in Fig. 4(b). The maximal strain energy density occurs around the corners. The total strain energies on different strain gradient orders can be found in Fig. 5. With increasing the number of particles, the energies of different levels converge. The strain energy is dominant while the higher order energies tend to decrease with increasing gradient orders. Indeed, the deformation under pure tension load is "uniform" for this numerical example.

2D plate with uniform deformation

The second example tests the influence of E n gradient elasticity subjected to a uniform load; E 0 , E 1 , E 2 elasticity theories are implemented. The geometry and boundary conditions are illustrated in Fig. 3(a). A plate with dimensions of 1×1 m 2 is discretized into 81 2 particles.

The material parameters are elastic modulus E = 30 GPa and Poisson's ratio ν = 0.3. The internal length scales are l 1 = l 2 = 0.05. The left side of the plate is fixed in all directions and the right side is subjected to a uniform tension load of p = 1 MPa/m. Fig. 6(a) shows that the displacement based on higher order elasticity theory is identical to conventional Energy (J) (d)

Figure 5: Total strain energy on level (a)

1 2 S : E; (b) 1 2 l 2 1 ∇S:∇E; (c) 1 2 l 4 2 ∇ 2 S :: ∇ 2 E and (d) 1 2 l 6 3 ∇ 3 S • (5) ∇ 3 E; N is the total number of particles.
elasticity for uniform deformations since the higher order strain components are quite small such that their contribution to the energy density can be neglected. However, the higher order terms make the deformation smoother as shown in Fig. 6(b). This indicates that the higher order gradient elasticity should be tested with in-homogeneous deformations.

2D plate subjected to point force

Let us test the capability of gradient theory for point loads. We adopt the dimensions of the plate and its material parameters from the previous subsections. However, one particle in the middle of the right side boundary of the plate is subjected to a point force of 1000 N. The geometry and boundary conditions are depicted in Fig. 3(b). The plate is discretized into 81x81 particles. E 1 -E 4 elasticity theories are considered. The deformations of the plate for E 1 , E 2 , E 3 , E 4 elasticity can be found in Fig. 7. Obviously, gradient elasticity can 'withstand' point loads. The higher order gradient elasticity has a smoother displacement field compared with gradient elasticity. This observation is consistent with the numerical analysis by FEM [START_REF] Jc Reiher | Finite-element analysis of polyhedra under point and line forces in second-strain gradient elasticity[END_REF] and by IGA [START_REF] Makvandi | Isogeometric analysis of first and second strain gradient elasticity[END_REF]. Comparisons of the displacement in x-direction of particles on the right side boundary (i.e.the red line in Fig. 3(b)) are plotted in Fig. 8.

The first-order and second-order derivatives of the displacement in x-direction of particles are shown in Fig. 9. The derivative of the displacement in E 1 elasticity changes sharply, in contrast to the smooth transition of the displacement gradient in higher order gradient elasticity.

Figure 9: The first-order and second-order derivative of deformation of all particles on the red line in Fig. 3(b).

displacement gradient in x-direction of E 2 elasticity can be found in Fig. 12. The larger l i , the smaller the strain field.

Ux (m)

Large deformation of 2D plate with a hole

Again, we adopt the material parameters and plate dimensions from the previous example and study the deformation of a 2D plate with a hole based on E n elasticity, with The displacements at step T = 1 are plotted in Fig. 13, where E 0 elasticity has the largest deformation and the deformations by E 2 , E 3 , E 4 , E 5 elasticity are similar. Fig. 14 shows that the displacement in y-direction of particle on the bottom line (e.g. the blue line in Fig. 3(c)). It can be seen that the higher order gradient theory has smaller deformation. The difference becomes smaller when the order of gradient elasticity increases. Contour plots of displacement gradients for E 0 , ..., E 5 elasticity are plotted in Fig. 15 and Fig. 16. Higher order elasticity exhibits a very smooth displacement gradient. The gradient of the displacement field for hyperelasticity (E 0 elasticity) is not smooth around the internal line. This is due to the fact that the first order NOM is used, which is continuous in the displacement but discontinuous in its derivative. Although the deformations for different elasticity theories at T = 1 are similar, the final

u ∂Z 2 , (b) ∂ 2 u ∂Y 2 , (c) ∂ 2 u ∂X 2 , (d) ∂ 2 v ∂Z 2 , (e) ∂ 2 v ∂Y 2 ,(f) ∂ 2 v ∂X 2 .

Conclusions

We have proposed an objective energy functional for finite deformation higher order gradient elasticity. The energy functional is based on the setting of the second Piola-Kirchhoff stress which is invariant under rigid body transformations. More specifically, the geometric nonlinear higher order gradient elasticity theory is formulated on the gradients of the right Cauchy Green tensor. The general form of higher order gradient elasticity may contain thousands of material parameters and we proposed a simplified version of gradient elasticity. Such simplification reduces the number of material parameters from 10 thousands to less than 10. A small number of material parameters can greatly simplify the experiment measurement and numerical implementation. The framework of gradient elasticity also allows for other forms of simplification of material parameters. We employed the nonlocal local operator method and Newton Raphson iteration method to find the numerical solution of higher gradient elasticity. The properties of gradient elasticity are studied by a series of numerical experiments. The numerical tests show that gradient elasticity can sustain point/line load without stress singularity. The mechanical response greatly depends on the internal length scales of gradient elasticity. Larger internal length scale induces a smaller and smoother deformation. Higher order gradient elasticity is numerically more stable and allows for larger ultimate load for the same structure. In the next stage, more physics-related research including the calibration of material parameters by experiments and numerical simulation, and the size effect, surface effect in metamaterials and gradient elasticity will be pursued. Some outlooks based on current research include, for example, 1. The higher order gradient elastoplasticity theory [START_REF] Reiher | Finite third-order gradient elastoplasticity and thermoplasticity[END_REF] and its numerical implementation.

Current research is restricted to elasticity with finite deformation and it cannot be applied to a dissipated system involving permanent deformation or irreversible process. The extension of higher order elasticity to higher order plasticity can broaden the range of plasticity theory.

2. More clear relationship between metamaterial and gradient elasticity is expected [START_REF] Khakalo | Modelling size-dependent bending, buckling and vibrations of 2d triangular lattices by strain gradient elasticity models: applications to sandwich beams and auxetics[END_REF][START_REF] Yang | Effective strain gradient continuum model of metamaterials and size effects analysis[END_REF].

One salient feature of gradient elasticity is the micro-structure, which is essential to the theory of metamaterials as well. Direct simulation of micro-structure requires tremendous computer power. Gradient elasticity circumvents these difficulties by introducing certain internal length scales, which however require more sophisticated measurement.

3. The wave propagation analysis of gradient elasticity. Gradient elasticity has the capability to account for interesting phenomena such as size effect, surface effect and nonlocal effect. These features may give rise to some abnormal physical phenomenon, which can be exploited to design some smart devices.

Appendix A. Symmetry of higher order tensor

For the 4th-order elasticity tensor, the symmetry can significantly reduce the number of material parameters. The symmetry of the Cauchy stress tensor (σ ij = σ ji and the generalized Hooke's laws (σ ij = C ijkl ε kl ) implies that C ijkl = C jikl . Similarly, the symmetry of the infinitesimal strain tensor implies that C ijkl = C ijlk . These symmetries are called the minor symmetries If in addition, since the displacement gradient and the Cauchy stress are work conjugate, the stress-strain relation can be derived from a strain energy density functional (U ), then

σ ij = ∂U ∂ε ij =⇒ C ijkl = ∂ 2 U ∂ε ij ∂ε kl . (A.1)
The arbitrariness of the order of differentiation implies that C ijkl = C klij .

The stiffness matrix C satisfies a given symmetry condition if it does not change when subjected to the corresponding orthogonal transformation, which may represent symmetry with respect to a point, an axis, or a plane.

According to [START_REF] Forte | Symmetry classes for elasticity tensors[END_REF][START_REF] Geymonat | Classes de symétrie des solides piézoélectriques[END_REF], the orthogonal transformation of a tensor of any order can be written as The general rotation matrix R can be written as

T (M ) := (Q M ) ...ijk... = • • • Q ip Q jq Q kr • • • M ...
R = R z (α) R y (β) R x (γ)
The coordinate transform of a vector in matrix and tensor notation is

v = Q • v and v i = λ ij v j (A.7)
The coordinate transform of a tensor in matrix and tensor notation is

σ = Q • σ • Q T and σ mn = λ mi λ nj σ ij (A.8)
The coordinate transform of a 4th-order tensor is Based on the symmetry calculation in Appendix A, the matrix form of isotropic gradient elasticity can be derived accordingly.

C = Q • Q • C • Q T • Q T , C ijkl =

Figure 1 :

 1 Figure 1: Domain with piecewise smooth surfaces; n is the outward unit normal direction of surface boundary ∂Ω, γ is the tangent direction of line boundary ∂ 2 Ω ij , m = γ × n is the outward normal unit vector of line boundary.

  1) indicate the Dirichlet and Neumann boundary conditions for u of different partial derivative orders.

Figure 2 :

 2 Figure 2: (a) Domain and notation. (b) Schematic diagram for support and dual-support, all shapes above are supports,S X = {X 1 , X 2 , X 4 }, S X = {X 1 , X 2 , X 3 }.

Figure 3 :

 3 Figure 3: The setup of the 2D plate and 3D plate and boundary conditions.

Figure 4 :

 4 Figure 4: (a) Deformation with unit of meters in x-direction, scaled by 10 times and (b) the distribution of total strain energy density with units of Joule per unit volume for discretization of 120 2 particles.

Figure 6 :

 6 Figure 6: (a) Displacement in x-direction and (b) strain in x-direction for particles on middle horizontal line (Blue line in Fig.3(a)).

Figure 7 :Figure 8 :

 78 Figure 7: Displacement with unit of meters in x-direction of the plate with deformation scaled by 10 7 times for (a) E 1 elasticity, (b) E 2 elasticity, (c) E 3 elasticity and (d) E 4 elasticity.

Figure 10 :Figure 11 :

 1011 Figure 10: The displacement in x-direction of particles on right side boundary of the plate based on E 1 elasticity and (b) E 2 elasticity.

Figure 12 :

 12 Figure 12: ∂u ∂X distribution based on E 2 elasticity for (a) l 1 = l 2 = 0.05, (b) l 1 = l 2 = 0.1, (c) l 1 = l 2 = 0.025.

Figure 13 :

 13 Figure 13: The deformations of the plate at T = 1 for (a) E 0 elasticity,(b) E 1 elasticity,(c) E 2 elasticity,(d) E 3 elasticity,(e) E 4 elasticity and (f) E 5 elasticity, respectively. Unit:meters.

Figure 14 :

 14 Figure 14: The displacement in y-direction of all particles on the right side boundary of the plate with load level T = 1, where the lines in (b) are magnified from (a).

Figure 15 :

 15 Figure 15: ∂u ∂Y at T = 1 for E 0 , ..., E 5 elasticity.

Figure 16

 16 Figure 16: ∂v ∂Y at T = 1 for E 0 , ..., E 5 elasticity.

18 5. 6 .Figure 17 :Figure 18 :Figure 19 : 35 -Figure 20 :Figure 21 :

 186171819352021 Figure 17: The converged final deformations of plate for (a) E 0 elasticity,(b) E 1 elasticity,(c) E 2 elasticity,(d) E 3 elasticity,(e) E 4 elasticity and (f) E 5 elasticity, respectively.

  pqr...(A.2) where Q is an orthogonal matrix given byO(n, R) = Q ∈ GL(n, R)|Q T Q = QQ T = I , (A.3)with GL(n, R) being the set of all real n × n matrices and I the identity matrix.The symmetry of certain orthogonal transformation of a tensor requiresQ M = M , with Q ∈ O(n, R). (A.4)The common orthogonal transformations in 3D include the reflection, rotation. The transformation matrices for symmetry planes are

  w ,Z , ..., v ,XY , ..., u ,XY Z , ..., w ,XXXX

	T
	102 terms

  1 D can be obtained in the same manner. The Neumann boundary condition on ∂Ω 0

  λ im λ jn λ ko λ lp C mnop (A.9) -→ ∇ε = (ε 111 , ε 122 , ε 133 , 2ε 123 , 2ε 113 , 2ε 112 , ε 211 , ε 222 , ε 233 , 2ε 223 , 2ε 213 , 2ε 212 , ε 311 , ε 322 , ε 333 , 2ε 323 , 2ε 313 , 2ε 312 ) (B.5) are the Voigt notations of jk and mn, respectively. Then the vectorial forms of couple stress and strain gradient can be written as -→ ∇ε = ( 11 , 12 , 13 , 2 14 , 2 15 , 2 16 , 21 , 22 , 23 , 2 24 , 2 25 , 2 26 , 31 , 32 , 33 , 2 34 , 2 35 , 2 36 ) (B.9) --→ ∇σ = (σ 11 , σ 12 , σ 13 , σ 14 , σ 15 , σ 16 , σ 21 , σ 22 , σ 23 , σ 24 , σ 25 , σ 26 , σ 31 , σ 32 , σ 33 , σ 34 , σ 35 , σ 36 ) (B.10)

	Based on Voigt rotations
	ij = 11 22 33 23, 32 13, 31 12, 21
	⇓	⇓ ⇓ ⇓	⇓	⇓	⇓	(B.6)
	α =	1 2 3	4	5	6
	we write the couple stress and strain gradient as
		iα =	∂ε α ∂x i	, σ iα =	∂σ α ∂x i		(B.7)
			h ijklmn → h iαlβ		(B.8)
	where α, β

The coordinate transform of a 6th-order tensor is

Solving Eqs .A.9,A.10 by Mathematica [START_REF] Wolfram | The MATHEMATICA® book, version 4[END_REF], we can obtain the independent variables in high order tensor.

There are 3 6 = 729 terms in H. The Minor symmetry reduces H into 171 independent terms.

The orthotropy requires H = H, C = C for three reflection symmetries A 1 , A 2 , A 3 .

The case of orthotropy (the symmetry of a brick) has 51 independent elements.

The isotropy property requires H = H, C = C for any rotation. This requirement reduces the number of independent terms in H from 171 to 5.

Appendix B. Matrix Form of strain gradient energy by Viogt notations

The tensor form of higher order tensor contains many repeated terms when symmetry property is considered. In terms of numerical implementation, it is more convenient to use the matrix form than to use tensor form. In conventional mechanics, the Voigt notation is an efficient method to formulate the matrix form. Let us take the strain gradient linear elasticity for an example. The other higher order tensor can be formulated in the same manner. The material constitutive for couple-stresses can be written as

The strain-gradient energy function is

where ε ijk , σ ijk are defined as

The vectorial forms of couple stress and strain gradient can be written as


Isabella Silva Barreto
email: isabella.silva_barreto@bme.lth.se

M Pierantoni

M Hammerman

E Törnquist

S Le Cann

A Diaz

J Engqvist

M Liebi

P Eliasson

H Isaksson

Characterization of collagen structural response to in situ loading of the rat Achilles tendon

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The mechanical properties of tendons are depending on a complex hierarchical design, where collagen is arranged into fibrils in a periodic, quarter staggered pattern at an intermolecular distance referred to as d-spacing. This study aimed to characterize the Achilles tendon's mechanical response to in situ loading at the nanoscale and how it relates to the resulting macroscopic tissue mechanical behavior.

Methods

Achilles tendons (n=16) were harvested from Sprague Dawley rats (10-14 weeks) and stored frozen in NaCl until measurements at the SAXS beamline (cSAXS, Swiss Light Source, PSI) (energy 12.4keV, beam size 150μm). The tendons were loaded using a custom-made uniaxial tensile device (Fig 1) [1] while SAXS patterns (q=0.05-1.45nm -1 ) were acquired simultaneously from the center of the tendon. A radiation damage test was conducted prior to testing which ensured that all analysis was conducted below 28kGy. The tendons were loaded (5mm/min) in ramp to failure (N=5), cyclic loading (N=5) or stress relaxation (N=6). From collagen peaks identified from the azimuthally integrated scattering intensity over the angular sector θ (Fig 1 .A), local fibril strain and collagen structural parameters were extracted [2,3]. All analysis was performed using MATLAB.

Results and Discussion

For all loading schemes, the fibril strain, strain heterogeneity and overlap followed the macroscopic load. During ramp to failure, the fibril strain increased linearly with tissue strain in the elastic region (

  Fig 2.A). The fibrils in the beam path yielded and broke (fibril strains 1.3±0.3%) prior to global tissue failure, which could be the result of an accumulation of nanoscale ruptures throughout the specimens. With increasing load, the fibril d-spacing increased and the overlap region decreased, which indicates fibril stretching and sliding simultaneously. During cyclic loading, the tissue hysteresis decreased with number of cycles (Fig 2.B, i). Interestingly, while the d-spacing at maximum load remained constant for all cycles (0.9±0.5%), it decreased sequentially at the end of each cycle (Fig 2.B, ii). This suggests that the initial decline in tissue hysteresis is accompanied by an increase in fibril recovery and stretchability. During stress relaxation, most of the collagen parameters experienced transient relaxation simultaneously as the tissue (Fig 2.C).ConclusionsDifferent loading scenarios enabled linking of the nano and macroscale in the Achilles tendon in terms of mechanics, such as accumulation of local fibril ruptures before global failure and adaptation of the nanostructure to the applied load.

Figure 1 .

 1 Figure 1. A) schematic of experiment and B) the applied loading schemes.

Figure 2 .

 2 Figure 2. Changes in fibril parameters during A) ramp to failure, B) cyclic loading and C) stress relaxation.

Acknowledgments

This project was funded by the European Union's Horizon 2020, No 731019 (EUSMI) through beamtime at cSAXS beamline, PSI, and by the Knut and Alice Wallenberg Foundation.