Grégoire Corre

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

Karam Sab

Mohammed Khalil Ferradi

Xavier Cespedes

Mohammed Khalil Ferradi

Higher-order beam model with eigenstrains: theory and illustrations

Keywords: Beam Model, Asymptotic Expansion, Eigenstrains, NURBS MSC (2010) 74K10, 74S05, 74B10

come

Introduction

Structural analysis must take into account eigenstrains such as thermal load, prestress or weakly coupled swelling phenomena. Furthermore, beam models are certainly the most widely used models by structural engineers because of their simplicity. This motivates the introduction of these prescribed strains into beam models. However, the kinematics of beam models is often based on very simple assumptions which are often unable to describe precisely local effects induced by these strains prescribed arbitrarily in the section. In the absence of richer kinematics, shell models or full 3D computations are necessary. Those computations being time-consuming, the simplicity and the great time efficiency of beam models motivates the development of adequately refined higher-order beam models.

Beam models have a long history and a vast variety of contributions may be found in the literature. The first beam models were based on ad-hoc assumptions on the 3D fields which motivated the denomination axiomatic. Most of the time axiomatic models rely on an educated guess on the 3D displacement field in a separated form between the longitudinal coordinate and the in-section coordinates. Then, straightforward application of the minimum of potential energy leads to 1D boundary value problems corresponding to the beam model. The Euler-Bernoulli beam model was the first suggestion. In this model, it is assumed that the cross-section of the beam is rigid in its own plane and that it remains orthogonal to the neutral axis of the beam. Hence, this model neglects the transverse shear strain and suffers from a kinematic contradiction: since the section is not allowed to deform in its plane, transverse Poisson's effect is precluded. The Timoshenko beam model allows an independent rotation of the section with respect to the neutral axis in order to take the transverse shear strain into account [START_REF] Timoshenko | On the transverse vibrations of bars of uniform cross-section[END_REF]. However, since the section rotates rigidly, the shear strain is uniform in the section which does not satisfy the free lateral boundary. Furthermore, the transverse Poisson's effect is still restrained. These contradictions, suggested to push further the kinematic enrichment in a polynomial form as illustrated in [START_REF] Carrera | Beam Structures[END_REF][START_REF] Giunta | A thermal stress finite element analysis of beam structures by hierarchical modelling[END_REF]. However, having a correct approximation requires a high number of kinematic degrees of freedom and raises the question of the sparsity of the approximation as will be recalled below.

It turns out that, at leading order in the slenderness, these kinematic contradictions were resolved quite early thanks to Saint Venant solution [START_REF] De Saint | De la torsion des prismes: avec des considérations sur leur flexion ainsi que sur l'équilibre des solides élastiques en flexion[END_REF]. Indeed, this solution was originally derived for an elastic beam with a homogeneous and isotropic section loaded at the extremities in a weak sense. Starting from static considerations, a full 3D solution was obtained where the 6 classical generalized stress vary linearly along the beam and the section is free to deform in its plane as well as out of its plane. More precisely, the 3D displacement field appears as the superposition of the classical rigid motion of the section and of additional displacements related to the generalized stress which correct the over-constrained rigid motion of the section. As consequence, the 3D stress is better approximated than in the preceding approaches and the traction, bending, torsional and shear force stiffnesses are correctly evaluated.

Whereas this solution is perfectly relevant for rather compact and homogeneous sections, the Saint Venant solution is not accurate enough when considering more general sections such as anisotropic or heterogeneous materials, open and closed thin-walled sections. A classical illustration involves torsion warping. Indeed, the latter is assumed uniform along the beam in the Saint Venant solution. Hence, clamped boundary conditions may not be strongly satisfied. In practical situations, clamping a thin-walled beam generates a kinematic frustration which propagates far from the extremity and needs correct estimation. A first successful extension of Saint Venant solution was made by Vlasov [START_REF] Vlasov | Thin-walled elastic beams[END_REF] who included the torsional warping as an independent static degree of freedom (so called bi-moment) followed by Benscoter [START_REF] Benscoter | A theory of torsion bending for multicell beams[END_REF] who further assumed that the corresponding kinematic degree of freedom is also independent. Both works rely on the torsional warping correction included in the Saint Venant solution.

The efficiency of these models to capture end effects encouraged further refinement in more general configurations. A noticeable contribution is from Iesan [START_REF] Iesan | Saint-Venant's problem for inhomogeneous and anisotropic elastic bodies[END_REF] who extended the Saint Venant solution to the case of fully anisotropic and heterogeneous section. In addition the beam was loaded with body forces and surface tractions in a separated form between in-section coordinates and the longitudinal coordinate. The longitudinal variation of the loads was assumed polynomial of a fixed order and the corresponding solutions for each order is obtained by a recursive process.

Almost identical results were obtained from the formal asymptotic expansion. Indeed, this method is based on a scaling of the original 3D problem so that it depends explicitly on a small parameter. In the case of beams, it is the ratio between the typical size of the cross section and the length of the beam. Then the solution is assumed following an asymptotic expansion with respect to the small parameter and inserted in the 3D equations. A collection of embedded in-section and longitudinal problems is obtained which is solved by induction. The monograph from Trabucho & Viano [START_REF] Trabucho | Mathematical modelling of rods[END_REF] presents the method and links the lowest orders of the asymptotic expansion with the Saint Venant solution as well as Vlasov beam model. Most of recent contributions questions the correct derivation of boundary conditions especially at higher order [START_REF] Buannic | Higher-order effective modeling of periodic heterogeneous beams. I. Asymptotic expansion method[END_REF][START_REF] Buannic | Higher-order effective modeling of periodic heterogeneous beams. II. Derivation of the proper boundary conditions for the interior asymptotic solution[END_REF][START_REF] Kim | An asymptotic analysis of composite beams with kinematically corrected end effects[END_REF][START_REF] Kim | On the asymptotic boundary conditions of an anisotropic beam via virtual work principle[END_REF]. Another family of enriched beam model was derived following the so-called "Variational Asymptotic Method" [START_REF] Hodges | Nonlinear Composite Beam Theory[END_REF][START_REF] Yu | Elasticity Solutions Versus Asymptotic Sectional Analysis of Homogeneous, Isotropic, Prismatic Beams[END_REF][START_REF] Yu | Variational asymptotic beam sectional analysis -An updated version[END_REF][START_REF] Yu | Validation of the Variational Asymptotic Beam Sectional Analysis[END_REF]. This approach is very similar to formal asymptotic expansion and most of its developments were made assuming large displacements and rotations of the section which makes higher-order developments much more involved.

In addition to the difficulties regarding boundary conditions, the main limitations of the asymptotic expansion approach are the very high regularity of the applied load which is required when going higher order and the embedded structure of the sequence of 1D problems to be solved. This makes the classical approach impractical for engineering applications. Remarkably, the same difficulty arises with the extended Saint Venant solution from Iesan [START_REF] Iesan | Saint-Venant's problem for inhomogeneous and anisotropic elastic bodies[END_REF] since only polynomial loadings are solution. For instance, concentrated loads commonly used in practice do not satisfy such regularity.

A solution to overcome this difficulty is to consider the whole family of in-section displacements or "modes" generated from formal asymptotic expansion as kinematic enrichment carried by independent generalized beam displacements. Indeed, the application of the minimum of the potential energy does not restrict severely the regularity of the longitudinal distribution of the applied load. The very first illustration of this approach is the model from Benscoter [START_REF] Benscoter | A theory of torsion bending for multicell beams[END_REF] which treats the torsional warping as an independent degree of freedom. In the case of a homogeneous and isotropic beam this idea was fully generalized by Miara & Trabucho [22] (also detailed in [START_REF] Trabucho | Mathematical modelling of rods[END_REF]) and so called "Galerkin spectral approximation". This work is based on the seminal idea from Vogelius & Babuska [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]. Two noticeable observations were made. 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. 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.

