domain: spi math

text:
Hussein Nassar
email: nassarh@missouri.edu

Arthur Lebée
email: arthur.lebee@enpc.fr

Emily Werner

Strain compatibility and gradient elasticity in morphing origami metamaterials

Keywords: Origami, Compliant shells, Poisson coefficient, Gradient elasticity, Geometric nonlinearity

The principles of origami design have proven useful in a number of technological applications. Origami tessellations in particular constitute a class of morphing metamaterials with unusual geometric and elastic properties. Although inextensible in principle, fine creases allow origami metamaterials to effectively deform non-isometrically. Determining the strains that are compatible with coarse-grained origami kinematics as well as the corresponding elasticity functionals is paramount to understanding and controlling the morphing paths of origami metamaterials. Here, within a unified theory, we solve this problem for a wide array of well-known origami tessellations including the Miura-ori as well as its more formidable oblique, non-developable and non-flat-foldable variants. We find that these patterns exhibit two universal properties. On one hand, they all admit equal but opposite in-plane and out-of-plane Poisson's ratios. On the other hand, their bending energy detaches from their in-plane strain and depends instead on the strain gradient. The results are illustrated over a case study of the self-equilibrium geometry of origami pillars.

Introduction

Origami has long surpassed its artistic vocation. Its early uses in architectural geometry date back to the 1920s [START_REF] Wingler | Bauhaus[END_REF][START_REF] Lebée | From folds to structures, a review[END_REF]. By now, origami techniques have been applied to the design of deployable structures in astronautics, robotics and biomedical engineering [START_REF] Kuribayashi | Self-deployable origami stent grafts as a biomedical application of Ni-rich TiNi shape memory alloy foil[END_REF][START_REF] Schenk | Folded Shell Structures[END_REF][START_REF] Zirbel | Accommodating thickness in Origami-based deployable arrays[END_REF][START_REF] Callens | From flat sheets to curved geometries: Origami and kirigami approaches[END_REF]. Origami tessellations in particular make up versatile morphing 2D metamaterials that can access various low-energy high-curvature configurations [START_REF] Seffen | Compliant shell mechanisms[END_REF][START_REF] Norman | Multistable and morphing corrugated shell structures[END_REF]. Although inextensible in principle, the fine crease patterns of origami metamaterials enable them to effectively deform non-isometrically. In this letter, we determine the strains that are compatible with coarse-grained origami kinematics. Strain compatibility permits to fully characterize the 3D-embedded geometry when suitable boundary conditions are prescribed. In most cases however, the geometry further solves an elastostatic equilibrium for which we calculate the appropriate elasticity functional.

Morphing in origami metamaterials can be triggered in response to different loading conditions such as bending, twisting and stretching all of which couple in a manner unparalleled in other plates and shells; see Fig. 1. For instance, standard pure bending kinematics of plates imply that the curvatures in the bending plane and in the plane orthogonal to it, κ 1 and κ 2 , are in a proportion equal and opposite to the in-plane Poisson's ratio ν, i.e., ν = -κ 2 /κ 1 . Thus auxetic plates tend to bend synclastically while anauxetic plates tend to bend anticlastically [START_REF] Naboni | Embedding auxetic properties in designing active-bending gridshells[END_REF]. Origami metamaterials do not [START_REF] Schenk | Origami folding: A structural engineering approach[END_REF]. In fact, in what follows, we demonstrate that they typically exhibit the exact opposite behavior, i.e., ν = +κ 2 /κ 1 . Other loads that induce morphing, such as "pinching", are unique to origami [START_REF] Grey | Strain Reversal in Actuated Origami Structures[END_REF]. Indeed, in plates and membranes, a pinch usually produces highly localized and singular deformations. By contrast, pinching an origami metamaterial brings on global and nearly-uniform deformations [START_REF] Seppecher | Linear elastic trusses leading to continua with exotic mechanical interactions[END_REF][START_REF] Filipov | Origami tubes assembled into stiff, yet reconfigurable structures and metamaterials[END_REF]. The consequence, as we shall prove, is that membrane energy penalizes, not strains, but their gradients, i.e., the Christoffel symbols [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Gray | Modern differential geometry of curves and surfaces with Mathematica, 3rd Edition[END_REF][START_REF] Ciarlet | An introduction to differential geometry with applications to elasticity[END_REF].

Figure 1: Two views of an origami-folded hyperboloid: A periodic "Mars" pattern is laser-etched into a PET sheet then hand-folded, circularly bent and fixated with adhesive tape. The pattern is auxetic but bends into a saddle. The letter explains the observed self-equilibrium geometry.

Previous contributions have succeeded in describing several aspects of origami mechanics, including calculations of Poisson's ratios and bending moduli [START_REF] Wei | Geometric mechanics of periodic pleated origami[END_REF][START_REF] Schenk | Geometry of Miura-folded metamaterials[END_REF][START_REF] Pratapa | Geometric mechanics of origami patterns exhibiting Poisson's ratio switch by breaking mountain and valley assignment[END_REF]. Hereafter, we derive an effective medium theory complete with its field equations of compatibility and equilibrium. The theory naturally accounts for the observations described above and shows how the different deformation modes couple kinematically and energetically. The subject of the theory is a class of origami patterns that we call "zigzag sums". These are generated by translating one zigzag line along another so as to sweep a periodic polyhedral surface Σ whose edges and panels we interpret as an origami crease pattern; see Fig. 2a. Equivalently, Σ is a periodic crease pattern with a unit cell composed of four parallelograms. Four quadrilaterals per unit cell is as simple as origami tessellations can be without being trivial; moreover, here, the quads are parallelograms. Nonetheless, zigzag sums span various well-known patterns (e.g., the Miura-ori, the eggbox, "Sakoda's staircase" and "Barreto's Mars" patterns [START_REF] Lang | Twists, Tilings, and Tessellations[END_REF]) and encompass their oblique, non-developable and non-flat-foldable variants by maintaining two basic requirements: periodicity and parallelism. Meanwhile, the zigzag generators provide a handful of geometric seeds that can readily tune their deformation paths as well as the underlying elastic behavior.

Geometric mechanics of zigzag sums

The deformations of a zigzag sum Σ are the outcome of two competing mechanisms [START_REF] Lechenault | Mechanical response of a creased sheet[END_REF]. The first is the folding and unfolding of Σ along its edges. The second is the isometric bending of its panels. Signature origami kinematics emerges when the bending stiffness of the panels is significantly larger than the folding stiffness of the edges. Accordingly, we can assume that individual panels bend infinitesimally; Σ, on the other hand, can still embrace finitely curved surfaces by accumulating small deflections across many unit cells. Formally, we characterize the deformations of Σ in terms of a continuous parametrization x = x(ξ) of its mid-plane, where x is the current position of vertex ξ, with ξ = (ξ 1 , ξ 2 ) being the effective curvilinear coordinate along the zigzags directions; see Fig. 2b-e. Consider then a nonelementary unit cell, i.e., a supercell, centered about vertex ξ and containing the neighboring vertices ξ (ij) placed at x (ij) . Then, the adopted origami kinematics imply that position

x (ij) = X (ij) + r 2 δx (ij) (1)
is attained by a fold-only motion ξ (ij) → X (ij) (Fig. 2b-c) followed by a bending-induced perturbation r2 δx (ij) that is small compared to the characteristic size of the panels denoted r (Fig. 2c-d). In the limit of a tessellation where Σ contains a large number of unit cells, i.e., for r → 0, the plane locally tangent to Σ is spanned by two vectors

x 1 ≡ ∂x ∂ξ 1 = lim r→0 x (10) -x ( 10) r = u -v, x 2 ≡ ∂x ∂ξ 2 = lim r→0 x (01) -x (0 1) r = s -t, (2)
where u, v, s and t are the crease vectors meeting at vertex ξ and normalized with respect to r. The apparent in-plane deformations of Σ can therefore be quantified in terms of an effective metric tensor g of components 7) decreases like r or, equivalently, like the reciprocal of the number of unit cells per row. The plot also shows that the maximum angular deflection rδ (ij) decreases in the same fashion confirming the fold-to-bend orders of magnitude adopted in (1) and that, in the limit r → 0, all panels are flat. The inset shows the unit cell adopted for this example. (b,c) Two iterates: snapshots of the way in which the zigzag sum converges to a smooth surface taken for 64 and 128 unit cells per row. Note that the zigzag sum is finitely and non-uniformly stretched and bent. Illustrated surfaces are constructed by "propagating" inextensibility constraints away from a given boundary; see [START_REF] Nassar | Fitting surfaces with the Miura tessellation[END_REF] for details.

g µν = ⟨x µ , x ν ⟩ = ∥u -v∥
It is noteworthy that g 12 = ⟨u, s⟩ + ⟨v, t⟩ -⟨u, t⟩ -⟨v, s⟩ is a motion constant (i.e., dg 12 ≡ 0) as it is equal to a combination of constant lengths and angles. Thus, if g 12 = 0 so that the zigzags ξ 1 -and ξ 2 -contours are initially orthogonal, then orthogonality is maintained in any subsequent motion (Fig. 2f). However, in oblique tessellations, with g 12 = ∥x 1 ∥ ∥x 2 ∥ cos θ ̸ = 0, stretching couples to shearing according to

tan θ dθ = d∥x 1 ∥ ∥x 1 ∥ + d∥x 2 ∥ ∥x 2 ∥ = dg 2g 11 g 22 (4)
with g = det g. Rather intuitively then, shearing, understood as aligning the zigzags (θ → 0 or π), always goes to reduce the pattern's effective area (Fig 2g Most important is the observation that the fold-only motion ξ (ij) → X (ij) maintains crease parallelism and periodicity: this is truly where the fact that Σ is made out of parallelograms is crucial. Hence, folding alone produces no curvature and, in that regard, it is bending that takes over now. This is most apparent in the expressions of the second-order derivatives which, in the limit r → 0, read (20) -2δx + δx ( 20) ,

x 11 = δx
x 22 = δx (02) -2δx + δx (0 2) ,
x 12 = δx (11) -δx ( 11) -δx (1 1) + δx ( 11 ) , [START_REF] Zirbel | Accommodating thickness in Origami-based deployable arrays[END_REF] where x µν ≡ ∂ 2 x/∂ξ µ ∂ξ ν . Therefore, the x µν at ξ are linear functions of the bending DOFs of the super cell centered at ξ. That being said, not all δx (ij) are compatible with the inextensibility constraints. Ultimately, it is possible to show that (x 11 , x 22 , x 12 ) belong to a linear subspace of R 3 × R

x 22 = ⟨u, t ∧ s⟩ ⟨v, t ∧ s⟩ ⟨t, u ∧ v⟩ ⟨s, u ∧ v⟩ x 11 . (6)
The tangent components of x 22 and x 11 describe the strain gradients in the mid-plane of Σ in direction x 1 and x 2 whereas the normal components describe the normal curvatures κ 1 and κ 2 of the mid-plane in the same directions.

The above identity states that these respective quantities necessarily occur in equal proportions. Specifically, the in-plane and out-of-plane Poisson's ratios, ν and -κ 2 /κ 1 , of a zigzag sum Σ are equal and opposite and are given by the current geometry of the folds as per

ν = κ 2 κ 1 = g 11 g 22 ⟨u, t ∧ s⟩ ⟨v, t ∧ s⟩ ⟨t, u ∧ v⟩ ⟨s, u ∧ v⟩ , (7)
with

κ µ = ⟨x µµ , n⟩ /g µµ (no sum), n = x 1 ∧ x 2 / ∥x 1 ∧ x 2 ∥ being the unit normal [23, App. C].
Identity [START_REF] Seffen | Compliant shell mechanisms[END_REF] is of an asymptotic nature: in the limit r → 0 of infinitesimal unit cells, the Poisson's ratios of a zigzag sum Σ converge to one another. In other, more practical, terms, the equality is accurate within an error proportional to r. See Fig 3 for a numerical illustration of the convergence process and the relevant error bounds. Identity [START_REF] Seffen | Compliant shell mechanisms[END_REF] is also universal: it holds whether Σ is developable, flat-foldable, symmetric, rectangular, or not. Therein, the mixed products (e.g., ⟨u, t ∧ s⟩) have their signs determined by the mountain-valley assignment of the concerned folds. This allows us to deduce the rather nice corollary: if a zigzag sum Σ has an odd number of mountain folds meeting at a vertex then: (i) it is necessarily auxetic; and (ii) it necessarily bends into a saddle. In particular, all developable zigzag sums are auxetic and bend into saddles. Note that these results are only valid for nondegenerate configurations, i.e., where all folds are partially folded. Indeed, when a fold is flat, be it open or closed, mountain and valley assignments coalesce and the Poisson's ratio could become null, infinite or multivalued [START_REF] Pratapa | Geometric mechanics of origami patterns exhibiting Poisson's ratio switch by breaking mountain and valley assignment[END_REF]. Note also that the proven identity features the elongations and curvatures in two specific material directions ξ 1 -and ξ 2 -contours, be them orthogonal or not; in particular, the identity bears no immediate consequences on the Poisson's coefficient defined for two general orthogonal directions.

Self-equilibrium of origami pillars

The equality of the Poisson's ratios greatly constrains the 3D geometries accessible to a zigzag sum Σ [START_REF] Nassar | Fitting surfaces with the Miura tessellation[END_REF][START_REF] Nassar | Curvature, metric and parametrization of origami tessellations: theory and application to the eggbox pattern[END_REF][START_REF] Nassar | Macroscopic deformation modes of origami tessellations and periodic pin-jointed trusses: the case of the eggbox[END_REF][START_REF] Zheng | A continuum field theory for the deformations of planar kirigami[END_REF]. Consider, for tractability, the "Mars" tessellation of Fig. 4a. The tessellation lacks mirror symmetry but is otherwise developable, flat-foldable and equilateral. We find that compatible metrics and curvatures satisfy [23, App. D]