Practical implementation of this concept were suggested by El Fatmi [START_REF] Fatmi | A refined 1D beam theory built on 3D Saint-Venant's solution to compute homogeneous and composite beams[END_REF] Lahmar et al. [START_REF] Lahmar | Thermo-mechanical analysis of composite beams[END_REF] and Ferradi et al. [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF]. In [START_REF] Fatmi | A refined 1D beam theory built on 3D Saint-Venant's solution to compute homogeneous and composite beams[END_REF], the enrichment was limited to the contribution of the Saint Venant modes and transverse loadings which are uniform in the longitudinal direction. However, longitudinal or higher-order enrichments were discarded and there remained an educated guess for enriching further the beam model. This approach was also extended to uniform thermal loads in [START_REF] Lahmar | Thermo-mechanical analysis of composite beams[END_REF]. In [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF], the enrichment related to any kind of load was introduced up to an arbitrary order. A closed-form solution of the higher-order beam model was derived and comparisons with full 3D calculations were performed. Even for a concentrated load arbitrarily located in the section, the approximated solution yielded surprisingly good results.

In the present paper, the same approach is applied to eigenstrains and illustrated with two cases study. Furthermore, the higher order beam model is implemented with NURBS finite elements in order to allow longitudinal variations of the applied load. There are fewer contributions related to applied eigenstrain in beam theories. In addition to those previously mentioned, the general case of periodic beams as well as thin walled beam was investigated by Kolpakov [START_REF] Kolpakov | Problem of the theory of beams with initial stresses[END_REF][START_REF] Kolpakov | Application of homogenization method to justification of 1-D model for beam of periodic structure having initial stresses[END_REF][START_REF] Kolpakov | Stressed Composite Structures: Homogenized Models for Thin-Walled Nonhomogeneous Structures with Initial Stresses[END_REF]. However, the formal asymptotic expansion was not carried out up to an arbitrary order.

x 2

x 1 x 3 L h ∂Ω t Ω S +

S -

Fig. 1 The beam configuration

This paper is also an opportunity to redevelop the formal asymptotic expansion procedure and 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 is organized as follows. First the formal asymptotic expansion procedure is applied to a beam with a prescribed eigenstrain and the higher-order beam theory is derived and solved numerically with NURBS finite elements in Section 2. Then, two elementary cases study are investigated. A bridge submitted to a local elevation of temperature in Section 3 and a prestressed beam in Section 4

Kinematic enrichment based on the formal asymptotic expansion

The 3D problem

We consider a beam occupying the prismatic domain Ω (Figure 1) with a length L and a cross-sectional typical size h. The boundary ∂Ω is the union of the lateral (free) surface ∂Ω t and the two end sections S ± (clamped). The longitudinal coordinate is x 3 and the section coordinates are x 1 and x 2 denoted as x α 1, the corresponding reference frame is denoted (O, e 1 , e 2 , e 3 ) where O is an arbitrary point of the plane x 3 = 0.

The constitutive material of the beam is only function of the section coordinates x α and invariant in the longitudinal direction. For convenience and without limitation, the corresponding fourth order stiffness tensor C(x α ) is assumed monoclinic with respect to a plane of normal e 3 :

C αβγ3 = C 333α = 0. ( 1
)
The beam is only subjected to eigenstrains ε * i j . The corresponding 3D linear boundary value problem writes as:

               div x σ = 0 on Ω, σ = C : (ε -ε * ) on Ω, ε = ∇ s x u on Ω, σ • n = 0 on ∂Ω t , u = 0 on S ± , ( 2
)
where n is the outer normal to ∂Ω t , ∇ s x is the symmetric part of the 3D gradient operator and div x is the 3D divergence operator. Casting the weak form of this boundary value problem reveals that, in addition to the classical regularity of u, the eigenstrain load needs to be square integrable.

Note that another way to introduce eigenstrains would be to turn it into a body force f * = -div x (C : ε * ) and a force per unit surface T * = (C : ε * ) • n on ∂Ω t . Then it would be possible to use the method presented in [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF]. This remark shows that, without considering fields regularity, the optimality result presented in [START_REF] Miara | A Galerkin spectral approximation in linearized beam theory[END_REF] may be adapted to the present situation. However, this approach involves the preliminary computation of the divergence div x (C : ε * ), which may be a source of numerical imprecision, especially with heterogeneous sections or eigenstrains with low regularity. This motivates the present direct formulation of the higher-order beam theory.

Scaling and variable separation

A new set of coordinates y i is defined from the global coordinates,

(x 1 , x 2 , x 3 ) = (hy 1 , hy 2 , L y 3 ) , (3)
which rewrites the derivation operator as:

∇ x = 1 L ∇ y 3 + 1 η ∇ y α , (4)
as well as the integration over the domain Ω as:

∫ Ω dΩ = ∫ L 0 ∫ S dx α dx 3 = L 3 η 2 ∫ 1 0 ∫ S 0 dy α dy 3 = L 3 η 2 ∫ 1 0 dy 3 , (5)
where η = h L 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 eigenstrain is assumed to have the following separated form:

ε * i j = η dij (y α )T(y 3 ), (6)
where dij (y α ) is the eigenstrain distribution in the section and T(y 3 ) is its longitudinal variation. In this section, capital letters denote functions of only the y 3 coordinate (except for C) and • denotes functions of only in-section coordinates y α .

Expansion

The asymptotic expansion method is a formal procedure in which all fields are assumed sufficiently smooth. It yields a cascade of in-section and longitudinal boundary value problems which are classically solved recursively. In the present case, only the in-section 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] Buannic | Higher-order effective modeling of periodic heterogeneous beams. I. Asymptotic expansion method[END_REF][START_REF] Sanchez-Palencia | Non-Homogeneous Media and Vibration Theory[END_REF][START_REF] Trabucho | Mathematical modelling of rods[END_REF][START_REF] Zhao | Justification of the asymptotic expansion method for homogeneous isotropic beams by comparison with the Saint-Venant solution[END_REF]:

u = L U 0 α (y 3 )e α + η u 1 + η 2 u 2 + ... , (7)
ε = ε 0 + ηε 1 + η 2 ε 2 + ..., (8)
σ = σ 0 + ησ 1 + η 2 σ 2 + ... (9)
and introduced in the 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 p of η is identified. 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. This motivates the scaling of the eigenstrain in Equation [START_REF] Carrera | Beam Structures[END_REF]. For p ∈ N, each compatibility equations, boundary conditions and constitutive equations for p and equilibrium equations for p -1 yield an auxiliary problem on the cross-section which splits in two uncoupled boundary value problems.

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

T p+1 :                          σ p αβ,β + σ p-1 α3,3 = 0 on S 0 , (10a)
σ p αβ = C αβγδ ε p δγ -δ p1 dδγ T + C αβ33 ε p 33 -δ p1 d33 T on S 0 , (10b)
σ p 33 = C 33αβ ε p βα -δ p1 dβα T + C 3333 ε p 33 -δ p1 d33 T on S 0 , (10c)
ε p αβ = u p+1 (α,β) , ε p 33 = u p 3,3 on S 0 , (10d)
σ p αβ n β = 0 on ∂S 0 . ( 10e
)
where σ -1 = 0 and δ 1p = 1 if p = 1 and δ 1p = 0 else. Transposing the results from [START_REF] Amrouche | On the characterizations of matrix fields as linearized strain tensor fields[END_REF][START_REF] Ciarlet | Another approach to linearized elasticity and Korn's inequality[END_REF], for a simply connected cross-section and regular enough C and loadings, this boundary value problem on the displacement u p+1 α is a pure traction problem which is well-posed provided that the applied load is globally self-equilibrating for in-section translations and rotation:

σ p-1
α3,3 = 0 and

y β βα σ p-1 α3,3 = 0, ( 11
)
where αβ is the permutation operator: 11 = 22 = 0, 12 = +1, 21 = -1. Then, under condition [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF], the solution is uniquely defined up to the following rigid motion of the section in its plane:

u R, p+1 α = U p+1 α (y 3 ) + y β βα Θ p+1 (y 3 ). ( 12
)
where U p+1 α is a transverse displacement and Θ p+1 a twist rotation.

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

W p+1 :                  σ p 3α,α + σ p-1 33,3 = 0 on S 0 , ( 13a
)
σ p α3 = C α3β3 2 ε p β3 -δ p1 dβ3 T on S 0 , ( 13b
)
2ε p α3 = u p+1 3,α + u p α,3
on S 0 , (13c)

σ p α3 n α = 0 on ∂S 0 . ( 13d
)
Again, for a simply connected cross-section and regular enough C and loadings, this boundary value problem on the displacement u p+1 3 is well-posed if the applied load is globally self-equilibrating for the longitudinal translation:

σ p-1 33,3 = 0. (14)
In this case, the solution is uniquely defined up to a uniform longitudinal displacement:

u R, p+1 3 = U p+1 3 (y 3 ). ( 15
)
Resultants and macroscopic equilibrium equations The rigid motion of the section suggests the following definition of the beam resultants at each order p ≥ 0:

N p 3 = σ p 33 , M p α = y α σ p 33 , M p 3 = y β βα σ p α3 and V p α = σ p α3 , (16)
where N p 3 is the normal traction, M p α are the bending moments2, M p 3 is the moment of torsion and V p α are the shear forces. These resultants must comply with the following beam equilibrium equations for each p ≥ 0:

               N p 3,3 = 0, (17a)
M p α,3 = V p+1 α , (17b)
M p 3,3 = 0, (17c)
V p α,3 = 0. ( 17d
)
Indeed, from the in-section equilibrium equations (10a) and (13a):

σ p+1 3α,α + σ p 33,3 = ∫ ∂S 0 σ p+1 α3 n α dl + N p 3,3 = 0, y β σ p+1 3α,α + σ p 33,3 = -V p+1 β + ∫ ∂S 0 y β σ p+1 α3 n α dl + M p β,3 = 0, y γ γα σ p+1 αβ,β + σ p α3,3 = -y γ,β γα σ p+1 αβ + ∫ ∂S 0 y γ γα σ p+1 αβ n β dl + M p 3,3 = 0, σ p+1 αβ,β + σ p α3,3 = ∫ ∂S 0 σ p+1 αβ n β dl + V p α,3 = 0. ( 18
)
Note that equilibrium equations (17a-c-d) are identical to conditions [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF][START_REF] Iesan | Saint-Venant's problem for inhomogeneous and anisotropic elastic bodies[END_REF]. Hence, satisfying beam equilibrium equations ensures that T p and W p have a unique solution up to the rigid motions [START_REF] Giunta | A thermal stress finite element analysis of beam structures by hierarchical modelling[END_REF][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].

2 It will appear that M p α is the working conjugate to the curvature U p α,33 and not the conventional bending moment. Indeed, the classical definition is m p α = αβ y β σ p 33 = αβ M p β . This choice is made for convenience.

Cascade resolution

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 u 1 α is a rigid motion and the corresponding stress is null:

u 1 α = U 1 α (y 3 ) + y β βα Θ 1 (y 3 ) and σ 0 αβ = 0, σ 0 33 = 0. ( 19
)
Here, Θ 1 appears as the leading order angle of twist and U 1 α as the next order macroscopic transverse displacement.

Longitudinal displacement

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

W 1 :                σ 0 3α,α = 0 on S 0 , ( 20a
)
σ 0 α3 = C α3β3 2ε 0 β3 on S 0 , ( 20b
)
2ε 0 α3 = u 1 3,α + U 0 α,3 on S 0 , (20c)
σ 0 α3 n α = 0 on ∂S 0 . ( 20d
)
The applied load is self-equilibrating and the solution of this boundary value problem writes as:

u 1 3 = U 1 3 + y α U 0 α,3 and σ 0 α3 = 0, (21)
where U 0 α,3 appears as the bending rotation and U 1 3 as the leading order longitudinal displacement.

Second-order problems

Transverse displacement The transverse displacement u 2 α is derived through:

T 2 :                          σ 1 αβ,β = 0 on S 0 , ( 22a
)
σ 1 αβ = C αβγδ ε 1 δγ -dδγ T + C αβ33 ε 1 33 -d33 T , on S 0 , ( 22b
)
σ 1 33 = C 33αβ ε 1 βα -dβα T + C 3333 ε 1 33 -d33 T on S 0 , ( 22c
)
ε 1 αβ = u 2 (α,β) , ε 1 33 = U 1 3,3 + y α U 0 α,33
on S 0 , (22d)

σ 1 αβ n β = 0 on ∂S 0 . (22e)
Again, the applied load is globally self-equilibrating. The solution of this boundary value problem parametrized by the elongation U 1 3,3 , the curvatures U 0 α,33 and the eigenstrain T writes as the linear superposition of each contribution:

u 2 α = ũe 3 α U 1 3,3 + ũχ 1 α U 0 1,33 + ũχ 2 α U 0 2,33 + ũT α T + U 2 α + y β βα Θ 2 , ( 23
)
where ũe 3 α , ũχ 1 α , ũχ 2 α are in-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. Finally, ũT α is a transverse Poisson's effect related to the eigenstrain. In order to be uniquely defined, the following constraints are applied to all these in-section displacements: ũα = 0 and

y β βα ũα = 0. ( 24
)
Longitudinal displacement The longitudinal displacement u 2 3 complies with:

W 2 :                  σ 1 3α,α = 0 on S 0 , (25a)
σ 1 α3 = C α3β3 2 ε 1 β3 -dβ3 T on S 0 , ( 25b
) 2ε 1 α3 = u 2 3,α + y β βα Θ 1 ,3 + U 1 α,3
on S 0 , (25c)

σ 1 α3 n α = 0 on ∂S 0 . ( 25d
)
Copyright line will be provided by the publisher The applied load is globally self-equilibrating. The solution of this boundary value problem parametrized by Θ 1 ,3 , U 1 α,3 and T writes as the linear superposition of each contribution:

u 2 3 = ũχ 3 3 Θ 1 ,3 + ũT 3 T + U 2 3 + y α U 1 α,3 , (26)
where ũχ 3 3 is the torsion warping illustrated in Figure 3a and ũT 3 a warping related to the eigenstrain. Indeed, the displacement ũχ 3 3 is exactly the solution of the Neumann problem for Saint Venant's torsion. Again, these warpings are constrained as follows:

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

           σ 1 αβ = σe 3 αβ U 1 3,3 + σχ 1 αβ U 0 1,33 + σχ 2 αβ U 0 2,33 + σT αβ T, (28a)