g 12 = 0, g 11 (4 -g 22 ) = 4(cos α + cos β) 2 , κ 2 κ 1 = ν = -4 (cos α + cos β) 2 g 11 g 22 . (8)
This is in fact a system of non-linear PDEs weighing on the admissible configurations x. For instance, we explore the configurations of a rectangular domain of flat dimensions L 1 × L 2 that are folded and wrapped around an axis of symmetry into an origami pillar [START_REF] Lang | Twists, Tilings, and Tessellations[END_REF]. These pillars are surfaces of revolution with a Cartesian parametrization of the form

x(ξ) = (ρ(ξ 2 ) cos(ωξ 1 ), ρ(ξ 2 ) sin(ωξ 1 ), z(ξ 2 )). (9)
Here, the axis of symmetry is the z-axis, ξ 1 -contours are meridian lines and ξ 2 -contours are parallel lines. Equations [START_REF] Norman | Multistable and morphing corrugated shell structures[END_REF] then simplify into a system of ODEs governing z and ρ. Its solutions constitute a 2-DOF family of onesheeted hyperboloids (Fig. 4b,c) [23, App. E]. The range of motion of each pillar is bounded by the maximally folded and unfolded states of the pattern. Thus, for increasing aspect ratio L 2 /L 1 , the pillar will have a smaller range of motion, until it jams, or even becomes impossible to form without tearing.

Among all accessible pillars, the ones that are in selfequilibrium exhibit minimum levels of strain energy ψ. Strain energy is composed of two contributions: crease folding and panel bending. Crease folding energy takes the form of a classical, however nonlinear, membrane energy. In what follows, we neglect this contribution and focus on the less-explored influence of panel bending [START_REF] Grey | Mechanics of Paper-Folded Origami : A Cautionary Tale[END_REF]. For "Mars"-folded origami pillars, centered about the equator, strain energy is exemplified on Fig. 4d. It appears that unfolding the pattern, by increasing the angle ∠(u, v), reduces bending energy. The unfolding flattens the pillar and reduces both curvatures κ 1 and κ 2 . That trend continues up to a point where energy is minimum (near state (iv)) and beyond which any further unfolding, or flattening, of the pattern actually causes a dramatic increase in bending energy. This suggests that bending energy ψ depends on some other deformation measures besides curvatures. As a matter of fact, the bending energy density b, where ψ = Ω bdξ, Ω being the reference domain of the pattern, is a quadratic form of the planarity defects δ (ij) of the panels. For their part, the defects are linear forms of the full, both out-of-plane and in-plane, components of the parametrization's second derivatives x µν . The out-of-plane components Γ 3µν ≡ ⟨n, x µν ⟩ are the coefficients of the second fundamental form; they quantify the curvatures and torsion of the embraced surface. The in-plane components Γ σµν ≡ ⟨x σ , x µν ⟩ are the Christoffel symbols; they quantify, not the strains, but their gradient, i.e., terms of the form g µν,σ . Accordingly, b = b(Γ; g) is a metric-dependent quadratic form of the curvatures, torsion, and strain gradient. Detailed derivations of the expression of b leading to Fig. 4d Back to the pillars, unfolding the pattern reduces both curvatures and strain-gradients. However, as the geometric rigidities associated with the curvatures (i.e., ∂ 2 b/∂Γ 3µν ) remain bounded, the ones associated with the strain gradient (i.e., ∂ 2 b/∂Γ σµν ) diverge for states that are close to being flat. Such states occur near the outer rims of sufficiently unfolded pillars; their presence further signals that the pillar has reached the boundary of the kinematically admissible domain, hence the energy blow-up observed on Fig. 4d. Strain-gradient energy further dominates the response of any plane, non-uniform, state such as the "ring" observed midway through Fig. 4b. Indeed, uniformly folded plane states do not engage panel bending and therefore have zero energy. By contrast, non-uniformly folded states, with gradients of folding angles, cannot be achieved without panel bending. When the state is plane, it has zero curvatures and torsion, and energy becomes function of the strain gradient alone [START_REF] Seppecher | Linear elastic trusses leading to continua with exotic mechanical interactions[END_REF]. More generally, it is noteworthy that ψ does not penalize strains or fold angles, however large, so long as they are uniform. In fact, ψ does not even refer to a specific natural state or any specific natural fold angles in reference to which strains should be measured [START_REF] Ciarlet | An introduction to differential geometry with applications to elasticity[END_REF]. Instead, it refers to a higher-order strain measure, namely Γ, which quantifies not how much the pattern is folded, but rather how far it is from being uniformly folded.

Conclusion

In conclusion, we proposed a continuum theory of origami-folded shells where the crease pattern is a zigzag sum. The theory accounts for geometric nonlinearities and accurately predicts the states that are accessible under isometric folding thanks to a universal identity relating inplane deformations and out-of-plane curvatures. The study further demonstrates that the energy density of origamifolded shells can heavily depend on the in-plane straingradient in a way that is unparalleled in classical shell theories. It should be possible to extend the proposed theory, in one form or another, to more general tessellations that possess periodic folding motions. By contrast, alternative approaches are likely to be needed to tackle tessellations that do not fold periodically (e.g., Huffman grids and Yoshimura-like patterns [START_REF] Evans | Rigidly foldable origami gadgets and tessellations[END_REF]). Finally, it would be of interest to investigate how the proposed theory can inform shape programming and morphing planning of origami and kirigami structures [START_REF] Tachi | Designing freeform origami tessellations by generalizing Resch's patterns[END_REF][START_REF] Dudte | Programming curvature using origami tessellations[END_REF][START_REF] Guseinov | Programming temporal morphing of self-actuated shells[END_REF][START_REF] Siéfert | Programming stiff inflatable shells from planar patterned fabrics[END_REF][START_REF] Jin | Kirigami-Inspired Inflatables with Programmable Shapes[END_REF][START_REF] Feng | The designs and deformations of rigidly and flat-foldable quadrilateral mesh origami[END_REF][START_REF] Li | Metamorphosis of three-dimensional kirigamiinspired reconfigurable and reprogrammable architected matter[END_REF][START_REF] Yang | Grasping with kirigami shells[END_REF][START_REF] Melancon | Multistable inflatable origami structures at the metre scale[END_REF] for robotics and 4D printing applications. motion and is perturbed by r 2 δx subsequent to the infinitesimal bending of the panels. Note that the second-order derivatives x µν are invariant by composition with linear motions including rigid body motions and periodic stretching and contraction. Accordingly, it is possible, without loss of generality, to default the displacements of the central four creases to zero. In other words, we set δx = δx (10) = δx (01) = δx ( 10) = δx (0 1) = 0.

(B.1) Furthermore, by linearity, the contributions of the remaining displacements can be investigated independently then superposed. So let δx ( 11) , δx ( 11 ) , δx (1 1) , δx ( 20) and δx (0 2) all be null for now and consider a nonzero displacement δx (11) : the latter must be orthogonal to u so as to preserve the length of x (11) x (01) as well as to s so as to preserve the length of x (11) x (10) , to leading order. Hence, δx (11) = δ (11)

u ∧ s (B.2)
for some planarity defect δ (11) . Similarly, we have

δx (02) = δ + v ∧ t, δx (20) = δ -v ∧ t. (B.3)
Expressions for δ ± are easily obtained by observing that the lengths of x (11) x (02) and x (11) x (20) are preserved. Indeed, infinitesimal inextensibility implies t + u, δ (11)

u ∧ s -δ + v ∧ t = 0, v + s, δ (11) u ∧ s -δ -v ∧ t = 0, (B.4) namely, δ + δ -= ⟨u ∧ s, t⟩ ⟨v ∧ t, s⟩ ⟨v ∧ t, u⟩ ⟨u ∧ s, v⟩ . (B.5)
As for the x µν , they are given by

x 11 = δx (20) = δ -v ∧ t, x 22 = δx (02) = δ + v ∧ t,
x 12 = δx (11) = δ (11) u ∧ s.

(B.6)
In particular,

x 22 = ⟨u ∧ s, t⟩ ⟨v ∧ t, s⟩ ⟨v ∧ t, u⟩ ⟨u ∧ s, v⟩ x 11 . (B.7)
The contributions of δx ( 11) , δx ( 11 ) and δx (1 1) (i.e., δ ( 11) , δ ( 11 ) and δ (1 1) , respectively) can be readily calculated in the same fashion and shown to yield exactly the same constraint on x 22 and x 11 . This can be predicted by appreciating the symmetric way in which the proportionality coefficient depends upon the crease vectors. In conclusion, the above constraint holds for any combination of admissible infinitesimal displacements of the super cell.

Last, note that ⟨x 12 , x 1 ⟩ = δ (11) u ∧ s, uv = -δ (11) ⟨u ∧ s, v⟩

= -δ -⟨v ∧ t, s⟩ = -δ -⟨v ∧ t, s -t⟩ = -⟨x 11 , x 2 ⟩ . (B.8)
That is: ∂ 1 ⟨x 1 , x 2 ⟩ = 0. In the same manner, one shows that ∂ 2 ⟨x 1 , x 2 ⟩ = 0 in order to recover dg 12 = 0 as claimed in the main text.

Appendix C. The in-and out-of-plane Poisson's coefficients are equal and opposite

Projecting the proportionality constraint of x 11 and x 22 along the normal n yields

κ 2 g 22 = ⟨u ∧ s, t⟩ ⟨v ∧ t, s⟩ ⟨v ∧ t, u⟩ ⟨u ∧ s, v⟩ κ 1 g 11 (C.1)
by the definition of the normal curvatures in directions 1 and 2. Projecting along x 1 yields

⟨x 22 , x 1 ⟩ = ⟨u ∧ s, t⟩ ⟨v ∧ t, s⟩ ⟨v ∧ t, u⟩ ⟨u ∧ s, v⟩ ⟨x 11 , x 1 ⟩ . (C.2) But ⟨x 22 , x 1 ⟩ = ∂ 2 ⟨x 2 , x 1 ⟩ -⟨x 2 , x 12 ⟩ = ∂ 2 ⟨x 2 , x 1 ⟩ - 1 2 ∂ 1 ⟨x 2 , x 2 ⟩ = - 1 2 ∂ 1 g 22 , (C.3)
since dg 12 = 0. Also, ⟨x 11 , x 1 ⟩ = ∂ 1 g 11 /2. Projecting over x 2 yields similar relations so that overall

- dg 22 2 = ⟨u ∧ s, t⟩ ⟨v ∧ t, s⟩ ⟨v ∧ t, u⟩ ⟨u ∧ s, v⟩ dg 11 2 . (C.4)
Hence, the in-plane Poisson's coefficient is

ν = - dg 22 /g 22 dg 11 /g 11 = g 11 g 22 ⟨u ∧ s, t⟩ ⟨v ∧ t, s⟩ ⟨v ∧ t, u⟩ ⟨u ∧ s, v⟩ = κ 2 κ 1 (C.5)
and matches, up to a sign, the out-of-plane Poisson's coefficient -κ 2 /κ 1 .

Appendix D. The Poisson's ratios of the "Mars" pattern

The "Mars" pattern is a developable flat-foldable zigzag sum. Let us suppose that the pattern is equilateral so that u = v = s = t. Now we have ⟨u, s⟩ = -⟨v, t⟩ = cos α, ⟨u, t⟩ = -⟨v, s⟩ = cos β.

(D.1)

We are at liberty to adopt a mountain-valley assignment convention: we let creases s, t and v be mountains and let crease u be a valley. Then, inspection of (A.3) reveals that s w = -t w . This allows to greatly simplify (A.5) ultimately into

g 22 = 4 -4 (cos α + cos β) 2 g 11 . (D.2)
The Poisson's coefficient is deduced by differentiation. It reads

ν = -4 (cos α + cos β) 2 g 11 g 22 (D.3)
and is automatically equal to κ 2 /κ 1 .

Recall that g 11 and g 22 are bounded. In the maximally unfolded state, cos(u, v) = -cos(α -β) so that, in any case, g 11 remains smaller than 2 + 2 cos(α -β). The other bounds are similarly deduced and, overall, we have

4 cos 2 α + β 2 ≤ g 11 ≤ 4 cos 2 α -β 2 , 4 sin 2 α -β 2 ≤ g 22 ≤ 4 sin 2 α + β 2 .
(D.4)

Appendix E. "Mars"-folded hyperboloids

Consider an axisymmetric origami pillar folded out of an equilateral "Mars" pattern with a parametrization x ∈ R 3 of the form x = (ρ(ξ 2 ) cos(qξ 1 ), ρ(ξ 2 ) sin(qξ 1 ), z(ξ 2 )).

(E.1)

We denote L 1 and L 2 the dimensions of the origami pattern in its flat, maximally unfolded, reference state. We then let ξ 1 ∈ [-S 1 /2, S 1 /2] and ξ 2 ∈ [-S 2 /2, S 2 /2] with

S 1 = L 1 2 cos((α -β)/2) , S 2 = L 2 2 sin((α + β)/2) . (E.2)
The pillar is obtained by folding the pattern then glueing together the ends ξ 1 = 0 and ξ 1 = S 1 . In other words, we let q = 2π/S 1 .

With these notations, it straightforward to see that

x 1 = (-qρ(ξ 2 ) sin(qξ 1 ), qρ(ξ 2 ) cos(qξ 1 ), 0),

x 2 = (ρ ′ (ξ 2 ) cos(qξ 1 ), ρ ′ (ξ 2 ) sin(qξ 1 ), z ′ (ξ 2 )),
x 11 = (-q 2 ρ(ξ 2 ) cos(qξ 1 ), -q 2 ρ(ξ 2 ) sin(qξ 1 ), 0),