σ 1 α3 = σχ 3 α3 Θ 1 ,3 + σT α3 T, (28b)
σ 1 33 = σe 3 33 U 1 3,3 + σχ 1 33 U 0 1,33 + σχ 2 33 U 0 2,33 + σT 33 T . ( 28c
)
Expressing the traction and bending moments leads to the following constitutive equations:

           N 1 3 = A 3 U 1 3,3 + S 1 U 0 1,33 + S 2 U 0 2,33 + N T 3 T, (29a)
M 1 1 = S * 1 U 1 3,3 + D 1 U 0 1,33 + D 12 U 0 2,33 + M T 1 T, (29b)
M 1 2 = S * 2 U 1 3,3 + D * 12 U 0 1,33 + D 2 U 0 2,33 + M T 2 T, (29c)
where:

         A 3 = σe 3 33 , S 1 = σχ 1 33 , S 2 = σχ 2 33 , N T 3 = σT 33 , S * 1 = y 1 σe 3 33 , D 1 = y 1 σχ 1 33 , D 12 = y 1 σχ 2 33 , M T 1 = y 1 σT 33 , S * 2 = y 2 σe 3 33 , D * 12 = y 2 σχ 1 33 , D 2 = y 2 σχ 2 33 , M T 2 = y 2 σT 33 . (30)
The modulus A 3 is the traction stiffness, D 1 and D 2 are the bending stiffnesses. The moduli N T 3 , M T 1 and M T 2 are the traction and bending moments induced by the eigenstrain. It is proved in Appendix (A) that S * 1 = S 1 , S * 2 = S 2 and D * 12 = D 12 . The stiffnesses S 1 and S 2 are related to the first moments of inertia. Furthermore, there is a position for O, the origin of the reference frame, such that S 1 = S 2 = 0 and rotating the reference frame with respect to e 3 , there is an angle such that D 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:

           N 1 3 = A 3 U 1 3,3 + N T 3 T, (31a)
M 1 1 = D 1 U 0 1,33 + M T 1 T, (31b)
M 1 2 = D 2 U 0 2,33 + M T 2 T . ( 31c
)
Similarly the torsion is expressed as function of the macroscopic displacements and the eigenstrain:

M 1 3 = D 3 Θ 1 ,3 + M T 3 T . ( 32
)
where the torsion stiffness and the torsion induced by the eigenstrain are:

D 3 = y γ γα σχ 3 α3 and M T 3 = y γ γα σT α3 . (33)
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 ( 31) and ( 32) would be simply coupled in such a case.

Third-order problems

Transverse displacement The transverse displacement u 3 α is derived through T 3 and loaded by Θ 1 ,33 , T ,3 , U 2 3,3 and U 1 α,33 . The applied load is self-equilibrating in translation. Indeed:

σ 1 α3 = V 1 α = M 0 α,3 = 0. ( 34
)
Furthermore, from the macroscopic equilibrium in torsion (17c) and the constitutive equation [START_REF] Zhao | Justification of the asymptotic expansion method for homogeneous isotropic beams by comparison with the Saint-Venant solution[END_REF], it is possible to express

Θ 1
,33 as function of T ,3 :

M 1 3,3 = D 3 Θ 1 ,33 + M T 3 T ,3 = 0. ( 35
)
Substituting this relation in T 3 ensures that it is equilibrated in rotation and leads to:

T 3 :                          σ 2 αβ,β + σT α3 - M T 3 D 3 σχ 3 α3 T ,3 = 0 on S 0 , ( 36a
)
σ 2 αβ = C αβγδ ε 2 δγ + C αβ33 ε 2 33 , on S 0 , ( 36b
)
σ 2 33 = C 33αβ ε 2 βα + C 3333 ε 2 33 on S 0 , ( 36c
)
ε 2 αβ = u 3 (α,β) , ε 2 33 = ũT 3 - M T 3 D 3 ũχ 3 3 T ,3 + U 2 3,3 + y α U 1 α,33
on S 0 , (36d)

σ 2 αβ n β = 0 on ∂S 0 . ( 36e
)
The solution of this boundary value problem writes as:

u 3 α = ũT∇ α T ,3 + ũe 3 α U 2 3,3 + ũχ 1 α U 1 1,33 + ũχ 2 α U 1 2,33 + U 3 α + y β βα Θ 3 with ũT∇ α = 0 and y β βα ũT∇ α = 0 (37)
and yields only one new transverse displacement localization related to the first-order variations of T.

Copyright line will be provided by the publisher

Longitudinal displacement

The longitudinal displacement u 3 3 complies with W 3 and is loaded by Θ 2 ,3 , T ,3 , U 1 3,33 and U 0 α,333 . The applied load is not globally self-equilibrating for each individual loading. From the constitutive equation (31a) and the equilibrium equation (17a) with p = 1 it is deduced that :

U 1 3,33 = - N T 3 A 3 T ,3 . (38)
Substituting this in W 3 ensures it is well-posed. Furthermore, it is also convenient to ensure that the load is also selfequilibrating in bending. This is obtained, from the constitutive equations (31b) and (31c) and the equilibrium equation (17b) for p = 1:

U 0 1,333 = V 2 1 -M T 1 T ,3 D 1 and U 0 2,333 = V 2 2 -M T 2 T ,3 D 2 . ( 39
)
Inserting these relations in problem W 3 leads to:

W 3 :                    σ 2 3α,α + σχ 1 33 V 2 1 D 1 + σχ 2 33 V 2 2 D 2 + σT 33 - N T 3 A 3 σe 3 33 - M T 1 D 1 σχ 1 33 - M T 2 D 2 σχ 2 33 T ,3 = 0 on S 0 , ( 40a
)
σ 2 α3 = C α3β3 2ε 2 β3 on S 0 , ( 40b
) 2ε 2 α3 = u 3 3,α + ũχ 1 α V 2 1 D 1 + ũχ 2 α V 2 2 D 2 + ũT α - N T 3 A 3 ũe 3 α - M T 1 D 1 ũχ 1 α - M T 2 D 2 ũχ 2 α T ,3 + y β βα Θ 2 ,3 + U 2 α,3 on S 0 , (40c)
σ 2 α3 n α = 0 on ∂S 0 . ( 40d
)
The solution is parametrized by the shear forces V 2 α the first-order variations of the eigenstrain T ,3 and higher-order displacements. It writes as the linear superposition of each contributions:

u 3 3 = ũV 1 3 V 2 1 + ũV 2 3 V 2 2 + ũT∇ 3 T ,3 + ũχ 3 3 Θ 2 ,3 + y α U 2 α,3 + U 3 3 . with ũ3 = 0 (41)
The longitudinal displacements ũV 1 3 and ũV 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 (40a) 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 beams.

Macroscopic constitutive equations

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

           σ 2 αβ = σe 3 αβ U 2 3,3 + σχ 1 αβ U 1 1,33 + σχ 2 αβ U 1 2,33 + σT∇ αβ T ,3 , (42a)
σ 2 α3 = σV 1 α3 V 2 1 + σV 2 α3 V 2 2 + σT∇ α3 T ,3 + σχ 3 α3 Θ 2 ,3 , (42b)
σ 2 33 = σe 3 33 U 2 3,3 + σχ 1 33 U 1 1,33 + σχ 2 33 U 1 2,33 + σT∇ 33 T ,3 . (42c)
Expressing the second-order traction and bending moments leads to:

           N 2 3 = A 3 U 2 3,3 + N T∇ 3 T ,3 , (43a)