x 22 = (ρ ′′ (ξ 2 ) cos(qξ 1 ), ρ ′′ (ξ 2 ) sin(qξ 1 ), z ′′ (ξ 2 )), (E.3)
with ρ ′ ≡ dρ/dξ 2 , ρ ′′ ≡ dρ ′ /dξ 2 and so on. Then,

g 11 = q 2 ρ 2 and g 22 = ρ ′2 + z ′2 (E.4)
depend on one another through the in-plane kinematical constraint

ρ ′2 + z ′2 = 4 -4 (cos α + cos β) 2 q 2 ρ 2 . (E.5)
The out-of-plane constraint, i.e., x 11 ∥ x 22 , immediately yields z ′′ = 0. (E.6)

Together, these two ODEs can be integrated into

z ′ = cst, ρ = 4c 2 (ξ 2 -ξ o ) 2 + (cos α + cos β) 2 q 2 c 2 , (E.7)
with c 2 = 1 -z ′2 /4. Hence, the origami pillar embraces a hyperboloid. The pattern does not cover the whole hyperboloid however and is limited to the band spanned by ξ 2 ∈ [-S 2 /2, S 2 /2]. Overall then, and up to rigid body motions, the folded pillar has two degrees of freedom. The first is c: it describes the shape of the hyperboloid; as it goes from 1 to 0, the hyperboloid changes from a "yoyo" into a cylinder. The second is ξ o : it specifies how the band of length L 2 that the pattern covers is centered on or offset away from the equator. Not all values (c, ξ o ) are admissible however since candidate solutions must further satisfy the upper and lower bounds weighing on the metric, namely (D.4). In the present case, these read

4 cos 2 ((α + β)/2) ≤ q 2 ρ 2 ≤ 4 cos 2 ((α -β)/2). (E.8)
For a centered band (ξ o = 0) in particular, these bounds reduce to

4 cos 2 α + β 2 ≤ (cos α + cos β) 2 c 2 , c 2 q 2 S 2 2 + (cos α + cos β) 2 c 2 ≤ 4 cos 2 α -β 2 , (E.9)
and place a maximum bound on the aspect ratio S 2 /S 1 = qS 2 /(2π) of origami pillars that can be formed. In the main text, the angle ∠(u, v) at the equator is used rather than c. These two are in a one-to-one correspondence

2 -2 cos(u, v) = g 11 = (cos α + cos β) 2 c 2 . (E.10)

Appendix F. Bending energy -General case

We have found that the admissible x µν satisfy five linear constraints: three describe the proportionality of x 11 and x 22 , and two reduce to dg 12 = 0. Thus, the admissible x µν belong to a four-dimensional linear subspace of (R 3 ) 3 space. Each dimension is spanned by a DOF that we attribute to one planarity defect δ (ij) , (i, j) ∈ {1, 1} 2 such that δx (11) = δ (11) u ∧ s, δx ( 11) = δ ( 11) s ∧ v, δx ( 11 ) = δ ( 11 ) v ∧ t,

δx (1 1) = δ (1 1) t ∧ u. (F.1)
Here too, we began by defaulting the same displacements as in (B.1) to zero. Next, we derive expressions for the δ (ij) in function of the x µν .

Generalizing equation (B.6), it is possible to derive full expressions for the x µν in terms of the δ (ij) . We find,

x 11 = ⟨u, s ∧ v⟩ ⟨t, s ∧ v⟩ δ (11) v ∧ t + ⟨v, u ∧ s⟩ ⟨t, u ∧ s⟩ δ ( 11) t ∧ u + ⟨v, t ∧ u⟩ ⟨s, t ∧ u⟩ δ ( 11 ) u ∧ s + ⟨u, v ∧ t⟩ ⟨s, v ∧ t⟩ δ (1 1) s ∧ v. (F.2)
As for x 12 , it is readily available and reads

x 12 = δ (11) u ∧ s -δ ( 11) s ∧ v + δ ( 11 ) v ∧ t -δ (1 1) t ∧ u. (F.3)

These two vector equations provide a system (with two redundant equations), which can be solved for the δ (ij) . For instance, projecting x 11 over u and x 12 over v leads to a 2 × 2 system ⟨x 11 , u⟩ = ⟨u, s ∧ v⟩ ⟨u, v ∧ t⟩ ⟨t, s ∧ v⟩ (δ (11) + δ (1 1) ),

⟨x 12 , v⟩ = ⟨u, s ∧ v⟩ δ (11) -⟨u, v ∧ t⟩ δ (1 1) , (F.4)
which can be solved for δ (11) and δ (1 1) to give

δ (11) = ⟨x 11 , u⟩ ⟨t, s ∧ v⟩ + ⟨x 12 , v⟩ ⟨u, s ∧ v⟩ ⟨u, s ∧ v⟩ ⟨t -s, u ∧ v⟩ , δ (1 1) = ⟨x 11 , u⟩ ⟨t, s ∧ v⟩ -⟨x 12 , v⟩ ⟨u, v ∧ t⟩ ⟨u, v ∧ t⟩ ⟨t -s, u ∧ v⟩ . (F.5)
Similar considerations lead to ⟨x 11 , v⟩ = ⟨v, u ∧ s⟩ ⟨v, t ∧ u⟩ ⟨t, u ∧ s⟩ (δ ( 11) + δ ( 11 ) ), ⟨x 12 , u⟩ = -⟨u, s ∧ v⟩ δ ( 11) + ⟨u, v ∧ t⟩ δ ( 11 ) , (F.6) and to δ ( 11) = ⟨x 11 , v⟩ ⟨t, u ∧ s⟩ -⟨x 12 , u⟩ ⟨v, u ∧ s⟩ ⟨v, u ∧ s⟩ ⟨t -s, u ∧ v⟩ , δ ( 11 ) = ⟨x 11 , v⟩ ⟨t, u ∧ s⟩ + ⟨x 12 , u⟩ ⟨v, t ∧ u⟩ ⟨v, t ∧ u⟩ ⟨t -s, u ∧ v⟩ .

(F.7)

Most importantly, the planarity defects δ (ij) are linear forms of the full, in-and out-of-plane, components of the parametrization's second derivatives. Now each panel in the unit cell, contributes a term g (ij) (r 2 δ (ij) ), (i, j) ∈ {1, 1} 2 , to the bending energy density of the tessellation equal to the bending energy of the relevant panel for a given planarity defect r 2 δ (ij) . The potential g (ij) can be nonlinear in principle, but given that the planarity defects are of order O(r 2 ), it is reasonable to linearize it in the vicinity of r → 0, that is while assuming, at the same time, that the plane is the natural state of the panels. In conclusion, the bending energy density of the tessellation takes the form b

= (i,j)∈{1, 1} 2 1 2 D (ij) (δ (ij) ) 2 (F.8)
where the D (ij) are the panels flexural rigidities normalized with respect to the area of a reference unit cell and where the δ (ij) have been shown to be configurationdependent (i.e., (u, v, s, t)-dependent) linear forms of the second derivatives x µν . Equivalently, bending energy density b = b(Γ iµν ; g αβ ) can be written as a metric-dependent quadratic form of the Christoffel symbols and of the curvatures. This is done next in the particular case of the "Mars" pattern.

Note that while we linearized b into a quadratic form in the limit r → 0 (i.e., for infinitesimal panel bending), we make no a priori assumptions regarding the macroscopic curvatures (i.e., the x µν ). These are free to take finite values. Last, the flexural rigidities D (ij) are material constants. Here, we suppose for simplicity that they are all equal to D.

Appendix G. Bending energy -The "Mars" pattern

Similar considerations to the ones that led to equation (A.3) show that the current configuration of the folds (u, v, s, t) can be determined from the tangent vectors (x 1 , x 2 ). Indeed, we have, in the case of an equilateral "Mars" pattern, (G.5)

u = 1 2 x 1 + cos α -cos β g 22 x 2 -u n n, v = - 1 2 x 1 + cos α -
In conclusion, the total strain energy of a "Mars" pattern is

ψ(x) = Ω b(Γ; g)dξ 1 dξ 2 (G.6)
where the bending strain energy density is b(Γ; g) = D 2 (δ (11) ) 2 + (δ ( 11) ) 2 + (δ (1 1) ) 2 + (δ ( 11 ) ) 2

(G.7) with the planarity defects δ (ij) being the metric-dependent linear forms given in (F.5) and (F.7). Therein, note how the rigidity coefficients in front of the in-plane strain-gradient are proportional to the term 1/u n which diverges when g 11 and g 22 approach their maximum and minimum values, i.e., at the boundary of the domain of kinematical admissibility. Meanwhile, the rigidity coefficients weighing the curvatures remain bounded.

Figure 2 :

 2 Figure 2: Zigzag sums. (a) A zigzag sum in a periodic reference configuration: translating one zigzag along another sweeps a polyhedral surface corresponding to a crease pattern complete with a natural assignment of mountain and valley folds. (b) A super-cell in an abstract reference state: indices (ij) are discrete coordinates for the vertices; the central node is at ξ ≡ ξ (00) . (c) A super-cell in a physical folded state: vectors (u, v, s, t) are ξ-dependent non-dimensional vectors aligned with the creases; they fully characterize fold-only configurations as well as the tangent vectors x 1 and x 2 . (d) A super-cell following an infinitesimal bending δx of the panels: note how the parallelograms are slightly distorted now. (e) A global view of a finitely deformed zigzag sum (Σ): each super cell is designated by its central node's coordinates ξ = (ξ 1 , ξ 2 ) (red dot) and is obtained by the process illustrated in panels (b-d). (f) An orthogonal zigzag sum: initially orthogonal zigzag directions remain orthogonal during folding (g 12 ≡ 0). (g) An oblique zigzag sum: initially non-orthogonal zigzag directions are sheared during folding; even then, g 12 remains constant.

Figure 3 :

 3 Figure 3: Convergence analysis. (a) Convergence speed: the maximum error in equation (7) decreases like r or, equivalently, like the reciprocal of the number of unit cells per row. The plot also shows that the maximum angular deflection rδ(ij) decreases in the same fashion confirming the fold-to-bend orders of magnitude adopted in (1) and that, in the limit r → 0, all panels are flat. The inset shows the unit cell adopted for this example. (b,c) Two iterates: snapshots of the way in which the zigzag sum converges to a smooth surface taken for 64 and 128 unit cells per row. Note that the zigzag sum is finitely and non-uniformly stretched and bent. Illustrated surfaces are constructed by "propagating" inextensibility constraints away from a given boundary; see[START_REF] Nassar | Fitting surfaces with the Miura tessellation[END_REF] for details.

  3 × R 3 spanned by 4 DOFs, each attributed to a planarity defect (δ (ij) ) of one of the four panels of the central unit cell (Fig 3a) [23, App. B]. By the rank-nullity theorem, the x µν satisfy 9 -4 = 5 linear compatibility equations. Straightforward calculations show that the first two of these equations are g 12,1 = g 12,2 = 0, i.e., they reproduce dg 12 = 0. The other three relations are far more interesting; they combine into the vector identity [23, App. B]

  can be found in the Supplemental Material [23, App. F, G].

Figure 4 :

 4 Figure 4: "Mars"-folded origami pillars are one-sheeted hyperboloids. (a) A unit cell of the "Mars" pattern: the supplementary opposite angles ensure developability and flat-foldability; an equilateral pattern is further rectangular. (b) Action of the first DOF on the shape of the pillar: given a small enough aspect ratio L 2 /L 1 , the pillar can be inverted inside-out; here, the pattern has N = 100 unit cells per meridian, and 8 cells per parallel, with a flat-unfolded aspect ratio L 2 /L 1 = 8%. (c) Action of the second DOF: the pattern flattens as the opening angle ∠(u, v) at a given parallel is increased; here, the pattern has N = 100 unit cells per meridian and 40 unit cells per parallel with L 2 /L 1 = 40%.(d) Bending energy: plot of normalized total bending energy as a function of the opening angle ∠(u, v) at the equator for origami pillars centered about the equator with L 2 /L 1 = 40%; the discretely summed energy for N = 100 is plotted against the continuum limit N → ∞; data points (i -iv) are the ones illustrated on panel (c); greyed out domain cannot be accessed with the current aspect ratio; the continuum theory satisfactorily predicts both energy and kinematical admissibility. Adopted values: α = π/3, β = π/2.

  cos β g 22x 2 -u n n, function of the metric tensor. Thus, the coefficients appearing in the expression of b can be written in terms of the metric tensor as well⟨s, u ∧ v⟩ = -⟨t, u ∧ v⟩ = u n √ g/2, ⟨u, s ∧ t⟩ = ⟨v, s ∧ t⟩ = -cos α + cos β g 11 √ gu n , (G.3)with g = det g = g 11 g 22 . As for the remaining terms involving the x µν , they can be expanded into combinations of the Christoffel symbols Γ σµν ≡ ⟨x µν , x σ ⟩ and of the coefficients of the second fundamental form Γ 3µν ≡ ⟨x µν , n⟩.Recall that, in the present case where the curvilinear coordinates are rectangular, the Christoffel symbols read

	s = t = u n = 1 -cos α + cos β g 11 cos α + cos β g 11 g 11 4 -x 1 + x 1 -(cos α -cos β) 2 1 2 x 2 , 1 2 x 2 , g 22 being a known Namely, we have with ⟨x 11 , u⟩ = 1 2 Γ 111 + cos α -cos β g 22 Γ 211 -u n Γ 311 , ⟨x 11 , v⟩ = -1 2 Γ 111 + cos α -cos β g 22 Γ 211 -u n Γ 311 , ⟨x 12 , u⟩ = 1 2 Γ 112 + cos α -cos β g 22 Γ Γ 111 = 1 2 g 11,1 , Γ 211 = -1 2 g 11,2 , Γ 112 = 1 2 g 11,2 , Γ 212 = 1 2 g 22,1 , Γ 122 = -1 2 g 22,1 , Γ 222 = 1 2 g 22,2 .	(G.1) (G.2)