M 2 1 = D 1 U 1 1,33 + M T∇ 1 T ,3 , (43b)
M 2 2 = D 2 U 1 2,33 + M T∇ 2 T ,3 , (43c)
where N T∇

M 2 3 = D 3 Θ 2 ,3 + y S α αβ V 2 β + M T∇ 3 T ,3 , (44)
where the torsion induced by the variations of eigenstrains is M T∇ 3 = y γ γα σT∇ α3 and the shear center of the beam is defined as:

y S 1 = -y α αβ σV 2 β3 and y S 2 = y α αβ σV 1 β3 . ( 45
)
When the section presents two axis of symmetry, the shear center is in O but this is not always true.

Copyright line will be provided by the publisher

Fourth-order and higher-order problems

The induction process may be pursued any higher order. Indeed, noticing that V p α,3 = 0 and following the same procedure as for T 3 and W 3 it appears that T 4 and W 4 are formally identical to T 3 and W 3 , incrementing the orders and depending on the second derivative of T. Hence, going higher order leads to the derivation of displacement localizations related to higher derivatives of T relevant for faster variations of T.

The use of the asymptotic expansion method is based on the scaling in equation ( 3). Hence the rescaled coordinates y i have been used in the expression of the auxiliary problems T p and W p . However, the distinction between the two sets of coordinates is no longer necessary in practice once the section modes are computed. The use of the coordinates y i is therefore dropped in all what follows and replaced by the use of the coordinates x i .

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 U p i , U p α,3 and Θ p 3 . They are respectively related to the following displacement modes:

ũU 1 = 1 0 0 , ũU 2 = 0 1 0 , ũU 3 = 0 0 1 , ũΘ 2 = 0 0 -x 1 , ũΘ 1 = 0 0 x 2 , ũΘ 3 = -x 2 x 1 0 . ( 46
)
Second, the six correctors related to the six beam resultants3 were derived: ũe 3 α , ũχ 1 α , ũχ 2 α , ũχ 3 3 , ũV 1 3 , ũV 2 3 . They are also referred to as Saint Venant's modes [START_REF] Fatmi | A refined 1D beam theory built on 3D Saint-Venant's solution to compute homogeneous and composite beams[END_REF][START_REF] Iesan | Saint-Venant's problem for inhomogeneous and anisotropic elastic bodies[END_REF]. Note that ũχ 3 3 is the warping used by Benscoter [START_REF] Benscoter | A theory of torsion bending for multicell beams[END_REF]. Third, exactly as in [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF], a family of modes related to the eigenstrain loading and its longitudinal variations was obtained: ũT , ũT∇ , ũT∇ 2 ... Finally, this suggests gathering all these modes in the following approximation for the 3D displacement:

u = n i=1 ũi (x α )X i (x 3 ) (47)
where n is the number of modes and X k (x 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 are linearly independent. But in the eventuality where the components of the eigenstrains would be described respectively with different longitudinal functions (T 1 , T 2 ...), redundancies may occur. In this case it is necessary to orthogonalize the basis of modes.

Numerical approximation of the higher-order beam model

The discretization of the section used for the resolution of the auxiliary problems is presented. The discretization of the longitudinal beam element and the numerical approximation of the total displacement are then exposed. This approximated expression of the displacement is used to formulate the expression of the minimum of potential energy, leading to the 1D boundary value problem.

Numerical resolution of the auxiliary problems

For each order p, the weak forms of T p and W p are expressed. The resolution of the formulations obtained yields the displacement modes. The numerical resolution is operated by a discretization of the section with finite elements. The elements chosen here are triangles, the interpolation being quadratic. More precisely, the isoparametric expressions of the interpolation functions are:

N 1 (a 1 , a 2 ) = a 1 (2a 1 -1) , N 2 (a 1 , a 2 ) = a 2 (2a 2 -1) , N 3 (a 1 , a 2 ) = (1 -a 1 -a 2 )(1 -2a 1 -2a 2 ) N 4 (a 1 , a 2 ) = 4a 1 a 2 , N 5 (a 1 , a 2 ) = 4a 2 (1 -a 1 -a 2 ) , N 6 (a 1 , a 2 ) = 4a 1 (1 -a 1 -a 2 ). ( 48
)
But there is no limitation in the choice of the type of finite elements. Stress and strain are computed at the three Gauss points of each element while displacement is computed at the nodes of the sectional mesh. Note that T p and W p are also loaded with eigenstrains which must be expressed at the Gauss points (in Equation (40c) for instance). Hence, for the resolution of T p and W p , the displacement loading the auxiliary problems must be interpolated at the Gauss points thanks to the interpolation functions.

3 From the traction, bending and torsion constitutive equations ( 31) and ( 32), ũe 3 α , ũ χ β α and ũ χ 3 3 may directly be expressed as function of N 3 , M β and M 3 .

Numerical approximation of the beam element

To solve the higher-order beam model, the one-dimensional longitudinal problem needs to be formulated. The kinematic variables X i are expressed with interpolation functions:

X i (x 3 ) = m k=1 N k (x 3 ) Ũi,k , (49)
with m the number of interpolation functions. Ũi,k is the generalized displacement associated to the mode i and the interpolation function k.

We use NURBS (Non-Uniform Rational B-Splines). NURBS are a generalization of B-Splines and Bézier curves). NURBS basis functions are defined by a degree p and a non-uniform knot vector Ξ = {ξ 1 , ξ 2 , ..., ξ l } with l the number of knots. This domain partition allows the definition of basis functions, defined recursively by the Cox-de Boor's formulas. 0 th order functions are piecewise constants:

N k,0 (ξ) = 1 if ξ k ≤ ξ < ξ k+1 , 0 otherwise. ( 50
)
Higher-order basis functions (p > 0) are then defined by

N k+1, p (ξ) = ξ -ξ k ξ k+p -ξ k N k, p-1 (ξ) + ξ k+p+1 -ξ ξ k+p+1 -ξ k+1 N k+1, p-1 (ξ). (51)
It defines m = l + p -1 NURBS basis functions for the order p, each function being C p-1 . It can be noticed that NURBS basis functions form a partition of unity:

m k=1 N k, p (ξ) = 1. ( 52
)
Moreover the basis functions are in general not interpolatory except at the ends of the beam. This motivates distinction between knots and "nodes". Each basis function is a polynomial with a compact support contained in the interval ξ k , ξ k+p+1 . Their use is therefore convenient to describe very fast variations of a field: refining the mesh defines more basis functions with short supports and affords a better localization of the field. This is an important advantage of NURBS over Lagrange polynomials: the support of a Lagrange polynomial is the whole domain [ξ 0 , ξ l ], and a Lagrange polynomial can show important oscillations on this interval in certain cases, a situation often called the Runge's phenomenon. The number of interpolation functions m depends on the refinement of the mesh and the chosen polynomial degree. In practice, in order to better describe fast variations of a mechanical field, the degree of the NURBS shouldn't be too high, since the higher the degree is, the larger the supports of the functions are. An illustration of NURBS basis functions with a variable refinement of the mesh is presented in Figure 4.

Copyright line will be provided by the publisher

Formulation of the beam problem

For a fixed degree p of NURBS interpolation, the total displacement can now be approximated as:

û = n i=1 ũi (x) m k=1 N k, p (x 3 ) Ũi,k . ( 53
)
The stationarity of the minimum of the potential energy based on equations (2) requires:

Find u ∈ K, ∫ Ω ε(u) : C : ε( û)dΩ = ∫ Ω ε * : C : ε( û)dΩ, ∀ û ∈ K, (54)
where K = {regular enough u | u = 0 on S ± } is the set of kinematically compatible fields. The injection of the approximation û (53) into equation ( 54) leads to the classical linear system:

K Ũ = F th , (55)
where Ũ = Ũi,k with 1 ≤ i ≤ n and 1 ≤ k ≤ m. The resolution of (55) yields the value of all the kinematic unknowns. Equation (53) then gives the total displacement.

2.6.4

Locking study NURBS offers many advantages but they are not yet free from locking [START_REF] Echter | Numerical efficiency, locking and unlocking of NURBS finite elements[END_REF]. A cantilever beam has been studied in order to investigate this phenomenon and define the range of applicability of the present discretization. The beam has a constant square section of S = 1 × 1 m 2 , is clamped for x 3 = 0 and is loaded for x 3 = L by a vertical force F(L) = F 0 /L 3 . The young's modulus and the Poisson's ratio are E = 35 GPa and ν = 0. The kinematics of the beam comprises the rigid and Saint Venant's modes, and also up to three modes associated to the load applied and its gradients as described in [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF] and so called force modes. This means that there are only 15 kinematic DOF per section. The beam element is interpolated with NURBS of various degrees, defined by 11 knots evenly distributed in [0, L], with ξ 1 = 0 and ξ 11 = L. Hence the total number of DOF is the same for all calculations.

The analytical expression of the vertical displacement at x 3 = L for Euler's model is well-known. Hence, the following relative error between Euler's model and the higher-order beam model is defined:

e Euler = u beam 2 -u Euler 2 u Euler 2 , with u Euler 2 = F L 3 3E I , and
u beam 2 = 1 S u 2 (x 1 , x 2 , L) , ( 56
)
where S is the area of the section. The Euler solution delivers a good indication on the numerical behavior of the present finite element when the slenderness increases since it is known to be the limit model for large slenderness. First, the influence of the order of the NURBS on the locking phenomenon is investigated for an enrichment limited to 2 force modes (the applied load and its first gradient). Figure 5a shows the relative error as function of the slenderness of the beam for several NURBS orders. For low slenderness, all relative errors are high. Indeed, for such slenderness, it is the Euler model which is not valid. When increasing the slenderness, all beam models are expected to converge towards the Euler solution. This is not the case when locking occurs. The first order NURBS basis functions are identical to the Lagrange polynomials classically used. Therefore the same locking phenomenon is observed and e Euler goes to 1 for high slenderness. The second order NURBS also suffers from locking after L/h = 20. Third order NURBS tends to the Euler's solution with a best match for L/h = 50 and e Euler < 1.10 -4 . Locking appears from L/h = 100. However, distance to Euler's solution remains only about 1% for L/h = 500.

Second, the influence of the number of force modes is investigated with second order NURBS. Figure 5b shows the relative error as function of the slenderness of the beam for several number of force modes. It shows that increasing the number of force modes does not influence the slenderness above which locking occurs.

Therefore, for common situations with L/h < 100, second-order or third order NURBS can be considered as adapted for the interpolation of the beam element. In case of higher slenderness, increasing the number of longitudinal knots will mitigate locking.

Application to a thermal load

Thermal load on a clamped beam

To illustrate the method, we consider a beam clamped at both ends, and loaded with a thermal load. The section of the beam is represented in Figure 6. The beam is 40 m long. The volume of the beam is denoted by Ω. A thermal load is applied on Copyright line will be provided by the publisher 6 and7. The thermal load generates an isotropic strain Ω T :

ε * i j = α T -T ref δ i j ( 57
)
where α is the dilatation coefficient, here fixed at 12.10 -6 , and T ref is the reference temperature The material considered is homogeneous and isotropic. We consider the following values of the Young's modulus and the Poisson's ratio: E = 210 GPa and ν = 0.3.