212 -u n Γ 312 , ⟨x 12 , v⟩ = -1 2 Γ 112 + cos α -cos β g 22 Γ 212 -u n Γ 312 .

(G.4)

⟨u -v, s -t⟩ ⟨u -v, s -t⟩ ∥s -t∥ 2 .(3)Note that bending-related contributions disappear in the limit r → 0 and can be disregarded for now.

Acknowledgements

HN acknowledges support by the NSF under CAREER award No. CMMI-2045881. AL acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant No. ANR-17-CE08-0039 (project ArchiMatHOS). EW acknowledges support by the NASA Missouri Space Grant Consortium.

Appendix A. The in-plane Poisson's ratio: General case

The metric g of a zigzag sum has the components g 11 = ∥u -v∥ 2 , g 22 = ∥s -t∥ 2 and g 12 = ⟨u -v, s -t⟩. Note that the scalar products ⟨u, s⟩, ⟨u, t⟩, ⟨v, s⟩ and ⟨v, t⟩ are combinations of lengths and angles that are isometrically preserved; they are motion constants and so is g 12 . It is somewhat more challenging to derive a relationship between g 11 and g 22 ; this is done next. First, note that s and t can be decomposed into

with

Components s w and t w are easily determined by considering the magnitudes of s and t. Indeed, we have

with u = ∥u∥, and so on. Therein, note that the only variables are

(A.4) and are both functions of g 11 . Last, we write

and accordingly deduce g 22 as a function of g 11 , as well as ν as a function of g 11 . It is worth mentioning that g 11 and g 22 have upper and lower bounds corresponding to some creases being maximally folded or unfolded. Beyond these bounds the pattern would penetrate itself or break apart. These situations can be avoided by enforcing s 2 w > 0 and t 2 w > 0.

Appendix B. The admissible second-order derivatives

Consider the vertices of a super cell placed at x (ij) = X (ij) + r 2 δx (ij) where the X (ij) describe a pre-folded state and the r 2 δx (ij) are bending-induced perturbations. Similarly, the central vertex ξ is at X subsequent to the fold-only


A • Lebée

G • Corre

M K • Ferradi

K • Sab

X Cespedes

Civil engineering applications of the Asymptotic Expansion Load Decomposition beam model, an overview

Keywords: Model reduction, Higher-order beam model, Asymptotic expansion, Plasticity, Structural analysis

The Asymptotic Expansion Load Decomposition higher-order beam model is based on the classical two-scale asymptotic expansion in the linear elasticity framework. It was successively extended to eigenstrains and to plasticity in small deformations in different papers. The present paper offers a comprehensive and consistent presentation of our approach applied to civil engineering applications.

Introduction

Until very recently, most structural engineering offices were designing structures based either on very simple beam models or on complex finite elements modeling. Thanks to the recent development of cloud-computing, a new business model is emerging: online structural analysis. Indeed, cloud computing gives structural engineers access to more powerful computing platforms. It also reduces maintenance operations for the provider since only one system is maintained instead of adapting the software to multiple platforms (Windows, Unix, multiple versions of them, etc.). This is the path chosen by the start-up Strains engineering co-founded by Xavier Cespedes in 2014.

The development of such tools requires interactive and real-time simulation softwares.

Hence, an efficient mechanical modeling of structural elements and connections is critical.

One of the main project of this start-up is the development of a higher-order beam model capable of capturing in more details the influence of the applied load or the boundary conditions for thin-walled and heterogeneous beams. Indeed, traditional design based on simple beam modeling presents the advantage to be easy to understand and to allow rapid verification of the results. However, it is not able to give detailed information about local distribution of stresses. In order to obtain more detailed results, complex finite element modeling is necessary. Though, this requires expensive and black-boxed software and is computationally more demanding.

The present paper summarizes and presents a critical review of the Asymptotic Expansion

Load Decomposition beam model developped and implemented for Strains engineering and published in several contributions [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF][START_REF] Corre | Higher-order beam model with eigenstrains: theory and illustrations[END_REF][START_REF] Corre | The Asymptotic Expansion Load Decomposition elastoplastic beam model[END_REF]. This model originates from the seminal work of Miara and Trabucho [START_REF] Miara | A Galerkin spectral approximation in linearized beam theory[END_REF] for linear elastic beams and was first directly applied by Ferradi et al [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF]. Because structural analysis must also take into account eigenstrains such as thermal load, prestress or weakly coupled swelling phenomena, this new type of loading was treated by Corre et al [START_REF] Corre | Higher-order beam model with eigenstrains: theory and illustrations[END_REF] in a linear elastic framework. Finally, since the total plastic strain may be considered as an eigenstrain, the Asymptotic Expansion Load Decomposition model was extended to elasto-plasticity assuming small perturbations in [START_REF] Corre | The Asymptotic Expansion Load Decomposition elastoplastic beam model[END_REF].

The paper is organized in two sections. First, the higher-order beam model is derived in the framework of linear elasticity for a given load distribution and some illustrations are briefly presented. Second, the higher-order beam model is extended to elasto-plasticity assuming small perturbations and is validated against a full 3D computation. Finally, limitations and outlooks of the approach are presented.

The linear elastic Asymptotic Expansion Load Decomposition beam model

Most beam models are based on ad-hoc assumptions on the 3D fields which motivated the denomination axiomatic. They rely on an educated guess on the 3D displacement field in a separated form between the longitudinal coordinate and the transverse coordinates. Then, straightforward application of the minimum of potential energy leads to 1D boundary value problems corresponding to the beam model. Numerous strategies were devised for building the kinematics of the beam. Remarkably, the extended Saint Venant solution from Iesan [START_REF] Iesan | Saint-Venant's problem for inhomogeneous and anisotropic elastic bodies[END_REF] and the formal asymptotic expansion method [START_REF] Trabucho | Mathematical modelling of rods[END_REF] delivers almost the same collection of cross-section displacements. Based on the seminal idea from Vogelius and Babuška [START_REF] Vogelius | On a Dimensional Reduction Method I. The Optimal Selection of Basis Functions[END_REF][START_REF] Vogelius | On a Dimensional Reduction Method II. Some Approximation-Theoretic Results[END_REF] which also originated the family of "hierarchical models" for plates and shells [START_REF] Actis | Hierarchic models for laminated plates and shells[END_REF], Miara and Trabucho [START_REF] Miara | A Galerkin spectral approximation in linearized beam theory[END_REF] made two noticeable observations regarding this beam kinematics. First, the formal asymptotic expansion delivers a free family of kinematic enrichment which is dense in the space of the 3D solution. This means that going sufficiently high in the expansion allows arbitrary refinement of the 3D solution. Second, the truncation of this family ensures that the corresponding beam model is asymptotically consistent except at the boundary. Indeed, for loads which do not generate boundary layers, Miara and Trabucho [START_REF] Miara | A Galerkin spectral approximation in linearized beam theory[END_REF] proved a higher-order convergence result. This means that the kinematic enrichment delivered by the formal asymptotic expansion is optimal in terms of approximation error far from the extremities of the beam.

In the following, the asympotic expansion procedure for deriving the AELD kinematics is recalled for a beam under body loads and its link with classical beam models is highlighted.

Then, the higher-order beam model is derived and two applications are presented.

The kinematic enrichment based on the formal asymptotic expansion

In this section, the formal asymptotic expansion procedure is recalled for a linear elastic beam loaded by an arbitrary body force. It is an opportunity to emphasize its close link with Saint-Venant solution as well as the extension from Iesan [START_REF] Iesan | Saint-Venant's problem for inhomogeneous and anisotropic elastic bodies[END_REF]. It appears that a series of 3D displacement fields distributed over the cross-section may be obtained as function of the classical beam generalized displacements as well as the longitudinal distribution of the applied load and its higher gradients. Such a basis of "sections modes" may be derived for any given applied load [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF] and was extended to prescribed eigenstrains by Corre et al [START_REF] Corre | Higher-order beam model with eigenstrains: theory and illustrations[END_REF].

The 3D problem

𝑥 2 𝑥 1 𝑥 3 𝐿 ℎ 𝜕Ω 𝑡 Ω S + S -

Fig. 1: The beam 3D configuration

We consider a beam occupying the prismatic domain Ω (Figure 1) with a length 𝐿 and a cross-sectional typical size ℎ. The boundary 𝜕Ω is the union of the lateral (free) surface 𝜕Ω 𝑡 and the two end sections S ± (clamped). The longitudinal coordinate is 𝑥 3 and the transverse coordinates are 𝑥 1 and 𝑥 2 denoted as 𝑥 𝛼 , the corresponding reference frame is

denoted (𝑂, 𝒆 -1 , 𝒆 -2 , 𝒆 -3
) where 𝑂 is an arbitrary point of the plane 𝑥 3 = 0.

The constitutive material of the beam, relating the stress tensor 𝝈 to the linearized deformation 𝜺, is only function of the transverse coordinates 𝑥 𝛼 and invariant in the longitudinal direction. For convenience and without limitation, the corresponding fourth order stiffness tensor 𝑪 (𝑥 𝛼 ) is assumed monoclinic with respect to a plane of normal 𝒆 -3 :

𝐶 𝛼𝛽𝛾3 = 𝐶 333𝛼 = 0. ( 1
)
The beam is only subjected to a body load 𝑓 𝑖 . The corresponding 3D linear boundary value problem writes as:

           div -(𝝈) + 𝒇 -= 0 on Ω, 𝝈 = 𝑪 : 𝜺 𝒖 - on Ω, 𝝈 • 𝒏 -= 0 on 𝜕Ω 𝑡 , 𝒖 -= 0 on S ± , (2)
where 𝒏 -is the outer normal to 𝜕Ω 𝑡 , div -is the classical divergence operator and 𝜺 𝒖 -is the symmetric part of the displacement gradient.

Scaling and variable separation A new set of coordinates 𝑦 𝑖 is defined from the global coordinates,

(𝑥 1 , 𝑥 2 , 𝑥 3 ) = (ℎ𝑦 1 , ℎ𝑦 2 , 𝐿𝑦 3 ) , (3)
which rewrites the derivation operator as:

∇ 𝒙 -= 1
𝐿 ∇ 𝑦 3 + 1 𝜂 ∇ 𝑦 𝛼 , (4)
as well as the integration over the domain Ω as:

∫ Ω 𝑑Ω = ∫ 𝐿 0 ∫ S 𝑑𝑥 𝛼 𝑑𝑥 3 = 𝐿 3 𝜂 2 ∫ 1 0 ∫ S 0 𝑑𝑦 𝛼 𝑑𝑦 3 = 𝐿 3 𝜂 2 ∫ 1 0 𝑑𝑦 3 , (5)
where 𝜂 = ℎ 𝐿 is the small parameter related to the slenderness of the beam, S 0 is the scaled cross-section, independent from 𝜂 and denotes integration on the scaled cross-section.

The body load is assumed to have the following separated form:

𝑓 𝛼 = 𝜂 2 𝐿 𝑓 𝛼 (𝑦 𝜂 )𝐹 (𝑦 3 ) and 𝑓 3 = 𝜂 𝐿 𝑓 3 (𝑦 𝜂 )𝐹 (𝑦 3 ) (6)
where 𝑓 𝑖 are cross-section distributions of the load which are modulated by a single longitudinal function 𝐹 (𝑥 3 ). In this section, capital letters denote functions of only the 𝑦 3 coordinate (except for 𝑪) and • denotes functions of only cross-section coordinates 𝑦 𝛼 .

Expansion The asymptotic expansion method is a formal procedure in which all fields are assumed sufficiently smooth. It yields a cascade of cross-section and longitudinal boundary value problems which are classically solved recursively. In the present case, only the crosssection problems are of interest in order to derive a collection of displacement modes.

The displacement, strain and stress variables are expanded as power series of the small parameter as follows [START_REF] Sanchez-Palencia | Non-Homogeneous Media and Vibration Theory[END_REF][START_REF] Trabucho | Mathematical modelling of rods[END_REF][START_REF] Buannic | Higher-order effective modeling of periodic heterogeneous beams. I. Asymptotic expansion method[END_REF]:

𝒖 -= 𝐿 𝑈 0 𝛼 (𝑦 3 )𝒆 -𝛼 + 𝜂𝒖 - 1 + 𝜂 2 𝒖 - 2 + ... , (7)
𝜺 = 𝜺 0 + 𝜂𝜺 1 + 𝜂 2 𝜺 2 + ..., (8)
𝝈 = 𝝈 0 + 𝜂𝝈 1 + 𝜂 2 𝝈 2 + ... (9)
and introduced in the rescaled equations of the 3D boundary value problem [START_REF] Amrouche | On the characterizations of matrix fields as linearized strain tensor fields[END_REF] where each power 𝑝 of 𝜂 is identified. In Equation 7 The problem being linear, the choice of the starting order has no incidence on the final formulation in terms of physical variable. Here the starting order is chosen so that the leading order of the displacement field is 0. The starting order of the other fields is chosen accordingly and motivates the scaling of the load in Equation [START_REF] Buannic | Higher-order effective modeling of periodic heterogeneous beams. I. Asymptotic expansion method[END_REF].

For 𝑝 ∈ N, each compatibility equations, boundary conditions and constitutive equations for 𝑝 and equilibrium equations for 𝑝 -1 yield an auxiliary problem on the cross-section which splits in two uncoupled boundary value problems.

Transverse displacement First, the cross-section displacement problems (transverse mode) T 𝑝+1 are gathered for 𝑝 ≥ 0:

T 𝑝+1 :                        𝜎 𝑝 𝛼𝛽,𝛽 + 𝜎 𝑝-1 𝛼3,3 + 𝛿 𝑝3 𝑓 𝛼 𝐹 = 0 on S 0 , ( 10a
)
𝜎 𝑝 𝛼𝛽 = 𝐶 𝛼𝛽𝛾 𝛿 𝜀 𝑝 𝛿𝛾 + 𝐶 𝛼𝛽33 𝜀 𝑝 33 on S 0 , ( 10b
)
𝜎 𝑝 33 = 𝐶 33𝛼𝛽 𝜀 𝑝 𝛽 𝛼 + 𝐶 3333 𝜀 𝑝 33 on S 0 , ( 10c
)
𝜀 𝑝 𝛼𝛽 = 𝑢 𝑝+1 ( 𝛼,𝛽) , 𝜀 𝑝 33 = 𝑢 𝑝 3,3 on S 0 , ( 10d
)
𝜎 𝑝 𝛼𝛽 𝑛 𝛽 = 0 on 𝜕S 0 . ( 10e
)
where 𝝈 -1 = 0 and 𝛿 3 𝑝 = 1 if 𝑝 = 3 and 𝛿 3 𝑝 = 0 else. Transposing the results from [START_REF] Ciarlet | Another approach to linearized elasticity and Korn's inequality[END_REF][START_REF] Amrouche | On the characterizations of matrix fields as linearized strain tensor fields[END_REF],

for a simply connected cross-section and regular enough 𝑪 and loadings, this boundary value problem on the displacement 𝑢 𝑝+1 𝛼 is a pure traction (Neumann) problem which is well-posed provided that the applied load is globally self-equilibrating for cross-section translations and rotation:

𝜎 𝑝-1 𝛼3,3 + 𝛿 𝑝3 𝑓 𝛼 𝐹 = 0 and 𝑦 𝛽 𝜖 𝛽 𝛼 𝜎 𝑝-1 𝛼3,3 + 𝛿 𝑝3 𝑓 𝛼 𝐹 = 0, (11)
where 𝜖 𝛼𝛽 is the permutation operator: 𝜖 11 = 𝜖 22 = 0, 𝜖 12 = +1, 𝜖 21 = -1. Then, under compatibility conditions [START_REF] Corre | A new higher-order elastoplastic beam model for reinforced concrete[END_REF], the solution is uniquely defined up to the following rigid motion of the section in its plane:

𝑢 𝑅, 𝑝+1 𝛼 = 𝑈 𝑝+1 𝛼 (𝑦 3 ) + 𝑦 𝛽 𝜖 𝛽 𝛼 Θ 𝑝+1 (𝑦 3 ). ( 12
)
where 𝑈 𝑝+1 𝛼 is a transverse displacement and Θ 𝑝+1 a twist rotation.

Longitudinal displacement Second, the longitudinal displacement problems (warping mode) W 𝑝 are obtained for 𝑝 ≥ 0:

W 𝑝+1 :                  𝜎 𝑝 3𝛼, 𝛼 + 𝜎 𝑝-1 33,3 + 𝛿 𝑝2 𝑓 3 𝐹 = 0 on S 0 , ( 13a
)
𝜎 𝑝 𝛼3 = 𝐶 𝛼3𝛽3 2𝜀 𝑝 𝛽3 on S 0 , ( 13b
) 2𝜀 𝑝 𝛼3 = 𝑢 𝑝+1 3, 𝛼 + 𝑢 𝑝 𝛼,3 on S 0 , ( 13c
)
𝜎 𝑝 𝛼3 𝑛 𝛼 = 0 on 𝜕S 0 . ( 13d
)
Again, for a simply connected cross-section and regular enough 𝑪 and loadings, this boundary value problem on the displacement 𝑢 𝑝+1 3 is well-posed if the applied load is globally self-equilibrating for the longitudinal translation:

𝜎 𝑝-1 33,3 + 𝛿 𝑝2 𝑓 3 𝐹 = 0. (14)
In this case, the solution is uniquely defined up to a uniform longitudinal displacement:

𝑢 𝑅, 𝑝+1 3 = 𝑈 𝑝+1 3 (𝑦 3 ). ( 15
)
Resultants and macroscopic equilibrium equations The rigid motion of the section suggests the following definition of the beam resultants at each order 𝑝 ≥ 0:

𝑁 𝑝 3 = 𝜎 𝑝 33 , 𝑀 𝑝 𝛼 = 𝑦 𝛼 𝜎 𝑝 33 , 𝑀 𝑝 3 = 𝑦 𝛽 𝜖 𝛽 𝛼 𝜎 𝑝 𝛼3 and 𝑉 𝑝 𝛼 = 𝜎 𝑝 𝛼3 , (16)
where 𝑁 𝑝 3 is the normal traction, 𝑀 𝑝 𝛼 are the bending moments , 𝑀 𝑝 3 is the moment of torsion and 𝑉 𝑝 𝛼 are the shear forces.

From the compatibility conditions (10a) and (13a) it is possible to prove that these resultants must comply with the following classical beam equilibrium equations for each 𝑝 ≥ 0:

             𝑁 𝑝 3,3 + 𝛿 𝑝1 𝑝 3 𝐹 = 0, (17a)
𝑀 𝑝 𝛼,3 -𝑉 𝑝+1 𝛼 + 𝛿 𝑝1 𝜇 𝛼 𝐹 = 0, ( 17b
)
𝑀 𝑝 3,3 + 𝛿 𝑝2 𝜇 3 𝐹 = 0, ( 17c
)
𝑉 𝑝 𝛼,3 + 𝛿 𝑝2 𝑝 𝛼 𝐹 = 0. ( 17d
)
It will appear that 𝑀 𝑝 𝛼 is the working conjugate to the curvature 𝑈 𝑝 𝛼,33 and not the conventional bending moment. Indeed, the classical definition is 𝑚 𝑝 𝛼 = 𝜖 𝛼𝛽 𝑦 𝛽 𝜎 𝑝 33 = 𝜖 𝛼𝛽 𝑀 𝑝 𝛽 . This choice is made for convenience.

where 𝑝 𝛼 = 𝑓 𝛼 are the transverse loads, 𝑝 3 = 𝑓 3 is the longitudinal load, 𝜇 𝛼 = 𝑦 𝛼 𝑓 3 are distributed bending couples and 𝜇 3 = 𝑦 𝛽 𝜖 𝛽 𝛼 𝑓 𝛼 is a distributed torsion load.

The series of problems are now solved order by order.

First-order problems

Transverse displacement The problem T 1 is not loaded. Consequently, the transverse displacement 𝑢 1 𝛼 is a rigid motion and the corresponding stress is null:

𝑢 1 𝛼 = 𝑈 1 𝛼 (𝑦 3 ) + 𝑦 𝛽 𝜖 𝛽 𝛼 Θ 1 (𝑦 3 ) and 𝜎 0 𝛼𝛽 = 0, 𝜎 0 33 = 0. ( 18
)
Here, Θ 1 appears as the leading order angle of twist and 𝑈 1 𝛼 as the second order macroscopic transverse displacement.

Longitudinal displacement

The longitudinal displacement problem (warping mode) W 1 writes as:

W 1 :                𝜎 0 3𝛼, 𝛼 = 0 on S 0 , ( 19a
)
𝜎 0 𝛼3 = 𝐶 𝛼3𝛽3 2𝜀 0 𝛽3 on S 0 , ( 19b
) 2𝜀 0 𝛼3 = 𝑢 1 3, 𝛼 + 𝑈 0 𝛼,3 on S 0 , ( 19c
)
𝜎 0 𝛼3 𝑛 𝛼 = 0 on 𝜕S 0 . ( 19d
)
The applied load is self-equilibrating and the solution of this boundary value problem writes as:

𝑢 1 3 = 𝑈 1 3 + 𝑦 𝛼 𝑈 0 𝛼,3 and 𝜎 0 𝛼3 = 0, (20)
where 𝑈 0 𝛼,3 appears as the bending rotation and 𝑈 1 3 as the leading order longitudinal displacement.

Second-order problems

Transverse displacement The transverse displacement 𝑢 2 𝛼 is derived through:

T 2 :                      𝜎 1 𝛼𝛽,𝛽 = 0 on S 0 , ( 21a
)
𝜎 1 𝛼𝛽 = 𝐶 𝛼𝛽𝛾 𝛿 𝜀 1 𝛿𝛾 + 𝐶 𝛼𝛽33 𝜀 1 33 , on S 0 , ( 21b
)
𝜎 1 33 = 𝐶 33𝛼𝛽 𝜀 1 𝛽 𝛼 + 𝐶 3333 𝜀 1 33 on S 0 , (21c)
𝜀 1 𝛼𝛽 = 𝑢 2 ( 𝛼,𝛽) , 𝜀 1 33 = 𝑈 1 3,3 + 𝑦 𝛼 𝑈 0 𝛼,33 on S 0 , (21d)

𝜎 1 𝛼𝛽 𝑛 𝛽 = 0 on 𝜕S 0 . (21e)
Again, the applied load is globally self-equilibrating. The solution of this boundary value problem parametrized by the elongation 𝑈 1 3,3 and the curvatures 𝑈 0 𝛼,33 writes as the linear superposition of each contribution:

𝑢 2 𝛼 = 𝑢 𝑒 3 𝛼 𝑈 1 3,3 + 𝑢 𝜒 𝛽 𝛼 𝑈 0 𝛽,33 + 𝑈 2 𝛼 + 𝑦 𝛽 𝜖 𝛽 𝛼 Θ 2 , ( 22
)
where

𝑢 𝑒 3 𝛼 , 𝑢 𝜒 𝛽
𝛼 are cross-section displacements related to transverse Poisson's effect under pure traction and pure curvatures which are illustrated for a square section in Figure 2. When the section is homogeneous, these correctors have a closed-form expression which is detailed in [START_REF] Zhao | Justification of the asymptotic expansion method for homogeneous isotropic beams by comparison with the Saint-Venant solution[END_REF] for instance. In order to be uniquely defined, the following constraints are applied to all these cross-section displacements: Fig. 3: Warping modes related to pure torsion and pure shear forces for a homogeneous square section with an isotropic material