T -T ref = T max -T ref = 40 • C on Ω T 0 • C on Ω \ Ω T ( 58

Reference solution

A 3D model computed with Code_Aster is taken as a reference solution for our study. A convergence study related to the refinement of the mesh is first carried out. 4 meshes are investigated. Each mesh is made of hexahedric elements, each element being quadratic. The refinements of the 4 meshes are described in Table 1.

In order to assess the convergence of the different computations, the values of the strains are observed in the middle of the eigenstrain area (x 3 = 15 m) and close to the strain discontinuity (x 3 = 9, 90 m). The results are extracted at the center of the slab's thickness (x 2 = 0.65 m) and presented in Figure 8.

The results quickly converge. The Figure 8f representing the strain ε 33 at x 3 = 9, 90 m shows a slower convergence. Figure 8f shows that convergence is reached for mesh 3, since the curves of mesh 3 and mesh 4 are almost overlapping. Results obtained with mesh 2 shows a gap with the converged values about 5%. Mesh 4 is chosen as the reference solution. In order to compare the results of the solutions presented in Figures 8a-8f we define the following L 2 -estimator for a given x 3 and for x 2 = 0.65 m:

e L 2 (ε i j ) = ∫ L 1 ε i j (x 1 ) -ε (4) i j (x 1 ) 2 dx 1 1/2 ∫ L 1 ε (4) i j (x 1 ) 2 dx 1 1/2 (59)
where [START_REF] Buannic | Higher-order effective modeling of periodic heterogeneous beams. II. Derivation of the proper boundary conditions for the interior asymptotic solution[END_REF] is the width of the considered section. The relative errors between meshes 1, 2 and 3 and mesh 4 are gathered in Table 2. The maximum relative errors of each solution give an estimate of the accuracy of the reference solution and will be compared to the error of the higher-order beam approximation.

L 1 = [-5,

Higher-order beam solution

For the computation of the present method, the section is meshed with 712 triangle elements, as shown in Figure 9. Each element is quadratic. The longitudinal and transverse components of the eigenstrain are assumed to vary accordingly to two different longitudinal functions T 1 and T 2 :

ε * α3 = 0 ε * αβ = d(x γ )T 1 (x 3 )δ αβ ε * 33 = d(x γ )T 2 (x 3 ) (60)
where d(x α ) is defined as:

d(x γ ) = α(T max -T ref ) if x γ ∈ S T 0 if x γ S T (61)
Copyright line will be provided by the publisher The use of two different functions T 1 and T 2 is an arbitrary choice based on experience. Since the components of the eigenstrain follow two different longitudinal functions, the 2D computation twice as many eigenstrain modes, which represents a kinematics richer than if a single function had been used. The first orthonormalized eigenstrain modes computed up to the second gradients of T 1 and T 2 are given in Figure 9.

The longitudinal mesh is composed of 42 knots: 41 knots are evenly distributed between 0 and L and an additional knot is placed at x 3 = 9, 90. This last knot is added in order to better compute displacements and strains close to the eigenstrain discontinuity. This longitudinal mesh is less refined than the one of the reference model which has 400 elements along its longitudinal axis.

The model is first computed without eigenstrain modes. Thus the only modes considered are rigid and Saint Venant's modes as presented in Section 2.5. We call this solution S 0 with 12 kinematic DOF per section. The model is then computed Table 3 Comparison of the solutions up to the second gradient of the eigenstrain u ∇T 2 i , i = 1, 2. This should yield a maximum of 6 modes (u 2 ). In our case of study u T 1 = u T 2 . We therefore get only 5 additional eigenstrain modes which are presented in Figure 9. It can be noticed than u T 1 = u T 2 does not imply u T 1 ∇ = u T 2 ∇ , since the recursive definition of transversal problems T p+1 and warping problems W p+1 is second order: they imply loading terms as function of displacements of the order p -1 and p.

T 1 , u T 2 , u ∇T 1 , u ∇T 2 , u ∇T 2 1 and u ∇T 2
We call this second solution S with 17 kinematic DOF per section. The comparison between S and S 0 highlights the contribution of the eigenstrain modes to the global response. The main features of the 3D solution and the solution S are presented in Table 3.

Comparison of the solutions

The deformed structure obtained with the 3D computation and with the higher-order beam model are presented in Figure 10 (the same scale is used for both figures). The color map shows the norm of the 3D displacement. Figure 10 illustrates the ability of the higher-order beam model to capture the cross-sectional displacements. The higher-order beam model can also satisfactorily capture the discontinuities of loads: the heated area is easily identifiable in Figure 10b.

A finer comparison of the present beam model to the reference model is now performed. Like with the reference solution, the results presented hereafter are all extracted at the mid surface of the slab (x 2 = 0.65 m), as defined in Figure 6. In order to compare the different solutions we define as previously the following L 2 -estimators:

e L 2 (u i ) = ∫ L 1 u i (x 1 ) -u 3D i (x 1 ) 2 dx 1 1/2 ∫ L 1 u 3D i (x 1 ) 2 dx 1 1/2 , e L 2 (ε i j ) = ∫ L 1 ε i j (x 1 ) -ε 3D i j (x 1 ) 2 dx 1 1/2 ∫ L 1 ε 3D i j (x 1 ) 2 dx 1 1/2 (62)

Results in the heated area

The first results presented are computed in the middle of the heated area at x 3 = 15 m. Displacement and strain are shown in Figures 11. Shear strains being less significant are not presented here.

The displacement u 1 matches well with the 3D solution, since e L 2 (u 1 (x 3 = 15 m)) = 1, 69.10 -2 . The discontinuity of the eigenstrain clearly appears in Figure 11a. Displacement obtained with S 0 shows that the beam modes are not sufficient to describe the global response to this specific load. The vertical displacement u 2 matches with the 3D model with an error e L 2 (u 2 (x 3 = 15 m)) = 2, 33.10 -2 , but a higher expansion order would bring even more satisfying results. Looking at the longitudinal displacement, the 3D solution is more rigid than the beam solution, and e L 2 (u 3 (x 3 = 15 m)) = 4, 00.10 -2 . The accuracy of the solution could be improved by refining the longitudinal mesh. According to Figures 11d to 11f, the axial strains computation shows very satisfying prediction. The strain ε 33 shows a little difference with the reference model e L 2 (ε 33 (x 3 = 15 m)) = 2, 45.10 -2 . This error can still be lowered by increasing the number of modes or by refining the longitudinal mesh. As presented in Table 3, the computation of this method is fast, and it provides satisfying results.

Discontinuity and boundary conditions

It is now interesting to compare the response provided by both models close to the strain discontinuity and close to one clamped end. The displacements computed at x 3 = 9, 90 m and x 3 = 1 m are shown in Figure 12. Close to the discontinuity, the prediction of solution S is still satisfying. The cross-sectional displacements matches with the 3D model: e L 2 (u 1 (x 3 = 9, 90 m)) = 1, 42.10 -2 and e L 2 (u 2 (x 3 = 9, 90 m)) = 1, 50.10 -2 . Even if the longitudinal displacement is less satisfying with a L 2 -error e L 2 (u 3 (x 3 = 9, 90 m)) = 6, 94.10 -2 , it appears that the model has no real difficulty in capturing the discontinuities. If needed, the results of u 3 could be improved by increasing the number of longitudinal knots.

The effect of the boundaries on the displacements appears in Figures 12d to 12f ).

Yet the estimation is suitable for engineering practice. It is consistent to find less accurate results close to the boundaries since the boundary conditions are not considered in the computation of the modes. The boundaries requires particular strains that could be captured by adding modes specifically computed for the boundary conditions under consideration. 4 Application to a prestressed beam 4.1 Prestressed cable in a cantilever beam We further illustrate the present model with a cantilever beam prestressed with a steel cable. The section of the beam is represented in Figure 13. The beam is 20 m long. The concrete domain is denoted by Ω c and the steel domain by Ω s . A constant eigenstrain ε * = ε * 33 e 3 ⊗ e 3 is applied in Ω s , with ε * 33 = 7.10 -3 , corresponding to 23 MN tension in the prestressed cable. Both materials are homogeneous and isotropic with:

• (E, ν) concrete = (35 GPa, 0.2)

• (E, ν) steel = (200 GPa, 0.3)
The beam is 20 m long. Note that, a real concrete beam would require additional reinforcement bars as well as a non-linear constitutive behavior. This simplified example is chosen here to illustrate the ability of eigenstrain modes to capture rather fast variations of the strain in the section.

Reference and beam solution

As for the previous example, a convergence study is carried out in order to choose the reference solution. The study is computed with Code_Aster with full 3D computations. The mesh of the chosen reference solution is made of 99680 pentahedric elements, each element being quadratic. This mesh has been constructed by meshing the end section with 2492 triangles, and then by extruding this 2D mesh on 40 section evenly distributed from x 3 = 0 to x 3 = L. For the computation with the higher-order beam model, the section is meshed with 1788 triangle elements. Each element is quadratic. The displacement modes are computed accordingly to this load, using the same decomposition of the eigenstrain as in the previous example. The sectional modes are computed up to the fourth gradient of the eigenstrain. The 5 first orthonormalized modes associated to the eigenstrain u T , u ∇T , u ∇T 2 , u ∇T 3 and u ∇T 4 are represented in Figure 14. These three cross-sectional modes and two warping modes clearly illustrate the action of the cable on the beam. The shift between purely transversal or purely warping modes at each order comes from the monoclinic symmetry of the constitutive material (1) and the absence of shear in the thermal loading of the present case study.

The longitudinal mesh is composed of 41 knots evenly distributed on the length of the beam. The longitudinal mesh used here is the same as the one of the reference model. A first solution only with the rigid and the Saint Venant's modes is computed and called S 0 with 12 kinematic DOF per section. The solution using also the 5 eigenstrain modes is called S and involves 17 kinematics DOF per section. The main features of the 3D model and the solution S are presented in Table 4. Fig. 14 The 5 eigenstrain modes related to the prestress in the steel cable used in case study 2.

Comparison of the solutions

The deformed structure obtained with the 3D computation and with the higher-order beam model are presented in Figure 15 (the same scale is used for both figures). As expected, the prestressed cable compresses and raises the beam. The higherorder beam model captures cross-sectional displacements such as the lowering of the edges of the table, and warping displacements such as the punching effect of the cable which can be observed at the end of the beam. A finer comparison of the models is presented now. Displacements and strains are computed at mid-span of the beam at x 3 = 10 m and close to the clamped end at x 3 = 0.5 m. The results presented hereafter are all extracted at the axis of symmetry of the section for x 1 = 0 m. The shear strains are not presented here because of their very small amplitude. Thanks to the enrichment of its kinematics, the beam model is able to capture the strains locally generated by the cable. However, the strains ε 11 and ε 22 of solution S does not perfectly match with the 3D solution in the prestressed area, leading to the errors e L 2 (ε 11 (x 3 = 10 m)) = 7, 67.10 -2 and e L 2 (ε 22 (x 3 = 10 m)) = 2, 69.10 -1 . The accuracy of these results could be improved by providing a more refined mesh in the prestressed area of the model.

Results close to support

The previous example showed that the results were less satisfactory near the boundary conditions. This is confirmed in the present case by the computations at x 3 = 0.5 m. The beam solution S still correctly matches with the 3D solution, but larger errors are observed: e L 2 (u 2 (x 3 = 0, 5 m)) = 2, 43.10 -1 and e L 2 (u 3 (x 3 = 0, 5 m)) = 5, 12.10 -2 . It can be noticed that even if the amplitude of the vertical displacement close to support is globally 1000 times smaller than in the middle of the beam at x 3 = 10 m, the detailed influence of the cable is still well captured by the beam model.

Regarding strains, the beam model does not perfectly match with the 3D model but remains suitable for engineering practice. The L 2 -error of the axial strains are e L 2 (ε 11 (x 3 = 0, 5 m)) = 1, 35.10 -1 , e L 2 (ε 22 (x 3 = 0, 5 m)) = 2, 48.10 -1 and e L 2 (ε 33 (x 3 = 0, 5 m)) = 1, 56.10 -1 . Solution S does not describe the variation of ε 33 around x 2 = 1, 80 m. This is explained by the influence of the boundary conditions on the displacement. This should be solved by the computation and the addition of new modes specific to the boundary conditions.

Conclusion

A numerical method based on the asymptotic expansion method was recently suggested by Ferradi et al. [START_REF] Ferradi | A model reduction technique for beam analysis with the asymptotic expansion method[END_REF]. The strength of this higher-order model is that the kinematics of the beam is enriched not only with 2D-modes related to the geometry of the section, but also according to the loads applied on the structure. This way, local effects produced by the application of the loads are directly captured by the beam element. Moreover, the extension of the kinematics does not require any a priori knowledge on the solution of the problem. Based on the use of the asymptotic expansion method, this model enables the user to enrich the kinematics of the beam until any expansion order n, and to refine the results subsequently.

Because of the presence of many inelastic phenomena in civil engineering applications, the current paper extends this method to the case of eigenstrains. The introduction of eigenstrains in the equilibrium equations leads to the computation of additional modes specific to the strains applied to the structure. These modes are added to the higher-order kinematics of the beam.

When compared to a 3D computation, the model presented here shows very satisfying results with a significantly reduced computational cost. Indeed, only very few additional modes and corresponding beam DOF (3 to 5) were required for capturing fairly well the applied eigenstrains. These time performances can still largely be improved by more advanced numerical techniques. As shown in this article, the model is able to capture strain discontinuities: the thermal discontinuity introduced in the first example is faithfully described in the results computed by the beam model. This example also illustrates the ability of the model to capture both transverse and longitudinal discontinuities. The second example highlights the ability of the model to render local behavior such as the punching of the prestress cable on the end section of the clamped beam. However, the very close vicinity of the boundaries sometimes seems more difficult to compute, as exposed in the first -5.00E-07 0.00E+00 example. In order to better describe the mechanical behavior next to boundary conditions, the introduction of new modes specific to these boundary conditions is already under investigation.

A Proof of the symmetry of the beam constitutive equation

In this section, the symmetry of the constitutive equation ( 29) is briefly sketched for S η = S * η .

Recall that:

Fig. 2 3 (b) Warping mode ũV 1 3 (c) Warping mode ũV 2 3 Fig. 3

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

  and bending moments induced by the longitudinal variations of the eigenstrain T. Similarly the second-order torsion is expressed as function of the macroscopic displacements and the eigenstrain:

Fig. 4 2

 4 Fig.[START_REF] Buannic | Higher-order effective modeling of periodic heterogeneous beams. I. Asymptotic expansion method[END_REF] 2 nd -order NURBS basis functions for the knot vector {0, 0.05, 0.1, 0.15, 0.2, 0.4, 0.7, 1}

Fig. 5 1 = 2 Fig. 6 Fig. 7

 51267 Fig.[START_REF] Buannic | Higher-order effective modeling of periodic heterogeneous beams. II. Derivation of the proper boundary conditions for the interior asymptotic solution[END_REF] Study of locking on a cantilever beam loaded at its end, NURBS interpolation

4 (e) Strain ε 22 , x 3 = 9 Fig. 8

 422398 Fig.[START_REF] De Saint | De la torsion des prismes: avec des considérations sur leur flexion ainsi que sur l'équilibre des solides élastiques en flexion[END_REF] Axial strains in the heated area and close to the temperature discontinuity for the 3D solution

( a ) 2 Fig. 9

 a29 Fig.9The 5 eigenstrain modes related to the heated area in the cross-section for case study 1.

  Fig. 10 Deformed structure under thermal load (amplification factor = 200)

Strain ε 33 Fig. 11

 3311 Fig. 11 Comparison of displacements and strains of the slab in the heated zone, x 3 = 15 m

Longitudinal displacement u 3 Fig. 12

 312 Fig. 12 Displacement of the heated slab close to the temperature discontinuity x 3 = 9, 90 m (a), (b), (c) and close to supports x 3 = 1 m (d), (e), (f)

( a )

 a Transverse mode u T (b) warping mode u ∇T (c) Transverse mode u ∇T 2 (d) warping mode u ∇T 3 (e) Transverse mode u ∇T 4

Fig. 15 4 . 3 . 1

 15431 Fig. 15 Deformed structure under prestress steel cable (amplification factor = 100)

  .

Strain ε 33 Fig. 16

 3316 Fig.[START_REF] Kim | An asymptotic analysis of composite beams with kinematically corrected end effects[END_REF] Displacement and strain at mid-length of the beam, x 3 = 10 m

Fig. 17

 17 Fig. 17Displacement and strain close to support, x 3 = 0, 5 m

  Fig. 17Displacement and strain close to support, x 3 = 0, 5 m

  problem 22 being linear, each localization field ũA , σ A satisfies the equations for a unit value of the macroscopic field A and vanishing the other macroscopic fields. From the longitudinal constitutive equation (22c):S η = C 33αβ ũχ η α,β + C 3333 x η and S * η = x η C 33αβ ũe 3 α,β + C 3333 ,(64)and the transverse constitutive equation (22b):S η = σe 3 αβ -C αβγδ ũe 3 δ,γ ũχ η α,β + C 3333 x η and S * η = σχ η αβ -C αβγδ ũχ η δ,γ ũe 3 α,β + x η C 3333(65)Because the in-plane stress must satisfy equilibrium equation (22a) as well as free boundary conditions (22e) we have: this and equation (65) it appears that S η = S * η . Similar arguments ensure D * 12 = D 12 .

Table 1

 1 Refinement of the meshes

	x 3 = 15 m	x 3 = 9, 90 m
	ε 11 mesh 1 2, 81 4, 38 7, 21 1, 34 5, 15 7, 53 ε 22 ε 33 ε 11 ε 22 ε 33 mesh 2 0, 40 0, 49 0, 02 0, 46 1, 85 5, 32 mesh 3 0, 98 0, 55 0, 01 0, 41 0, 63 0, 56

Table 2

 2 Strain relative errors e L 2 in percentage of meshes 1, 2 and 3 compared to mesh 4

Table 4

 4 Main features of the 3D solution and solution S

		5m
	0.2m	Ω s	0.14m
	1.80m	Ω c
			y 2
		1m y 3	y 1
	Fig. 13 Section of the prestressed beam
	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

In the following, Greek indices α, β, γ = 1,

denote in-section dimensions and Latin indices i, j, k, l = 1, 2,

3, all three dimensions. Einstein summation convention on repeated indices is used.Copyright line will be provided by the publisher

Copyright line will be provided by the publisher