𝑢 𝛼 = 0 and 𝑦 𝛽 𝜖 𝛽 𝛼 𝑢 𝛼 = 0. ( 23

Longitudinal displacement

The longitudinal displacement 𝑢 2 3 complies with:

W 2 :                𝜎 1 3𝛼, 𝛼 = 0 on S 0 , ( 24a
)
𝜎 1 𝛼3 = 𝐶 𝛼3𝛽3 2𝜀 1 𝛽3 on S 0 , ( 24b
) 2𝜀 1 𝛼3 = 𝑢 2 3, 𝛼 + 𝑦 𝛽 𝜖 𝛽 𝛼 Θ 1 ,3 + 𝑈 1 𝛼,3 on S 0 , ( 24c
)
𝜎 1 𝛼3 𝑛 𝛼 = 0 on 𝜕S 0 . ( 24d
)
The applied load is globally self-equilibrating. The solution of this boundary value problem parametrized by Θ 1 ,3 and 𝑈 1 𝛼,3 writes as the linear superposition of each contribution:

𝑢 2 3 = 𝑢 𝜒 3 3 Θ 1 ,3 + 𝑈 2 3 + 𝑦 𝛼 𝑈 1 𝛼,3 , (25)
where 𝑢 𝜒 3

3 is the torsion warping illustrated in Figure 3a. Indeed, the displacement 𝑢

𝜒 3 3 is
exactly the solution of the Neumann problem for Saint Venant's torsion. Again, this warping is constrained as follows:

𝑢 3 = 0. ( 26
)
Macroscopic constitutive equations From the solution of second order problems, the first order stress may be written as:

         𝜎 1 𝛼𝛽 = 𝜎 𝑒 3 𝛼𝛽 𝑈 1 3,3 + 𝜎 𝜒 1 𝛼𝛽 𝑈 0 1,33 + 𝜎 𝜒 2 𝛼𝛽 𝑈 0 2,33 , (27a)
𝜎 1 𝛼3 = 𝜎 𝜒 3 𝛼3 Θ 1 ,3 , (27b)
𝜎 1 33 = 𝜎 𝑒 3 33 𝑈 1 3,3 + 𝜎 𝜒 1 33 𝑈 0 1,33 + 𝜎 𝜒 2 33 𝑈 0 2,33 . (27c)
Expressing the traction and bending moments leads to the following constitutive equations:

         𝑁 1 3 = 𝐴 3 𝑈 1 3,3 + 𝑆 1 𝑈 0 1,33 + 𝑆 2 𝑈 0 2,33 , (28a)
𝑀 1 1 = 𝑆 * 1 𝑈 1 3,3 + 𝐷 1 𝑈 0 1,33 + 𝐷 12 𝑈 0 2,33 , (28b)
𝑀 1 2 = 𝑆 * 2 𝑈 1 3,3 + 𝐷 * 12 𝑈 0 1,33 + 𝐷 2 𝑈 0 2,33 , (28c)
where:

         𝐴 3 = 𝜎 𝑒 3 33 , 𝑆 1 = 𝜎 𝜒 1 33 , 𝑆 2 = 𝜎 𝜒 2 33 , 𝑆 * 1 = 𝑦 1 𝜎 𝑒 3 33 , 𝐷 1 = 𝑦 1 𝜎 𝜒 1 33 , 𝐷 12 = 𝑦 1 𝜎 𝜒 2 33 , 𝑆 * 2 = 𝑦 2 𝜎 𝑒 3 33 , 𝐷 * 12 = 𝑦 2 𝜎 𝜒 1 33 , 𝐷 2 = 𝑦 2 𝜎 𝜒 2 33 . (29)
The modulus 𝐴 3 is the traction stiffness, 𝐷 1 and 𝐷 2 are the bending stiffnesses. It is possible to prove that 𝑆 * 1 = 𝑆 1 , 𝑆 * 2 = 𝑆 2 and 𝐷 * 12 = 𝐷 12 . The stiffnesses 𝑆 1 and 𝑆 2 are related to the first moments of inertia. Furthermore, there is a position for 𝑂, the origin of the reference frame, such that 𝑆 1 = 𝑆 2 = 0 and rotating the reference frame with respect to 𝒆 -3 , there is an angle such that 𝐷 12 = 0. When the section is homogeneous, this choice of reference frame corresponds to the centroid of the section oriented along one of the principal axis of the second moments of inertia. This is assumed in the following:

         𝑁 1 3 = 𝐴 3 𝑈 1 3,3 , (30a)
𝑀 1 1 = 𝐷 1 𝑈 0 1,33 , (30b)
𝑀 1 2 = 𝐷 2 𝑈 0 2,33 . (30c)
Similarly, the torsion is expressed as function of the macroscopic displacements and the eigenstrain:

𝑀 1 3 = 𝐷 3 Θ 1 ,3 , (31)
where the torsion stiffness is:

𝐷 3 = 𝑦 𝛾 𝜖 𝛾 𝛼 𝜎 𝜒 3 𝛼3 .
(32) Whereas the uncoupling between traction and bending moments may always be satisfied with a proper choice of the reference frame, the uncoupling between torsion is obtained here because of the symmetry assumption (1). This assumption may be released without limiting the approach presented here. Indeed, constitutive equations ( 30) and ( 31) would be simply coupled in such a case and would make the derivation more involved.

ũ - 𝜒 3
3 is the warping used by Benscoter [START_REF] Benscoter | A theory of torsion bending for multicell beams[END_REF].

Third-order problems

Transverse displacement The transverse displacement 𝑢 3 𝛼 is derived through T 3 and loaded by Θ 1 ,33 , 𝑈 2 3,3 and 𝑈 1 𝛼,33 . Noticing that the constitutive equation ( 31) and the torsion equilibrium equation (17c) for 𝑝 = 1, enforce Θ 1 ,33 = 0, it appears that T 3 is exactly the same problem as T 2 , incrementing the orders .

Longitudinal displacement

The longitudinal displacement 𝑢 3 3 complies with W 3 and is loaded by Θ 2 ,3 , 𝑈 1 3,33 and 𝑈 0 𝛼,333 . The applied load is not globally self-equilibrating for each individual loading. From the constitutive equation (30a) and the equilibrium equation (17a)

with 𝑝 = 1 it is deduced that :

𝑈 1 3,33 = - 𝑝 3 𝐴 3 𝐹. ( 33
)
Substituting this in W 3 ensures it is well-posed. Furthermore, it is also convenient to ensure that the load is also self-equilibrating in bending. This is obtained, from the constitutive equations (30b) and (30c) and the equilibrium equation (17b) for 𝑝 = 1:

𝑈 0 1,333 = 𝑉 2 1 -𝜇 1 𝐹 𝐷 1 and 𝑈 0 2,333 = 𝑉 2 2 -𝜇 2 𝐹 𝐷 2 . ( 34
)
Inserting these relations in problem W 3 leads to:

W 3 :                          𝜎 2 3𝛼, 𝛼 + 𝜎 𝜒 1 33 𝑉 2 1 𝐷 1 + 𝜎 𝜒 2 33 𝑉 2 2 𝐷 2 + 𝑓 3 -𝑝 3 𝐴 3 𝜎 𝑒 3 33 -𝜇 1 𝐷 1 𝜎 𝜒 1 33 -𝜇 2 𝐷 2 𝜎 𝜒 2 33 𝐹 = 0 on S 0 , (35a)
𝜎 2 𝛼3 = 𝐶 𝛼3𝛽3 2𝜀 2 𝛽3 on S 0 , (35b)
2𝜀 2 𝛼3 = 𝑢 3 3, 𝛼 + 𝑢 𝜒 1 𝛼 𝑉 2 1 𝐷 1 + 𝑢 𝜒 2 𝛼 𝑉 2 2 𝐷 2 -𝑝 3 𝐴 3 𝑢 𝑒 3 𝛼 + 𝜇 1 𝐷 1 𝑢 𝜒 1 𝛼 + 𝜇 2 𝐷 2 𝑢 𝜒 2 𝛼 𝐹 + 𝑦 𝛽 𝜖 𝛽 𝛼 Θ 2 ,3 + 𝑈 2 𝛼,3 on S 0 , (35c)
𝜎 2 𝛼3 𝑛 𝛼 = 0 on 𝜕S 0 . ( 35d
)
The solution is parametrized by the shear forces 𝑉 2 𝛼 the applied load 𝐹 and higher-order beam displacements. It writes as the linear superposition of each contributions:

𝑢 3 3 = 𝑢 𝑉 𝛼 3 𝑉 2 𝛼 + 𝑢 𝐹 3 𝐹 + 𝑢 𝜒 3 3 Θ 2 ,3 + 𝑦 𝛼 𝑈 2 𝛼,3 + 𝑈 3 3 . with 𝑢 3 = 0 (36)
The longitudinal displacements 𝑢 𝑉 1 3 and 𝑢 𝑉 2 3 are warpings related to shear forces illustrated in Figure (3b) and (3c). Indeed, considering the whole problem W 3 loaded exclusively with shear forces, one can identify the corresponding Neumann problems in Saint Venant's beam theory. Furthermore, the equilibrium equation (35a) considered with only the shear forces loading and integrated on a partial section is actually Jouravskii's Formula [START_REF] Jouravskii | Remarques sur la résistance d'un corps prismatique et d'une pièce composée en bois ou en tôle de fer à une force perpendiculaire à leur longueur[END_REF] which gives a fair estimate of shear stress in thin-walled beams. Finally, 𝑢 𝑉 1 3 and 𝑢 𝑉 2 3 are related to the shear lag phenomenon in homogeneous thin-walled beams: when the beam is not enough slender, the longitudinal stress 𝜎 33 is not exactly the linear distribution derived at leading order (Equation (27c)).

This simplification comes from the symmetry assumption (1) and the scaling chosen for 𝑓 𝛼 such that bending deflection is of order 0.

Macroscopic constitutive equations

From the solution of third order problems, the secondorder stress may be written as:

           𝜎 2 𝛼𝛽 = 𝜎 𝑒 3 𝛼𝛽 𝑈 2 3,3 + 𝜎 𝜒 1 𝛼𝛽 𝑈 1 1,33 + 𝜎 𝜒 2 𝛼𝛽 𝑈 1 2,33 , (37a)
𝜎 2 𝛼3 = 𝜎 𝑉 𝛽 𝛼3 𝑉 2 𝛽 + 𝜎 𝐹 𝛼3 𝐹 + 𝜎 𝜒 3 𝛼3 Θ 2 ,3 , (37b)
𝜎 2 33 = 𝜎 𝑒 3 33 𝑈 2 3,3 + 𝜎 𝜒 1 33 𝑈 1 1,33 + 𝜎 𝜒 2 33 𝑈 1 2,33 (37c)
The traction and bending constitutive equations obtained from this stress distribution are identical to equations ( 30) incrementing the order. The second-order torsion is expressed as function of the macroscopic displacements and the applied load:

𝑀 2 3 = 𝐷 3 Θ 2 ,3 + 𝑦 𝑆 𝛼 𝜖 𝛼𝛽 𝑉 2 𝛽 + 𝑀 𝐹 3 𝐹, (38)
where the torsion induced by the non-uniform distribution of 𝑓 3 in the section is

𝑀 𝐹 3 =
𝑦 𝛾 𝜖 𝛾 𝛼 𝜎 𝐹 𝛼3 and the shear center of the beam is defined as:

𝑦 𝑆 1 = -𝑦 𝛼 𝜖 𝛼𝛽 𝜎 𝑉 2 𝛽3
and

𝑦 𝑆 2 = 𝑦 𝛼 𝜖 𝛼𝛽 𝜎 𝑉 1 𝛽3 . ( 39
)
When the section presents two axes of symmetry, the shear center is in 𝑂 but this is not always true.

Fourth-order and higher-order problems

Transverse displacement The transverse displacement 𝑢 4 𝛼 is derived through T 3 and loaded by 𝑉 2 𝛼 , 𝐹, 𝐹 ,3 , Θ 2 ,33 , 𝑈 3 3,3 and 𝑈 2 𝛼,33 . Taking into account equilibrium equations (17c) and (17d) and constitutive equation ( 38), the problem may be recast only as function of 𝐹 and 𝐹 ,3 , 𝑈 3 3,3 and 𝑈 2 𝛼,33 . This substitution ensures also that the compatibility conditions [START_REF] Corre | A new higher-order elastoplastic beam model for reinforced concrete[END_REF] are satisfied. This leads to:

T 4 :                                            𝜎 3 𝛼𝛽,𝛽 + 𝜎 𝐹 𝛼3 -𝜎 𝜒 3 𝛼3 𝑀 𝐹 3 𝐷 3 𝐹 ,3 + 𝑓 𝛼 -𝜎 𝑉 𝛽 𝛼3 𝑝 𝛽 - 𝜎 𝜒 3 𝛼3 𝐷 3 𝜇 3 -𝑦 𝑆 𝛽 𝜖 𝛽𝛾 𝑝 𝛾 𝐹 = 0 on S 0 , ( 40a
)
𝜎 3 𝛼𝛽 = 𝐶 𝛼𝛽𝛾 𝛿 𝜀 3 𝛿𝛾 + 𝐶 𝛼𝛽33 𝜀 3 33 , 𝜎 3 33 = 𝐶 33𝛼𝛽 𝜀 3 𝛽 𝛼 + 𝐶 3333 𝜀 3 33 on S 0 , ( 40b
)
𝜀 3 𝛼𝛽 = 𝑢 4
( 𝛼,𝛽) on S 0 , (40c)

𝜀 3 33 = 𝑢 𝐹 3 -𝑢 𝜒 3 3 𝑀 𝐹 3 𝐷 3 𝐹 ,3 -𝑢 𝑉 𝛽 3 𝑝 𝛽 + 𝑢 𝜒 3 3 𝐷 3 𝜇 3 -𝑦 𝑆 𝛽 𝜖 𝛽𝛾 𝑝 𝛾 𝐹 +𝑈 3 3,3 + 𝑦 𝛼 𝑈 2 𝛼,33
on S 0 , (40d)

𝜎 3 𝛼𝛽 𝑛 𝛽 = 0 on 𝜕S 0 . ( 40e
)
The solution of this boundary value problem parametrized by the load 𝐹, its first gradient 𝐹 ,3 , the elongation 𝑈 3 3,3 and the curvatures 𝑈 2 𝛼,33 . It writes as:

𝑢 4 𝛼 = 𝑢 𝐹 𝛼 𝐹 + 𝑢 𝐹∇ 𝛼 𝐹 ,3 + 𝑢 𝑒 3 𝛼 𝑈 3 3,3 + 𝑢 𝜒 𝛽 𝛼 𝑈 2 𝛽,33 + 𝑈 4 𝛼 + 𝑦 𝛽 𝜖 𝛽 𝛼 Θ 4 , (41)
The cross-section displacement 𝑢 𝐹 𝛼 is related to the non-uniform distribution in the section of the applied load 𝑓 𝛼 . Similarly, 𝑢 𝐹∇ 𝛼 is related to the variations of 𝐹 and the non-uniform distribution of 𝑓 3 in the section.

Higher orders The induction process may be pursued any higher order. This leads to the derivation of new displacement localizations only related to higher derivatives of 𝐹 relevant for faster variations of 𝐹.

Note that, the use of the asymptotic expansion method is based on the scaling in Equation [START_REF] Argyris | Finite element analysis of two-and three-dimensional elasto-plastic frames-the natural approach[END_REF]. Hence the rescaled coordinates 𝑦 𝑖 have been used in the expression of the auxiliary problems T 𝑝 and W 𝑝 . However, the distinction between the two sets of coordinates is no longer necessary in practice once the section modes are defined. The use of the coordinates 𝑦 𝑖 is therefore dropped in the following and replaced by the use of the coordinates 𝑥 𝑖 .

The higher-order beam model

Families of kinematic enrichment In the asymptotic expansion procedure, three families of kinematic enrichment emerged. First, the rigid motion of the section was carried by the six macroscopic variables 𝑈 𝑝 𝑖 , 𝑈 𝑝 𝛼,3 and Θ 𝑝 3 . They are respectively related to the following displacement modes:

𝒖 - 𝑈 1 = 1 0 0 , 𝒖 - 𝑈 2 = 0 1 0 , 𝒖 - 𝑈 3 = 0 0 1 , 𝒖 - Θ 2 = 0 0 -𝑥 1 , 𝒖 - Θ 1 = 0 0 𝑥 2 , 𝒖 - Θ 3 = -𝑥 2 𝑥 1 0 .
(42) Second, the six correctors related to the six beam resultants were derived:

𝑢 𝑒 3 𝛼 , 𝑢 𝜒 1 𝛼 , 𝑢 𝜒 2 𝛼 , 𝑢 𝜒 3 3 , 𝑢 𝑉 1 3 , 𝑢 𝑉 2 3 .
They are also referred to as Saint Venant's modes [START_REF] Iesan | Saint-Venant's problem for inhomogeneous and anisotropic elastic bodies[END_REF][START_REF] Fatmi | A refined 1D beam theory built on 3D Saint-Venant's solution to compute homogeneous and composite beams[END_REF]. This collection of 12 modes (rigid and Saint-Venant) is denoted B SV . Third, a family of modes related to the body force and its longitudinal variations was obtained:

𝒖 - 𝐹 , 𝒖 - 𝐹∇ , 𝒖 - 𝐹∇ 2 .
.., which is denoted

B 𝑓 .
This suggests gathering all these modes in the following kinematic approximation for the 3D displacement:

𝒖 -= 𝑛 ∑ 𝑘=1 𝒖 - 𝑘 (𝑥 𝛼 ) 𝑋 𝑘 (𝑥 3 ) (43)
where 𝑛 is the number of modes and 𝑋 𝑘 (𝑥 3 ) are longitudinal unknown fields. It is demonstrated in [START_REF] Miara | A Galerkin spectral approximation in linearized beam theory[END_REF] that the modes generated for a single prescribed body load distribution 𝑓 𝑖 in the cross-section is a Hilbert basis of the solution space for 𝒖 -.

Load superposition With proper scaling, it is possible to consider other loadings in a separate form between longitudinal coordinate and cross-section coordinates such as surface traction on the lateral boundary of the cross-section as well as prescribed eigenstrains. Furthermore, whereas most loadings in engineering practice comply with the separate form requested in equation ( 6) the latter may be too constraining in some specific cases. The problem being linear, it is possible to superpose loads. In that case, the collection of modes 𝒖 -𝑖 may become linearly dependent and is orthogonalized.

In practice, it appears that decomposing a cross-section load distribution into several distinct loadings improves the quality of the solution. For instance, in Equation ( 6), one may modulate 𝑓 3 (𝑥 𝜂 ) with 𝐹 (𝑥 3 ) separately from 𝑓 𝛼 (𝑥 𝜂 ) with another function 𝐺 (𝑥 3 ). However, this doubles the number of kinematic degrees of freedom related to the applied load (third family). The best trade off between the load decomposition and the accuracy of the solution

has not yet been explored in details.

From the traction, bending and torsion constitutive equations ( 30) and [START_REF] Vogelius | On a Dimensional Reduction Method II. Some Approximation-Theoretic Results[END_REF], 𝑢

Minimum potential energy principle

The higher-order beam model is obtained, inserting the approximated kinematics (43) in the minimum potential energy principle of the 3D problem [START_REF] Amrouche | On the characterizations of matrix fields as linearized strain tensor fields[END_REF]. Indeed, the variable separation between longitudinal and cross-section coordinates allows the separation of the corresponding integrals. Since all cross-section modes are known there remains only a boundary value problem related to the longitudinal coordinate.

The corresponding equilibrium equations write as:

∀𝑙, 𝑋 𝑘 𝐴 𝑘𝑙 + 𝑋 𝑘 ,3 (𝐵 𝑘𝑙 -𝐵 𝑙𝑘 ) -𝑋 𝑘 ,33 𝐶 𝑘𝑙 -𝐹 𝑙 = 0 ( 44
)
where: (45)

𝐴 𝑘𝑙 = 𝑢 𝑘 𝑖,
Boundary conditions In the 3D problem (2), the beam extremities were assumed fully clamped. Considering the approximated kinematics (43), this is achieved enforcing 𝑋 𝑖 = 0 at extremities. Other boundary conditions may be applied. Indeed, the first six modes corresponds to the rigid motion of the section. Restraining only these degrees of freedom is actually the boundary condition classically used in structural mechanics: warping and transverse displacements are let free.

Applications

Two case study are briefly presented. The first case study illustrates the advantage of introducing displacement modes related to the distribution in the cross-section of the applied load. The second case study was extended to eigenstrains and illustrates the possibility of describing more accurately the effect of prestress in a beam. In both cases, the auxiliary problems [START_REF] Corre | The Asymptotic Expansion Load Decomposition elastoplastic beam model[END_REF] and ( 13) generating the kinematic enrichment were solved with 2D quadratic Lagrange finite elements. Two different strategies were used for the higher-order beam solution.

Cantilever box girder beam under a concentrated load

A cantilever beam of span 10 m with a box section 3 m wide is investigated (Figure 4). The width and the height of the box are 1 m. The thickness of the walls is 5 cm. The beam is fully clamped at its first extremity and loaded by a concentrated force in the middle of the section at the other extremity (Figure 4).

In addition to the six rigid modes (42) and the six Saint Venant modes (B SV ), 3 warping and 3 transverse displacement modes (Figure 5) are included in the kinematics (B 𝑓 ). These modes clearly illustrate the local influence of the distribution in the section of the concentrated load.

Because they lead to ordinary differential equations, beam models are suited for closedform solutions. The homogeneous part of equations (44) have solutions of the form 𝑒 𝑘 𝑥 3 which may be found explicitly solving a quadratic eigenvalue problem [START_REF] Tisseur | The Quadratic Eigenvalue Problem[END_REF]. This technique was used in [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF] in order to find the solution of the beam loaded at its extremity. It presents the advantage to yield a solution extremely rapidly. However, because it may include sharp decaying exponential functions, its implementation requires specific treatment. Furthermore, it is restricted to linear problems with specific longitudinal distribution of the applied load.

Figure 4 shows the reconstructed beam solution. The local settlement related to the application of the concentrated force is very well captured.

Prestressed cable in a cantilever beam

The higher-order beam model is further illustrated with a cantilever concrete beam prestressed with a steel cable. The section of the beam is represented in Figure 6. The beam is 20 m long.

The concrete domain is denoted by Ω 𝑐 and the steel domain by Ω 𝑠 . A constant eigenstrain The higher-order beam model was discretized with NURBS shape functions. Contrary to the closed-form solution previously mentioned, this allows to apply various load distributions.

𝜺 * = 𝜀 * 33 𝒆 -3 ⊗ 𝒆 -3 is applied in Ω 𝑠 ,
Additionally, NURBS shape functions are known to mitigate locking phenomena. This was confirmed in a detailed study of the present implementation in [START_REF] Corre | Higher-order beam model with eigenstrains: theory and illustrations[END_REF].

The reconstructed solution is presented in Figure 8. As expected, the prestressed cable compresses and raises the beam. The higher-order beam model captures transverse displacements illustrated by the lowering of the edges of the table. It also captures warping displacements illustrated by the punching effect of the cable which can be observed at the end of the beam. A 3D reference solution was also computed by means of Code_Aster on a very similar mesh. The deformed beam is also presented in Figure 8 and the computation data is summarized in Table 1. It appears that, the model presented here shows very satisfying results with a significantly reduced computational cost. Indeed, only very few additional modes (3 to 5) and their corresponding beam DOF were required for capturing fairly well the applied eigenstrain. These time performances can still largely be improved by more advanced implementations and numerical techniques.

However, detailed analysis shows that in the vicinity of the boundaries, the quality of the solution is not as satisfying as in the bulk of the beam (see details in [START_REF] Corre | Higher-order beam model with eigenstrains: theory and illustrations[END_REF]). This is not surprising since higher-order convergence results has only been obtained in configurations were boundary layers related to the beam extremities are not present. In order to better describe the mechanical behavior close to boundary conditions, the approach developed independently by Strains Engineering consists in the introduction of new modes specific to these boundary conditions in the reduced basis.

The Asymptotic Expansion Load Decomposition elasto-plastic beam model

Introducing elasto-plastic behavior is more challenging. The inherent non-linearity of plasticity and the incremental nature of plastic analysis makes the definition of a relevant kinematics more difficult. Two main approaches are followed when solving an elasto-plastic beam problem: 1D elasto-plastic beam model based on a priori cross-section analysis or 3D elasto-plastic beam models based on a 3D beam kinematics.

The first natural approach is to express the plastic flow in terms of generalized beam variables and to solve an elasto-plastic 1D problem. This requires the elasto-plastic analysis of the cross-section for pure or combined generalized stresses and the derivation of the corresponding yield surface. The cross-section analysis may be incremental or based on limit analysis but assumes a uniform distribution of generalized stresses in the longitudinal direction: normal force, shear forces, bending moments and torque. A key difficulty is the derivation of a yield surface directly function of the beam generalized stress taking into account correctly their possible interactions as well as hardening. There were recent improvements in this direction, approximating the yield surface with facets or ellipsoids [START_REF] Bleyer | Yield surface approximation for lower and upper bound yield design of 3D composite frame structures[END_REF].

Once the yield surface is defined, there remains to compute the elasto-plastic response of the beam, either with closed form solutions [START_REF] Štok | Analytical solutions in elasto-plastic bending of beams with rectangular cross section[END_REF], limit analysis [START_REF] Olsen | Rigid plastic analysis of plane frame structures[END_REF] or by means of finite element approximations [START_REF] Argyris | Finite element analysis of two-and three-dimensional elasto-plastic frames-the natural approach[END_REF][START_REF] Papadrakakis | A computationally efficient method for the limit elasto plastic analysis of space frames[END_REF]. This approach has the advantage to present fast computation time, since only a 1D elasto-plastic problem needs to be solved. However, its accuracy remains limited by the beam theory assumptions. First, it cannot handle local phenomena related to the distribution of the applied load as well as to the boundary conditions. Second, it provides only an averaged description of the actual stress in the beam.

In order to improve the accuracy of the beam model, the second classical approach consists in setting a beam kinematics expressing the 3D displacement field in a separate form between the transverse coordinates and the longitudinal coordinate. This kinematics may be defined a priori or may evolve during the incremental procedure. For a fixed increment of the generalized displacements, the corresponding 3D stress is computed and the yield criterion is expressed locally. A local algorithm such as the radial return is processed on the whole body to compute the local plastic state of the beam. This locally admissible stress state is integrated on each cross-section yielding the corresponding longitudinal distribution of the beam generalized stresses. Finally, the beam global equilibrium is ensured with a standard Newton-Raphson procedure. The main difficulty lies again in the definition of a relevant kinematics able to describe the displacement related to plastic flow.

Most approaches where the kinematics is fixed a priori rely on the ones already discussed in linear elasticity such as Euler-Bernoulli, Timoshenko kinematics or even Saint-Venant solution, eventually with non-linear geometric corrections. Once the kinematics is defined, there remains to choose the number of integration points in the cross-section in order to compute precisely the local plastic flow. Multiplying integration points improves the accuracy of the results at the price of a higher computation time of the cross-sections integrals. This is the spirit of multi-fiber beam models (see for instance [START_REF] Mazars | Using multifiber beams to account for shear and torsion. Applications to concrete structural elements[END_REF]). Another direction is to enrich arbitrarily the section kinematics with degrees of freedom not necessarily related to classical cross-section displacements. This concept was formalized extensively by Carrera et al [START_REF] Carrera | Beam Structures[END_REF] and co-workers.

Because plastic flow may not be easily known a priori a natural improvement of the preceding methods is to update the beam kinematics during the load increments. This is the direction followed here using the asymptotic expansion.

3.1 Adaptation of the higher-order beam model to the elasto-plastic behavior

The linear higher-order beam model presented in the previous section is extended to elastoplasticity in the small strains framework. Hence, the boundary value problem expressed in Equation ( 2) is now considered with a 𝐽 2 elasto-plastic constitutive law and will be solved classically following a Newton-Raphson incremental scheme.

As for linear elasticity, a 3D approximated kinematics is formulated in a separated form in order to perform the dimension reduction. The first collection of modes to take into account is the basis of Saint Venant modes B SV as well as the basis of force modes B 𝑓 related to the applied load described previously. Nevertheless, this kinematics is not sufficient for describing accurately the possibly discontinuous plastic flow which may occur in the cross section. Therefore, the plastic strain computed at a given iteration of the Newton-Raphson procedure is taken into account for enriching the kinematics of the following iteration. This is possible using the procedure described in the previous section for a fixed plastic strain 𝜺 p distribution in the cross-section. Considering now the whole beam, 𝜺 p is not longitudinally uniform. Hence, several chosen cross-sections may be used for taking snapshots of the plastic strain in order to sufficiently enrich the kinematics of the model. These new plastic modes are computed and added to the kinematics on the fly. The basis of modes specific to a crosssection plastic strain distribution is denoted by B 𝜀 . Finally, the kinematics of the model is evolving during the Newton-Raphson procedure and is the union of the basis B SV , B 𝑓 and B 𝜀 . This union of basis is orthonormalized to form the total basis B with a total number of modes 𝑛 mod .

This approach presents two major advantages. First, it does not require additional elastoplastic computations in the cross-section. Second, the number of beam degrees of freedom remains very limited (about 20) thanks to the sparsity of the kinematics. From the optimality result proved in [START_REF] Miara | A Galerkin spectral approximation in linearized beam theory[END_REF], this approach is expected to be more efficient than arbitrary kinematic refinements.

Note that, contrary to Nonuniform Transformation Field Analysis [START_REF] Michel | Nonuniform transformation field analysis[END_REF][START_REF] Michel | Computational analysis of nonlinear composite structures using the nonuniform transformation field analysis[END_REF][START_REF] Fritzen | Reduced basis hybrid computational homogenization based on a mixed incremental formulation[END_REF] where a basis of plastic strains is introduced with the corresponding plastic multipliers, in the present approach, displacement plastic modes are added to the total 3D displacement approximation and plasticity is processed at each integration point of the 3D body.

The elasto-plastic algorithm

The implementation of the general framework introduced in the previous section is now detailed. This requires first the definition of the numerical approximation of the 3D body.

Then, the incremental resolution of the elasto-plastic problem is adapted so that processing the local constitutive equations remains standard whereas the global equilibrium iterations are performed with the reduced basis.

Numerical approximation of the higher-order beam model

The approach suggested in the previous section requires the definition of a 3D mesh of the beam composed of cross-sections meshes positioned along the longitudinal direction (Figure 9). Indeed, these cross-sections will be the domain of integration of the constitutive law.

Longitudinal discretization A longitudinal discretization of the beam is defined for the generalized displacements 𝑋 𝑚 (𝑥 3 ) introduced in Equation ( 43). The same collection of NURBS basis functions for each 𝑋 𝑚 (𝑥 3 ) is chosen:

𝑋 𝑚 (𝑥 3 ) = 𝑛 NURBS ∑ 𝑖=1 𝑁 𝑖 (𝑥 3 ) 𝑋 𝑚,𝑖 (46)
Fig. 9: Discretization of a square beam where 𝑁 𝑖 (𝑥 3 ) are the NURBS interpolation function and 𝑋 𝑚,𝑖 are the corresponding degrees of freedom. A set of 𝑁 𝑠 > 𝑛 NURBS longitudinal integration points is also defined for the integration of the interpolation functions. It is natural to place the cross-section meshes at the positions of these longitudinal integration points.

Cross-section discretization

The cross-section mesh used for the computation of the modes is the same as the one used in Section 2.3.2: the modes are computed by means of quadratic Lagrange triangle elements:

ũ - 𝑚 (𝑥 𝛼 ) = 𝑛 sec ∑ 𝑗=1 𝑳 - 𝑗 (𝑥 𝛼 ) 𝑢 𝑚, 𝑗 (47)
where 𝑛 sec is the number of nodes in the section, 𝑳 -𝑗 (𝑥 𝛼 ) are Lagrange interpolation functions and 𝑢 𝑚, 𝑗 are the nodal values of the displacement.

The local state variables, Δ𝜺, Δ𝜺 p and Δ𝝈 are computed in the 𝑁 𝑠 cross-sections meshes at the Gauss points of the quadratic triangle elements.

Plastic-mode cross-section During the computation, sections where a plastic strain is not zero are collected. All cross-sections where plasticity occurs could be used for the computation of the plastic modes: for each plastic strain distribution in each cross section, one or several modes could be computed. However, it would excessively increase the number of generalized displacement degrees of freedom 𝑛 dof and also increase the computation time dedicated to the corresponding modes. In order to limit the number of plastic modes to a few, only one cross-section called plastic-mode cross-section is chosen for taking snapshots of the plastic strain distribution. As a first approach, this choice is based on an educated guess.

Adaptation of the Newton-Raphson procedure

Formulation of the tangent stiffness of the beam model For a given basis of 𝑛 mod displacement modes, the numerical approximation of the total displacement may be written as follows:

𝒖 -(𝒙 -) = 𝑛 mod ∑ 𝑚=1 𝑛 NURBS ∑ 𝑖=1 𝑛 sec ∑ 𝑗=1 𝑳 - 𝑗 (𝑥 𝛼 )𝑁 𝑖 (𝑥 3 ) 𝑢 𝑚, 𝑗 𝑋 𝑚,𝑖 = 𝑛 NURBS ∑ 𝑖=1 𝑛 sec ∑ 𝑗=1 𝑳 - 𝑗 (𝑥 𝛼 )𝑁 𝑖 (𝑥 3 )𝑢 𝑖, 𝑗 (48)
where the 𝑛 3D = 𝑛 NURBS × 𝑛 sec local displacement degrees of freedom are:

𝑢 𝑖, 𝑗 = 𝑛 mod ∑ 𝑚=1 𝑢 𝑚, 𝑗 𝑋 𝑚,𝑖 . (49)
Considering the second form of (48), it appears that the 3D strain 𝜀(𝒖 -) and consequently the local constitutive equation may be directly computed from 𝑢 𝑖, 𝑗 without the need to specify the basis B. Hence, provided the finite element solution is stored as the collection of the 𝑢 𝑖, 𝑗 , the local integration of the constitutive equations as well as the computation of the local elasto-plastic tangent stiffness remain unchanged compared to a 3D formulation.

For a fixed basis B, injecting the numerical approximation of the kinematics (48) into the principle of virtual work corresponding the elasto-plastic 3D problem (2) leads to the expression of the residual expressed in terms of the increment of the 𝑛 dof = 𝑛 mod × 𝑛 NURBS generalized displacement degrees of freedom Δ𝑋 𝑚,𝑖 and the corresponding test degrees of freedom. The standard Newton-Raphson procedure is used in order to cancel this residual which leads to the following formal reduced equilibrium equation:

[𝑲 ep ] B 𝛿 𝑿 -B = {𝑹} B , (50)
where 2 is the global tangent stiffness matrix of the beam model expressed in basis B, computed with the local elasto-plastic tangent stiffness, 𝛿 𝑿 -B is the finite element vector of the generalized degrees of freedom 𝛿𝑋 𝑚,𝑖 and {𝑹} B ∈ R 𝑛 dof is the residual vector.

[𝑲 ep ] B ∈ (R 𝑛 dof )
Because, the tangent stiffness as well as the residual both depend on the choice of the basis B, they need to be updated each time the basis is changed.

Significant computational time is gained because 𝑛 dof which sets the size of the tangent stiffness is much smaller than the rather large 𝑛 3D which is required for a sufficiently detailed description of the fields in the cross-section.

Description of the numerical procedure The basis of modes B is first initialized and is composed B SV and 𝑛 fAE modes associated to the applied load B 𝑓 . While increments do not generate plastic flow, the basis B remains unchanged. The global tangent stiffness corresponds to the elastic one and each increment is solved in one iteration.

Let assume that increment 𝑛 generates a non-vanishing plastic strain Δ𝜺 p 𝑛 at the first iteration 𝑘 = 1. In this case, before starting the second iteration, the basis of modes is enriched and orthonormalized with 𝑛 pAE modes computed from the plastic strain distribution Δ𝜺 p 𝑛 observed in the plastic-mode cross-section, following [START_REF] Corre | Higher-order beam model with eigenstrains: theory and illustrations[END_REF]. This also requires the update of the residual.

It has been noticed from experience that plastic modes computed at subsequent iterations of the increment were very similar. Therefore the basis B used at iteration 𝑘 = 2 is kept until the convergence of the increment is reached. However, the converged plastic strain of the increment Δ𝜺 p 𝑛 may have changed. Hence, at the first iteration 𝑘 = 1 of the following increment 𝑛 + 1, the basis B is updated, replacing only plastic modes with new ones. Again, at the second iteration 𝑘 = 2 the basis is updated and then remains fixed until the convergence of the increment. This choice of updating the plastic modes only at the first two iterations of the increment remains valid as long as the load increments are not too important.

Application to a cantilever beam

To illustrate the efficiency of the model presented, a steel beam clamped at one end and loaded on its free end is investigated. The beam is a wide flange beam HE600M. This cross-section is class 1 in Eurocode 3, meaning that the beam reaches its limit of elasticity with no risk of local buckling. The geometry of the 6 m long beam is detailed in Figure 10. A load is applied with eccentricity at the top edge of the free end of the beam, as represented in Comparison with the 3D reference solution The beam solution is 80 times faster to compute.

The corresponding deformed structure is presented in Figure 14. The torque due to the eccentricity of the load induces a longitudinal rotation of the cross-section, and the transverse part of the load induces a bending of the beam. the 3D solution induces a difference between the two residual displacements observed at 𝑡 20 . For the solution 𝑆 1D considered here, the relative distance between the two solutions at 𝑡 20 is of 8.51%. This relative distance originates from the gap observed at 𝑡 10 which is maintained between 𝑡 11 and 𝑡 20 since the unloading is elastic. A residual displacement closer to the 3D reference could be obtained with more refined longitudinal meshes and with higher interpolation order for the interpolation functions. But in regards with the very low computational time offered by the beam solution, the results obtained can be considered satisfying for engineering applications.

In order to assess the accuracy of the beam solution when compared to the reference solution the six components of the plastic strain computed by the solution 𝑆 1D and by the reference solution at 𝑥 3 = 0, 5 m for the given displacement 𝑢 𝐴 = 0.8 m of point A, are presented in Figure 16. All the variables presented for a fixed displacement of point A are obtained by linear interpolations between the increments defined in Section 3.3 and figured by markers in Figure 15. The axial plastic strain presented in Figures 16a to 16c shows that nearly all parts of the section have reached the elastic limit. As expected, the eccentricity of the load on the free extremity of the beam creates a slightly uneven progression of plasticity in the section.

Therefore the highest values of the plastic strain components are observed at the top left and the bottom left of the section where the absolute values of the stresses are largest. The plastic strain computed by 𝑆 1D is slightly lower than the plastic strain computed by 𝑆 ref for each component. A late detection of plasticity due to the longitudinal refinement is suspected to originate this phenomenon.

The non-axial components presented in Figures 16d to 16f seem less satisfying but their amplitude is about 10 times lower than the axial components. Thus, the beam solution presented here shows satisfying results with a good comparison with the 3D solution.

Discussion and outlooks

Several applications of the AELD higher-order elasto-plastic beam model were presented.

The model is based on an enrichment of the reduced kinematics updated on the fly and does not need any a priori knowledge on the solution of the problem.

Both in linear elasticity framework and in 𝐽 2 elastoplasticity, the beam model required a computation about 80 times shorter than for the reference 3D solution computed on Code_Aster with comparable meshes. These gains in computation time are already very interesting and could be even larger by means of parallelization, in particular for the computation of the tangent stiffness matrix. The numerical differences between the beam solution and the reference solution for the examples presented are low, and could still be lowered. Indeed, the computation of the displacement modes is optimal only far from boundary conditions. At the clamped extremity the local stress varies rapidly. The model could therefore be improved if displacement modes specific to the boundary conditions were computed and added to the kinematics of the model. Furthermore, the choice of the plastic-mode cross-section was arbitrary and could be automated based for instance on the location of the first occurence of plasticity in the beam.

We conclude that the present beam model offers already interesting perspectives for engineering applications on the cloud. Nevertheless, it is meant to be validated on more complex case-studies.

Several extensions of the linear elastic AELD model are straightforward. The low and medium frequency dynamic behavior is easy to derive including inertia contributions from the assumed kinematics [START_REF] Taig | An unconstrained dynamic approach for the Generalised Beam Theory[END_REF]. Buckling analysis may also be adapted [START_REF] Taig | GBT pre-buckling and buckling analyses of thinwalled members under axial and transverse loads[END_REF], as the eigenstrains enrichment is a suitable basis for representing the pre-stress in the beam and deriving the geometric stiffness operator. Similarly, other weak couplings (thermal, humidity, creep etc.) or non standard physics [START_REF] Lurie | On the formulation of elastic and electroelastic gradient beam theories[END_REF] may be incorporated in the model. . Furthermore, the elastoplastic beam model has been presented here with an isotropic material and a 𝐽 2 yield criterion.

Its adaption to different yield criteria is easy and the extension to more complex materials like reinforced concrete was achieved with a Rankine yield criterion in [START_REF] Corre | A new higher-order elastoplastic beam model for reinforced concrete[END_REF].

Finally, as long as strains remains small, the kinematic enrichment provided by the asymptotic expansion remains relevant. Hence, the authors are considering a co-rotational formulation of the AELD beam model in order to introduce geometric non-linearities related to large displacements and rotations.

  , Einstein's implicite summation convention on repeated indices is used. Note that, Greek indices denotes transverse dimensions: 𝛼, 𝛽... = 1, 2 and Latin indices denotes all three dimensions: 𝑖, 𝑗 ... = 1, 2, 3.

2 𝛼Fig. 2 : 3 (

 223 Fig. 2: Transverse modes related to pure traction and pure curvatures for a homogeneous square section with an isotropic material

  related to 𝑁 3 , 𝑀 𝛽 and 𝑀 3 .

Fig. 4 :Fig. 5 :

 45 Fig.4: Cantilever box girder beam under applied load[START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF]

Fig. 6 :

 6 Fig. 6: Section of the prestressed beam

4 Fig. 7 :

 47 Fig. 7: The 5 eigenstrain modes related to the prestress in the steel cable used in case study 2.

Fig. 8 :

 8 Fig.8: Deformed structure under pre-tensioned steel cable (amplification factor = 100),[START_REF] Corre | Higher-order beam model with eigenstrains: theory and illustrations[END_REF]

Figures 10 and 11 .Fig. 10 : 2 Fig. 11 :Fig. 12 : 3 Fig. 13 :

 10111021112313 Fig. 10: Dimensions (mm) of the HE600M section, mesh and applied load

Fig. 14 :Fig. 15 :

 1415 Fig. 14: Deformed shape of the beam after the 10 load increments (solution 𝑆 1D )

Fig. 16 :

 16 Fig. 16: Plastic strain in solutions 𝑆 1D and 𝑆 ref close to the clamped extremity at 𝑥 3 = 0.5 m for 𝑢 𝐴 = 0.8 m

  𝛼 𝐶 𝑖 𝛼 𝑗 𝛽 𝑢 𝑙 𝑗,𝛽 , 𝐵 𝑘𝑙 = 𝑢 𝑘 𝑖 𝐶 𝑖3 𝑗 𝛽 𝑢 𝑙 𝑗,𝛽 , 𝐶 𝑘𝑙 = 𝑢 𝑘 𝑖 𝐶 𝑖3 𝑗3 𝑢 𝑙 𝑗 , and 𝐹 𝑘 = 𝑢 𝑘 𝑖 𝑓 𝑖 .

Table 1 :

 1 Main features of the 3D solution and solution 𝑆

	type of elements number of elements CPU computation time	reference model 15-nodes pentahedron 6-nodes triangle + 41 longitudinal knots higher-order beam model 99680 1788 + 1 1805 s 24 s