Fabien Treyssède
email: fabien.treyssede@ifsttar.fr

Patrice Cartraud

A two-dimensional formulation for the homogenization of helical beam-like structures under bending loads

Keywords: homogenization, beam, helical, finite element, spring, strand, cable, bending

In this paper, a two-dimensional formulation is proposed for modeling the mechanical behavior of helical beamlike structures subjected to bending loads. Helical structures include multi-wire cables, which are widely used in engineering applications. An accurate representation of their mechanical behavior typically requires large-size threedimensional finite element models. The proposed formulation, written on the cross-section only, offers a tremendous reduction of computational cost. Based on the asymptotic expansion method, the formulation is derived within the framework of homogenization theory. The first-order approximation of the three-dimensional problem is obtained from the solution of two successive problems: a microscopic three-dimensional problem and a macroscopic onedimensional problem. The latter corresponds to the equilibrium equations of a straight Navier-Bernoulli-Saint Venant beam, which effective elastic properties can be post-processed from the solution of the microscopic problem, rewritten in a helical curvilinear coordinate system (twisting system). Thanks to this coordinate system, we demonstrate that the microscopic problem can be reduced to the cross-section and solved by a two-dimensional finite element analysis. In the twisting system, it is shown that bending loads yet depend on the axial coordinate and require a specific treatment leading to separate variable solutions in the axial variable of complex type. Therefore, this paper advances one step further than previous papers in which the microscopic problem was reduced to a two-dimensional formulation with the drawback that only axial loads (extensional or torsional) could be considered. Numerical results are presented for cylinders, springs, and seven-wire strands. Good agreement with analytical solutions is achieved. Interestingly, the formulation allows an accurate analysis of mechanical contact effects on the homogenized properties.

Introduction

Helical structures are widely used in engineering applications. They can also be encountered, at the microscopic level, in biological or synthetic materials [START_REF] Grason | Colloquium: Geometry and optimal packing of twisted columns and filaments[END_REF]. Typical examples are given by springs, in mechanical systems, and multi-wire cables (strands and ropes), ubiquitous in civil engineering. Cables are usually constituted by helically wound wires and have the ability to support huge axial loads with relatively small bending stiffness. Strands and ropes are essential in many engineering applications such as bridges, power transmission lines, prestressing of concrete, mooring lines, lift machines, etc. Such structures are subjected to environmental degradations (fatigue, corrosion, temperature, wind). Predicting their mechanical behavior is essential to maintain their safety.

In this paper, we restrict our attention to continuously screw-symmetric structures. This includes springs or strands composed of a stack of helical wires wrapped with the same twisting rate around a straight axis (e.g. seven-wire strands). Screw symmetry of discrete type, occurring for instance in multi-layered strands with layers twisted in opposite directions or in stranded ropes (i.e. involving double-helix wires), is left for further works.

The static behavior of strands can be described by simplified models relying on beam theory. This kind of problem can be solved analytically, considering axial loads (extensional or torsional) [START_REF] Costello | Theory of Wire Rope[END_REF][START_REF] Jolicoeur | A numerical comparison of current mathematical models of twisted wire cables under axisymmetric loads[END_REF][START_REF] Ghoreishi | Validity and limitations of linear analytical models for steel wire strands under axial loading, using a 3d fe model[END_REF] as well as bending loads [START_REF] Cardou | Mechanical models of helical strands[END_REF][START_REF] Papailiou | On the bending stiffness of transmission line conductors[END_REF], or numerically with beam finite elements (FE) [START_REF] Nawrocki | A finite element model for simple straight wire rope strands[END_REF][START_REF] Lalonde | Modeling multilayered wire strands, a strategy based on 3d finite element beam-to-beam contacts -part i: Model formulation and validation[END_REF]. The limitation inherent to beam assumptions, such as wire flattening and stress concentration due to contact [START_REF] Ménard | Solid and 3D beam finite element models for the nonlinear elastic analysis of helical strands within a computational homogenization framework[END_REF], yet restricts the application of these simplified models.

To remove such a limitation, three-dimensional FE models have been proposed. In order to obtain a good representation of the mechanical behavior, and in particular, to minimize the edge effects at boundaries, the axial length of the FE model must be large enough. In practice, the length must be of the order of the pitch of the helix. This often leads to large-size and computationally expensive models [START_REF] Ghoreishi | Validity and limitations of linear analytical models for steel wire strands under axial loading, using a 3d fe model[END_REF][START_REF] Cartraud | Computational homogenization of periodic beam-like structures[END_REF][START_REF] Erdonmez | A finite element model for independent wire rope core with double helical geometry subjected to axial loads[END_REF][START_REF] Stanova | Computer modelling of wire strands and ropes part ii: Finite element-based applications[END_REF][START_REF] Zhang | Finite element solutions to the bending stiffness of a single-layered helically wound cable with internal friction[END_REF][START_REF] Foti | Analytical and finite element modelling of the elastic-plastic behaviour of metallic strands under axial-torsional loads[END_REF], because the helix pitch is usually large compared to the transverse dimension of the structure.

A tremendous reduction can be achieved by exploiting the continuous screw symmetry property of the geometry of the helical structure. It has been shown that the computational domain can be restricted to a basic periodic cell, corresponding to a thin three-dimensional FE slice of arbitrary thickness, thanks to specific displacement constraints between the two boundaries of the slice [START_REF] Jiang | A concise finite element model for simple straight wire rope strand[END_REF][START_REF] Jiang | Statically indeterminate contacts in axially loaded wire strand[END_REF]. This reduction can also be achieved from a different derivation following a beam-like homogenization method [START_REF] Ménard | Solid and 3D beam finite element models for the nonlinear elastic analysis of helical strands within a computational homogenization framework[END_REF][START_REF] Messager | Homogenization of helical beam-like structures: application to single-walled carbon nanotubes[END_REF]. Owing to its continuous nature, the screw symmetry can be even considered more efficiently by reducing the computational domain to a two-dimensional model written on the cross-section only. This further reduction has been proposed in Refs. [START_REF] Frikha | Mechanical modeling of helical structures accounting for translational invariance. part 1: Static behavior[END_REF][START_REF] Karathanasopoulos | Two dimensional modeling of helical structures, an application to simple strands[END_REF], later extended to non-linear problems in Refs. [START_REF] Karathanasopoulos | Finite element modeling of the elastoplastic axial-torsional response of helical constructions to traction loads[END_REF][START_REF] Filotto | Nonlinear planar model for helical structures[END_REF]. Nevertheless, these reduced models are currently limited to external loads of extensional and rotational types, excluding bending loads. For a seven-wire strand, note that a length reduction to the sixth of the pitch can be applied thanks to the discrete translational symmetry of the structure. This allows pure bending loads to be considered in a straightforward manner [START_REF] Jiang | A concise finite element model for pure bending analysis of simple wire strand[END_REF], but the computational cost remains relatively significant since a three-dimensional solution is required.

This paper proposes a two-dimensional formulation reduced to the cross-section for modeling the mechanical behavior of continuously screw-symmetric beam-like structures including bending loads. Following Refs. [START_REF] Messager | Homogenization of helical beam-like structures: application to single-walled carbon nanotubes[END_REF][START_REF] Frikha | Mechanical modeling of helical structures accounting for translational invariance. part 1: Static behavior[END_REF], the helical symmetry is accounted for within the framework of homogenization theory. The homogenization approach is based on the asymptotic expansion method and gives the first-order approximation of the three-dimensional problem from the solution of two successive problems: a microscopic three-dimensional problem and a macroscopic onedimensional problem. As shown in this paper, the former can be reduced to the cross-section and solved by a twodimensional FE analysis. The latter corresponds to the equilibrium equations of a Navier-Bernoulli-Saint Venant beam, which effective elastic properties are obtained from the solution of the microscopic problem. The reduction of the microscopic problem to the cross-section requires to formulate the homogenization theory in a particular helical coordinate system, called twisting system. In the twisting system, it is shown that bending loads yet depend on the axial coordinate and require a specific treatment leading to separate variable solutions in the axial variable of complex type. Therefore, this paper advances one step further our previous paper [START_REF] Frikha | Mechanical modeling of helical structures accounting for translational invariance. part 1: Static behavior[END_REF] in which the microscopic problem was reduced to a two-dimensional formulation with the drawback that only axial loads (extensional or torsional) could be considered.

This paper is organized as follows. First, Sec. 2 briefly recalls the main results of the homogenization of periodic beam-like structures following an asymptotic expansion method. Section 3 introduces the twisting curvilinear coordinate system. Thanks to this system, Sec. 4 shows how to reduce the initial three-dimensional microscopic problem to a two-dimensional problem posed on the cross-section. Two types of external loads have to be distinguished. Extensional and rotational loads, axially constant, can be treated together. Bending loads, which are shown to depend on the axial coordinate in the twisting system, require a separate treatment. The microscopic problem is then solved with a two-dimensional FE method. From the microscopic solution, Sec. 5 shows how to post-process the macroscopic stiffness of the homogenized beam. Finally, Sec. 6 presents some numerical results and validation for cylinders, springs, and seven-wire strands.

Background: homogenization of beam-like structures

This section recalls the main results concerning the homogenization of periodic beam-like structures based on an asymptotic expansion method. Details can be found in Ref. [START_REF] Buannic | Higher-order effective modelling of periodic heterogeneous beams -Part 1: Asymptotic expansion method[END_REF] for instance. More about asymptotic expansion methods for slender structures may be found in [START_REF] Sanchez-Hubert | Introduction aux méthodes asymptotiques et l'homogénéisation[END_REF][START_REF] Trabucho | Mathematical modelling of rods, Handbook of Numerical Analysis IV[END_REF][START_REF] Kalamkarov | Analysis, Design and Optimization of Composite Structures[END_REF].

A periodic beam-like structure exhibits a small parameter corresponding to the inverse of the slenderness ratio (i.e. the ratio between the width of the cross-section and the total length of the structure), and a small parameter corresponding to the ratio of the length of the periodic heterogeneity to the total length. These two parameters will be embedded in a single parameter, denoted as ε, which amounts to assume that they simultaneously tend to zero [START_REF] Buannic | Higher-order effective modelling of periodic heterogeneous beams -Part 1: Asymptotic expansion method[END_REF].

Two scales are hence considered: the macroscopic 1D-variable X 3 (slow scale) and the microscopic variables Y 1 , Y 2 , Y 3 (fast scale), which are related through:

(Y 1 , Y 2 , Y 3 ) = X 1 ε , X 2 ε , X 3 ε (1)
where upper-case letters have been used to denote Cartesian coordinates (rectangular coordinates). Accordingly, upper-case letters will be used to denote the unit vectors of the Cartesian basis, (E 1 , E 2 , E 3 ). Lower-case letters will be left for curvilinear coordinates (helical coordinates), which will be defined later. The beam asymptotic solution depends on X 3 , the axial slow variable, and Y i (i = 1, 2, 3), the fast variables, leading to the following change of variables:

(X 1 , X 2 , X 3 ) -→ (X 3 , Y 1 , Y 2 , Y 3 ) (2)
Then, the displacement solution is searched under the asymptotic expansion form:

u(X 3 , Y) = u 0 (X 3 , Y) + εu 1 (X 3 , Y) + ε 2 u 2 (X 3 , Y) + ... (3)
where the kth-order displacement u k (X 3 , Y) is periodic in Y 3 , with a period denoted as L 3 , the length of the cell at the microscopic scale. Introducing the expansion (3) in the equilibrium equation of linear elasticity and considering X 3 and Y as independent coordinates, the initial three-dimensional elasticity problem is replaced with a family of problems by equating the terms of the same order ε. These problems are posed on the periodic cell, denoted as V Y , thanks to the periodicity of u k , which plays the role of boundary conditions. The P k cell problem consists in finding u Y 3 -periodic such that ∀δu is kinematically admissible (which means smooth enough and periodic in Y 3 with a period L 3 ):

V Y ∇ s Y δu : σ k+1 dV Y = - V Y δu • (∇ X 3 • σ k )dV Y ( 4
)
where σ k+1 = C : ǫ k+1 and ǫ k+1 = ∇ s Y u k+2 + ∇ s X 3 u k+1 are the stress-strain and strain-displacement relationships at the order k + 1, C being the elasticity tensor. The notation ∇ s Y (•) is used for the symmetric strain operators with respect to Y. The notation ∇ X 3 •(•) is used for the divergence operators with partial differentiation with respect to the only variable X 3 . For conciseness, we assume that the peripheral boundary conditions are free so that no boundary condition terms appear in Eq. (4) (the reader may refer to [START_REF] Buannic | Higher-order effective modelling of periodic heterogeneous beams -Part 1: Asymptotic expansion method[END_REF] to include them).

It can be shown that the 0th-order displacement corresponds to a transverse deflection (the transverse macroscopic displacement):

u 0 (X 3 , Y) = U 0 α (X 3 )E α (5)
It has no axial component, which results from the property that the bending stiffness is much lower than the axial stiffness for slender structures. A proof of this result may be found in [START_REF] Trabucho | Mathematical modelling of rods, Handbook of Numerical Analysis IV[END_REF] for homogeneous beams and in [START_REF] Kolpakov | Calculation of the characteristics of thin elastic rods with a periodic structure[END_REF] for periodic beams. Solving the P -1 cell problem yields a 1st-order displacement of the form:

u 1 (X 3 , Y) = -Y α ∂U 0 α (X 3 ) ∂X 3 E 3 + U 1 3 (X 3 )E 3 + Φ 1 (X 3 )(Y 1 E 2 -Y 2 E 1 ) (α = 1, 2) (6)
which involves a transverse deflection term as well as a rigid-body displacement, which is a composition of a translation U 1 3 and a rotation Φ 1 about the axis E 3 . Note that Eq. ( 6) should also include the additional term U 1 α (X 3 )E α . However, this additional term can be omitted for simplicity because it leads to trivial solutions of rigid-body type at the upper order [START_REF] Buannic | Higher-order effective modelling of periodic heterogeneous beams -Part 1: Asymptotic expansion method[END_REF].

From Eq. ( 4) and noticing that σ 0 = 0, the P 0 cell problem consists in solving numerically the following weak variational form:

V Y ∇ s Y δu : C : ∇ s Y u 2 dV Y = - V Y ∇ s Y δu : C : ∇ s X 3 u 1 dV Y ∀δu (7)
The right-hand side is written in term of ∇ s X 3 u 1 , which from Eq. ( 6) depends on the macroscopic strains:

E E = ∂U 1 3 ∂X 3 , E T = ∂Φ 1 ∂X 3 , E C 1 = - ∂ 2 U 0 2 ∂X 2 3 , E C 2 = ∂ 2 U 0 1 ∂X 2 3 ( 8
)
E E is the macroscopic extension, E T is the macroscopic torsion rotation, E C 1 is the macroscopic curvature about the X 1 -axis and E C 2 is the macroscopic curvature about X 2 . The solution to Eq. ( 7) can hence be viewed as a linear combination of these four macroscopic strains, which we will denote as follows for the stress tensor:

σ 1 (X 3 , Y) = σ 1E (Y)E E (X 3 ) + σ 1T (Y)E T (X 3 ) + σ 1C 1 (Y)E C 1 (X 3 ) + σ 1C 2 (Y)E C 2 (X 3 ) (9)
Finally, by considering particular trial displacement fields of rigid-body type, it can be shown that the P 1 cell and P 2 cell problems enable Eq. ( 4) to be written as compatibility equations leading to the following first-order homogenized problem (macroscopic):

                     N 1 ,3 = 0 T 2 α,3 = 0 T 2 α + M 1 α,3 = 0 M 1 3,3 = 0 (α = 1, 2) (10)
where (•) ,3 = ∂/∂X 3 and with the homogenized constitutive law:

               N 1 M 1 3 M 1 1 M 1 2                = [K hom ]                E E E T E C 1 E C 2                (11)
and with the notations:

N 1 (X 3 ) =< Σ 1 33 >, M 1 3 (X 3 ) =< -Y 2 Σ 1 13 + Y 1 Σ 1 23 >, M 1 1 (X 3 ) =< Y 2 Σ 1 33 >, M 1 2 (X 3 ) =< -Y 1 Σ 1 33 >, T 2 α (X 3 ) =< Σ 2 α3 > (12)
where < • >= (1/L 3 ) V Y (•)dV Y and Σ k i j denotes the stress components in the Cartesian frame (E 1 , E 2 , E 3 ):

σ k = Σ k i j E i ⊗ E j ( 13
)
Note that Eqs. [START_REF] Cartraud | Computational homogenization of periodic beam-like structures[END_REF] correspond to the classical equilibrium equation in beam theory (discarding external forces here for simplicity). This homogenized problem has to be fully written in the Cartesian frame to represent the equivalent behavior of a straight beam. The effective stiffness matrix [K hom ] is determined from Eq. ( 9), the solution of the P 0 cell problem, and can be expressed as:

[K hom ] =                 < Σ 1E 33 > < Σ 1T 33 > < Σ 1C 1 33 > < Σ 1C 2 33 > < -Y 2 Σ 1E 13 + Y 1 Σ 1E 23 > < -Y 2 Σ 1T 13 + Y 1 Σ 1T 23 > < -Y 2 Σ 1C 1 13 + Y 1 Σ 1C 1 23 > < -Y 2 Σ 1C 2 13 + Y 1 Σ 1C 2 23 > < Y 2 Σ 1E 33 > < Y 2 Σ 1T 33 > < Y 2 Σ 1C 1 33 > < Y 2 Σ 1C 2 33 > < -Y 1 Σ 1E 33 > < -Y 1 Σ 1T 33 > < -Y 1 Σ 1C 1 33 > < -Y 1 Σ 1C 2 33 >                 (14)
It can be shown that [K hom ] is symmetric. The P 0 cell problem must be computed numerically. For helical structures, this is usually computationally expensive because the periodic cell is three-dimensional and its period L 3 is equal to the helix pitch, usually large compared to the cross-section dimension. The next sections demonstrate how to reduce this computation to the cross-section.

Introduction of the coordinate system for screw symmetry

As explained in the introduction, the reduction of the microscopic problem to the cross-section requires to formulate the homogenization theory in a particular helical coordinate system, called twisting system. This section introduces this system. For an introduction to the use of general curvilinear coordinate systems, tensor calculus, and differential geometry, the reader may refer to [START_REF] Synge | Tensor Calculus[END_REF][START_REF] Chapelle | The Finite Element Analysis of Shells-Fundamentals[END_REF][START_REF] Gray | Modern Differential Geometry of Curves and Surfaces with Mathmatica, 3rd Edition, Chapman and hall[END_REF] for instance.

Twisting coordinate frame

The relationship between the rectangular coordinates (Y 1 , Y 2 , Y 3 ) and the twisting coordinates (y 1 , y 2 , y 3 ) can be written as [START_REF] Frikha | Mechanical modeling of helical structures accounting for translational invariance. part 1: Static behavior[END_REF]:

           Y 1 = y 1 cos τy 3 -y 2 sin τy 3 Y 2 = y 1 sin τy 3 + y 2 cos τy 3 Y 3 = y 3 ( 15
)
where τ = 2π/L 3 is the torsion of the coordinate system and L 3 denotes the helix pitch of peripheral wires. A twisting system indeed corresponds to a particular case of helical coordinate system, having non-zero torsion but zero curvature. With this kind of system, the cross-section plane (y 1 , y 2 ) remains perpendicular to the straight axis Y 3 but rotates around this axis by following the helical nature of the geometry (see Fig. 1 for the example of a seven-wire strand). Therefore, the cross-section of a screw-symmetric structure will remain invariant in this system. Note that this kind of system has also been used for the analysis of helical waveguides in elastodynamics [START_REF] Onipede | Propagating waves and end modes in pretwisted beams[END_REF][START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF][START_REF] Treyssède | Mode propagation in curved waveguides and scattering by inhomogeneities: application to the elastodynamics of helical structures[END_REF] and optics [START_REF] Nicolet | Finite element analysis of helicoidal waveguides[END_REF].

A key point here is that the change of variables from rectangular to twisting coordinates is applied to the fast variables Y i instead of the slow variables X i . This means that we allow the torsion to operate at the microscopic scale. Conversely, replacing Y i with X i into Eq. ( 15) would have led to the implicit assumptions that τ is small (as adopted in Ref. [START_REF] Nicolet | High order asymptotic analysis of twisted electrostatic problems[END_REF] for instance).

To summarize, from Eqs. ( 2) and ( 15), two changes of variables are applied in this paper:

(X 1 , X 2 , X 3 ) -→ (X 3 , Y 1 , Y 2 , Y 3 ) -→ (X 3 , y 1 , y 2 , y 3 ) (16)
An important remark is that the rectangular coordinate X 3 is kept for the slow scale because the effective behavior targeted by the homogenization process is that of a straight beam. Therefore, the final change of variable given by Eq. ( 16) combines helical coordinates (for the fast scale) with a rectangular coordinate (for the slow scale). Formally, applying the chain rule leads to the differentiation from the initial to the final variables:

∂ ∂X i = ∂ ∂X 3 δ i3 + 1 ε J i j ∂ ∂y j (i = 1, 2, 3) (17)
where δ i j denotes the Kronecker symbol and J i j = ∂y j /∂Y i denotes the Jacobian of the twisting system transformation. From Eq. ( 3) and ( 16), the asymptotic solution will now depend on the slow rectangular scale X 3 and the twisting fast scale

y i (i = 1, 2, 3): u(X 3 , y) = u 0 (X 3 , y) + εu 1 (X 3 , y) + ε 2 u 2 (X 3 , y) + ... ( 18
)
where the kth-order displacement u k (X 3 , y) is periodic in y 3 (with period equal to L 3 since y 3 = Y 3 ). Furthermore, the Serret-Frenet basis associated with the twisting system, denoted as (e 1 , e 2 , e 3 ), is given in terms of the rectangular basis (E 1 , E 2 , E 3 ) by [START_REF] Frikha | Mechanical modeling of helical structures accounting for translational invariance. part 1: Static behavior[END_REF][START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF][START_REF] Treyssède | Mode propagation in curved waveguides and scattering by inhomogeneities: application to the elastodynamics of helical structures[END_REF]: [START_REF] Karathanasopoulos | Two dimensional modeling of helical structures, an application to simple strands[END_REF] and conversely:

           e 1 = E 1 cos τy 3 + E 2 sin τy 3 e 2 = -E 1 sin τy 3 + E 2 cos τy 3 e 3 = E 3
           E 1 = e 1 cos τy 3 -e 2 sin τy 3 E 2 = e 1 sin τy 3 + e 2 cos τy 3 E 3 = e 3 ( 20
)
The unit vectors e 1 , e 2 and e 3 are the normal, binormal and tangent vectors respectively [START_REF] Gray | Modern Differential Geometry of Curves and Surfaces with Mathmatica, 3rd Edition, Chapman and hall[END_REF]. The vector and tensor components expressed in the twisting coordinate frame will be denoted using lower-case symbols (upper-case symbols being left for components in the Cartesian frame):

u k = u k i e i , ǫ k = ǫ k i j e i ⊗ e j , σ k = σ k i j e i ⊗ e j ( 21
)
In particular, from Eqs. ( 13), ( 21) and ( 19), the following identities hold:

Σ 1 33 = σ 1 33 , Σ 1 13 = σ 1 13 cos τy 3 -σ 1 23 sin τy 3 , Σ 1 23 = σ 1 13 sin τy 3 + σ 1 23 cos τy 3 ( 22
)
These relations will be useful later in Sec. 5. At this stage, the formulation of the microscopic problem is still three-dimensional. As shown in the rest of the paper, the solution of this problem will turn out to be separable in terms of the y 3 -variable, hence leading to a two-dimensional model.

Key property: separation of variables

The key property of the microscopic problem considered in this paper is the separation of mechanical fields in terms of the y 3 -variable. Let us recall that the separation of a variable is applicable when the coefficients of the governing partial differential equations, including boundary conditions, are independent of that variable. This implies that the material properties have to be independent of y 3 , which will be assumed throughout this paper.

Furthermore, to apply a separation of the y 3 variable in our problem, the equilibrium equations have to be written in an appropriate curvilinear coordinate system capable of tracking the geometry of the structure (otherwise, the coefficients of the boundary conditions will depend on y 3 ), that is, the twisting frame here. The consequence of this change of coordinate system is that the coefficients of the partial differential operators become also dependent on the geometry through the Jacobian (see Eq. ( 17)). These coefficients can be conveniently written in terms of the so-called Christoffel symbols (see e.g. [START_REF] Synge | Tensor Calculus[END_REF][START_REF] Chapelle | The Finite Element Analysis of Shells-Fundamentals[END_REF]).

First, let us calculate the metric tensor, defined by (g) i j = g i • g j where the vectors (g i ) k = ∂Y k /∂y i form the covariant basis. From Eq. ( 15), the following expression is obtained for the metric tensor of the twisting coordinate system:

g =           1 0 -τy 2 0 1 τy 1 -τy 2 τy 1 τ 2 (y 2 1 + y 2 2 )           . ( 23
)
The Christoffel symbols, denoted as Γ k i j , are defined as follows [START_REF] Gray | Modern Differential Geometry of Curves and Surfaces with Mathmatica, 3rd Edition, Chapman and hall[END_REF]:

Γ k i j = ∂g i ∂x j • g k (24)
They can be expressed as a function of the metric tensor solely by [START_REF] Gray | Modern Differential Geometry of Curves and Surfaces with Mathmatica, 3rd Edition, Chapman and hall[END_REF]:

Γ k i j = 1 2 g kl ∂g jl ∂x i + ∂g il ∂x j - ∂g i j ∂x l ( 25
)
where g i j = g i • g j is the contravariant metric tensor, equal to the inverse of the covariant metric tensor (g i j = (g -1 ) i j ).

As noticed from Eq. ( 23), the independence of g on y 3 is fulfilled by the twisting coordinate system. The coefficients of the differential operators expressed in the twisting frame will depend only on the cross-section variables (y 1 , y 2 ) and the torsion τ (see e.g. the operator L y , given by Eq. ( 28)). The separation of the y 3 variable in the twisting frame is hence justified, somehow a posteriori, from the calculation of the metric tensor.

This key property will be exploited in the next section. It will be shown that the microscopic strains, representing the external forces (right-hand side) of the microscopic problem, admit a specific separation of the y 3 -variable, thus allowing the solution itself to be of the same form.

Solution of the microscopic problem

Transformation of the microscopic problem

The P 0 cell problem given by Eq. ( 7), written with rectangular variables, is transformed into the twisting coordinate system.

For convenience, the notations u k , ǫ k and σ k will now represent the components of fields in the orthonormal basis of the twisting system, (e 1 , e 2 , e 3 ), using Voigt notations for the second-order tensors, so that:

u k =          u k 1 u k 2 u k 3          , ǫ k =                          ǫ k 11 ǫ k 22 ǫ k 33 2ǫ k 12 2ǫ k 13 2ǫ k 23                          , σ k =                          σ k 11 σ k 22 σ k 33 σ k 12 σ k 13 σ k 23                          (26)
The constitutive law is now written as σ k = Cǫ k , where C is the 6-by-6 constitutive matrix of material properties.

Based on these notations, it can be shown that the microscopic strain ∇ s y u 2 and the macroscopic strain ∇ s X 3 u 1 are expressed in the twisting frame as:

∇ s y u 2 = L y u 2 , ∇ s X 3 u 1 = L X 3 ∂u 1 ∂X 3 ( 27
)
with:

L y =                              ∂ ∂y 1 0 0 0 ∂ ∂y 2 0 0 0 τy 2 ∂ ∂y 1 -τy 1 ∂ ∂y 2 + ∂ ∂y 3 ∂ ∂y 2 ∂ ∂y 1 0 τy 2 ∂ ∂y 1 -τy 1 ∂ ∂y 2 + ∂ ∂y 3 -τ ∂ ∂y 1 τ τy 2 ∂ ∂y 1 -τy 1 ∂ ∂y 2 + ∂ ∂y 3 ∂ ∂y 2                              , L X 3 =                          0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0                          (28)
Details about the calculation of the operator L y and the matrix L X 3 are given in Appendix A. Therefore, the transformation of Eq. ( 7) into the twisting coordinate system leads to the weak variational form:

V y (L y δu) T CL y u 2 dV y = - V y (L y δu) T CL X 3 ∂u 1 ∂X 3 dV y ∀δu (29)
for any kinematically admissible displacement δu. Note that from Eq. ( 23), one has det g = 1, so that: dV Y = det g dV y = dV y , where dV Y = dY 1 dY 2 dY 3 and dV y = dy 1 dy 2 dy 3 .

The right-hand-side of ( 29) is an external force representing the action of prescribed macroscopic strains, involved in the vector ∂u 1 ∂X 3 . As shown in Appendix B, this vector can be expressed in the twisting frame as:

∂u 1 ∂X 3 =          -y 2 E T y 1 E T E E + (y 2 cos τy 3 + y 1 sin τy 3 )E C 1 + (y 2 sin τy 3 -y 1 cos τy 3 )E C 2          (30)
In the above expression, two kinds of terms can be distinguished: those independent of y 3 , related to E E,T , and those dependent of y 3 , related to E C α (α = 1, 2). As shown in the next subsections, the treatment of these terms is different.

Solution without bending

Without bending (E C α = 0), the vector ∂u 1 ∂X 3 does not depend on y 3 :

∂u 1 ∂X 3 =          -y 2 E T y 1 E T E E          (31)
Owing to the invariance of the geometry in the twisting frame (outlined in Sec. 3.2), the mechanical fields do not longer depend on the microscopic coordinate y 3 . Therefore, the problem becomes is fully y 3 -independent so that V y (•)dV y = S y (•)dS y ( dy 3 ), where dS y = dy 1 dy 2 and the weak variational form (29) can be reduced to the twodimensional cross-section S y .

A finite element approximation of the form u 2 = N (e) U 2(e) is introduced, where N (e) is the matrix of twodimensional shape functions and U 2(e) the nodal displacements with three degrees of freedom at each node. The weak variational form (29) reduced on the cross-section yields the following linear matrix system:

KU 2 = -F 1 ( 32
)
with the elementary matrices:

K (e) = S (e) y (L ⊥ y N (e) ) T CL ⊥ y N (e) dS y F 1(e) = S (e) y (L ⊥ y N (e) ) T CL X 3          -y 2 E T y 1 E T E E          dS y ( 33
)
where T denotes the matrix transpose, the FE solution U 2 corresponds to the discretized microscopic displacement u 2 and where L ⊥ y denotes the operator L y discarding the derivatives ∂/∂y 3 . The solution σ 1 to this problem is written:

σ 1 (X 3 , y 1 , y 2 ) = σ 1E (y 1 , y 2 )E E (X 3 ) + σ 1T (y 1 , y 2 )E T (X 3 ) (34)
which does not depend on y 3 either. Without bending, it can be checked that the formulation of Refs. [START_REF] Frikha | Mechanical modeling of helical structures accounting for translational invariance. part 1: Static behavior[END_REF][START_REF] Karathanasopoulos | Two dimensional modeling of helical structures, an application to simple strands[END_REF], limited to extensional and torsional loads, is recovered. The goal of this paper is to extend this two-dimensional formulation to loads of bending type.

Solution with bending

With bending loads (E C α 0), the external force ∂u 1 ∂X 3 depends on y 3 . This configuration must be treated separately from the extensional and torsional strains E E and E T . The latter are therefore discarded, so that:

∂u 1 ∂X 3 =          0 0 (y 2 cos τy 3 + y 1 sin τy 3 )E C 1 + (y 2 sin τy 3 -y 1 cos τy 3 )E C 2          (35)
At first sight, the dependence on y 3 prevents the problem to be reduced on the cross-section.

The trick is to take advantage of the separation of variables allowed by the twisting system, as demonstrated in Sec. 3.2, to build solutions having an analytical y 3 -dependence in cos τy 3 and sin τy 3 . In order to get a more compact form, solutions are looked for in terms of the complex exponential e iτy 3 . This kind of dependence is standard in wave problems [START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF][START_REF] Nicolet | Finite element analysis of helicoidal waveguides[END_REF]. It can be justified here in elastostatics due to the helical nature of the geometry, which induces a form of anisotropy in the equilibrium equations. Anisotropy hinders the uncoupling between pure sin(•) and pure cos(•) components of fields (typically only applicable in the full isotropic case).

First, let us rewrite the last component of Eq. [START_REF] Nicolet | High order asymptotic analysis of twisted electrostatic problems[END_REF] as:

(y 2 cos τy 3 + y 1 sin τy 3 )E C 1 + (y 2 sin τy 3y

1 cos τy 3 )E C 2 = Re (y 2 -iy 1 )(E C 1 -iE C 2 )e iτy 3 ( 36
)
This shows that, in the complex space, ∂u 1 /∂X 3 obeys a separation of the y 3 -variable in e iτy 3 . Thanks to the consideration of Sec. 3.2, this separation is also applicable for u 2 . Furthermore, in order to reduce the problem to the cross-section, a convenient choice for the arbitrary test field δu is to adopt a dependence in e -iτy 3 . To summarize, we have:

u 1 (y 1 , y 2 , y 3 ) = u 1 (y 1 , y 2 )e iτy 3 , u 2 (y 1 , y 2 , y 3 ) = u 2 (y 1 , y 2 )e iτy 3 , δu(y 1 , y 2 , y 3 ) = δu(y 1 , y 2 )e -iτy 3 (37)
so that L y becomes:

L y =                              ∂ ∂y 1 0 0 0 ∂ ∂y 2 0 0 0 τy 2 ∂ ∂y 1 -τy 1 ∂ ∂y 2 + iτ ∂ ∂y 2 ∂ ∂y 1 0 τy 2 ∂ ∂y 1 -τy 1 ∂ ∂y 2 + iτ -τ + ∂ ∂y 1 τ τy 2 ∂ ∂y 1 -τy 1 ∂ ∂y 2 + iτ ∂ ∂y 2                              (38)
From Eqs. ( 36)- [START_REF] Chandrasekaran | Finite element analysis of Hertz contact problem with friction[END_REF], the discretization of the weak variational form [START_REF] Chapelle | The Finite Element Analysis of Shells-Fundamentals[END_REF] can now be reduced to the cross-section (note that the integrands do no longer depend on y 3 thanks to the choice of Eq. ( 37)), yielding a linear system KU 2 = -F 1 with:

K (e) = S (e) y (L y N (e) ) * CL y N (e) dS y , F 1(e) = S (e) y (L y N (e) ) * CL X 3          0 0 (y 2 -iy 1 )(E C 1 -iE C 2 )          dS y ( 39
)
where * denotes the matrix transpose conjugate operation. As opposed to Eq. ( 34), the stress σ 1 now depends on y 3 through the separate variable form:

σ 1 (X 3 , y 1 , y 2 , y 3 ) = σ 1C 1 (y 1 , y 2 )E C 1 (X 3 )e iτy 3 + σ 1C 2 (y 1 , y 2 )E C 2 (X 3 )e iτy 3 (40)
From Eqs. [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF], the solution to the external force ( 35) can be obtained by taking the real part of the complex solution. The microscopic displacement field solution to E C α 0 is hence given by Re(U 2 e iτy 3 ).

Interestingly, let us assume

E C 1 = E C 2 = E C .
Noticing that Re (y 2iy 1 )(-iE C )e iτy 3 = Im (y 2iy 1 )E C e iτy 3 and introducing it into Eq. ( 36), it can be deduced that the imaginary part of the complex solution, Im(U 2 e iτy 3 ), computed for (E C 1 = E C , E C 2 = 0) indeed corresponds to the real part of the solution computed for (E C 1 = 0, E C 2 = E C ). Therefore, the solution to the bending problem can be computed only once by exploiting both the real and the imaginary parts of the complex solution U 2 e iτy 3 . In particular, it can be checked that the following relationship holds: Im{σ 1C 1 (y 1 , y 2 )e iτy 3 } = Re{σ 1C 2 (y 1 , y 2 )e iτy 3 }, Im{σ 1C 2 (y 1 , y 2 )e iτy 3 } = -Re{σ 1C 1 (y 1 , y 2 )e iτy 3 } (41)

Determination of the macroscopic stiffness

Stress relationships between Cartesian and curvilinear components

Once the microscopic problem is solved, the stresses are computed thanks to Eqs. [START_REF] Nicolet | Finite element analysis of helicoidal waveguides[END_REF] or [START_REF] Treyssède | Dispersion curve veering of longitudinal guided waves propagating inside prestressed seven-wire strands[END_REF]. These stresses are given as components in the curvilinear coordinate system and must be transformed to Cartesian components in order to compute the homogenized stiffness matrix [K hom ] defined in Sec. 2 and given by Eq. [START_REF] Erdonmez | A finite element model for independent wire rope core with double helical geometry subjected to axial loads[END_REF].

First, notice that:

< • >= 1 L 3 V Y (•)dV Y = 1 L 3 L 3 0 S y (•)dS y dy 3 ( 42
)
Since τ = 2π/L 3 , the integral L 3 0 (•)dy 3 vanishes except for the following integrands:

L 3 0 dy 3 = L 3 , L 3 0 cos 2 τy 3 dy 3 = L 3 2 , L 3 0 sin 2 τy 3 dy 3 = L 3 2 (43)
Owing to Eqs. ( 15), ( 22), ( 34), ( 40) and (43), it can be checked that the components of [K hom ] are given by:

< Σ 1E,T 33 >= S y σ 1E,T 33 dS y , < Σ 1C α 33 >= 0, < -Y 2 Σ 1E,T 13 + Y 1 Σ 1E,T 23 >= S y (-y 2 σ 1E,T 13 + y 1 σ 1E,T 23 )dS y , < -Y 2 Σ 1C α 13 + Y 1 Σ 1C α 23 >= 0, < Y α Σ 1E,T 33 >= 0, < Y 2 Σ 1C α 33 >= 1 2 Re S y (y 2 + iy 1 )σ 1C α 33 dS y , < -Y 1 Σ 1C α 33 >= -1 2 Im S y (y 2 + iy 1 )σ 1C α 33 dS y (α = 1, 2) (44)
It can be noticed that there is no coupling between bending and stretch-torsion (K hom 13 = K hom 14 = K hom 23 = K hom 24 = 0). Furthermore, using Eqs. ( 41), the following identities can be obtained:

< -Y 1 Σ 1C 2 33 >=< Y 2 Σ 1C 1 33 >, < -Y 1 Σ 1C 1 33 >= -< Y 2 Σ 1C 2 33 > (45)
This proves that, for any continuously screw-symmetric structure, the homogenized bending stiffnesses in both directions are equal (K hom 33 = K hom 44 ), and that no coupling occurs in bending (K hom 34 = -K hom 43 = 0 owing to the symmetry of [K hom ]). The homogenized stiffness matrix has therefore the following structure:

[K hom ] =               K hom 11 K hom 12 0 0 K hom 22 0 0 sym. K hom 33 0 0 K hom 33               (46)

Post-processing the resultant forces

Following Ref. [START_REF] Cartraud | Computational homogenization of periodic beam-like structures[END_REF], the components of [K hom ] can be readily post-processed from the total strain energy per unit length:

E p = 1 2                E E E T E C 1 E C 2                T [K hom ]                E E E T E C 1 E C 2                (47)
First, let us consider the case without bending (E C α = 0). The strain energy can be written as:

E p = 1 2 S y (ǫ 1 ) T σ 1 dS y ( 48
)
where

σ 1 = Cǫ 1 , ǫ 1 = (L ⊥ y u 2 + L X 3 ∂u 1
X 3 ) and ∂u 1 ∂X 3 is defined by Eq. ( 31). Let us introduce the FE vector Φ 1 defined at the element level from:

Φ 1(e) = S (e) y (L X 3 N (e) ) T CL X 3 ∂u 1 ∂X 3 dS y (49)
Accounting for the equilibirum KU 2 = -F 1 (F 1 being defined by Eq. ( 33)), the strain energy (48) simplifies into:

E p = 1 2 (F 1 ) T U 2 + (D 1 ) T Φ 1 (50)
where D 1 is the FE vector artificially built from the field

∂u 1 X 3 = [-y 2 E T y 1 E T E E ] T (i.e. at node (n): D 1(n) = [-y (n) 2 E T y (n) 1 E T E E ] T ).
In case of bending (E C α 0), the fields are complex and the following relation holds:

1 L 3 L 3 0
Re(ǫ 1 e iτy 3 ) T Re(σ 1 e iτy 3 )dy 3 = 1 2 Re(ǫ 1 * σ 1 ) (51)

The strain energy becomes:

E p = 1 4 Re(F 1 * U 2 + D 1 * Φ 1 ) (52)
where F 1 is now defined by Eq. ( 39), the vectors Φ 1 and D 1 are now built from the field ∂u 1 X 3 = [0 0 (y 2iy 1 )(E C 1 -iE C 2 )] T at every node.

As an example, considering the solution for E C 1 = 1 (E E,T = E C 2 = 0), the strain energy as calculated by Eq. (52) leads to K hom 33 /2.

Results

Since the bending behaviour constitutes the main originality of the paper, the numerical results of this paper are focused on bending loads (we recall that in Sec. 4.2, it has been shown that the formulation obtained under extensional or torsional loads is equivalent to that of Ref. [START_REF] Frikha | Mechanical modeling of helical structures accounting for translational invariance. part 1: Static behavior[END_REF], where numerical results can be found). Furthermore, the homogenized bending stiffness has been proved to be equal in both directions (see Sec. 5). Therefore, the curvature E C 2 can be set to zero for simplicity (E C 1 0).

In this paper, the FE meshes are generated by the free software Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF] using six-node triangles (quadratic interpolation). Isotropic materials are considered. The Young's modulus and Poisson ratio are denoted by E and ν respectively.

Validation

Two test cases are considered to validate the twisting formulation of Sec. 4.3. The Poisson ratio is set to ν = 0.3.

Circular straight beam

The first test case corresponds to the bending of a straight circular beam, that is, a cylinder of circular cross-section. The moment-curvature relationship for this problem is well-known and leads to the following bending stiffness:

K hom 33 = EI c ( 53
)
where I c = πr 4 c /4 (r c being the radius of the cross-section). The microscopic displacement can be expressed analytically by [START_REF] Trabucho | Mathematical modelling of rods, Handbook of Numerical Analysis IV[END_REF]:

u 2 (X 3 , Y) = -νY 1 Y 2 E 1 + ν 2 (Y 2 1 -Y 2 2 )E 2 E C 1 + ν 2 (Y 2 1 -Y 2 2 )E 1 + νY 1 Y 2 E 2 E C 2 (54)
With the cross-section centered at the origin (y 1 , y 2 ) = (0, 0), the problem can be analyzed with a twisting system of arbitrary torsion τ because the geometry of the structure remains invariant under the twisting transformation.

The component K hom 33 of the homogenized stiffness matrix, computed from Secs. 4.3 and 5.2, is given in Table 1 for several torsion values. As expected, the computed homogenized bending stiffness is found to be independent of τ. The relative error remains lower than 10 -4 %, demonstrating the accuracy of the two-dimensional formulation for bending loads.

Figure 2 shows the real and imaginary parts of the normalized stress field σ 1 33 /E as well as the superimposed cross-section FE mesh (4743 degrees of freedom), computed for r c = 1, τ = 1, E C 1 = 1 (E C 2 = 0). As theoretically demonstrated by Eq. ( 41), it can be observed that the imaginary part is equal to the bending solution about the y 2 axis, i.e. for E C 2 = 1 (E C 1 = 0).

Figure 3a shows the real part of the microscopic displacement in the y 2 direction. The displacement agrees with the analytical solution shown in Fig. 3b, built from Eq. (54). Figure 3c gives the imaginary part of the microscopic displacement in the y 1 direction: this displacement field indeed corresponds to the solution (real part) obtained for a macroscopic curvature imposed in the other direction (E C 2 = 1), in accordance with the theoretical considerations of Sec. 4.3.

Spring

The second validation test case corresponds to the bending of a spring. The analytical expression of the bending stiffness of a spring can be found in Ref. [START_REF] Costello | Theory of Wire Rope[END_REF] and is given by: where I h = πr 4 h /4 for a circular cross-section of radius r h and φ denotes the helix lay angle of the spring, defined from:

K hom 33 = EI h 2 cos φ 2 + ν sin 2 φ ( 55
) (a) (b) (c)
tan φ = R h τ (56)
with R h the radius of the helix centerline (and τ = 2π/L 3 ). In the numerical model, it is noteworthy that the cross-section lies in the (y 1 , y 2 ) plane of the twisting frame and hence corresponds to the Y 3 = 0 plane cut. Let (x, y) = (x(t), y(t)) be a curve parametrization of the crosssection normal to the helix centerline, t denoting the curve parameter. For a circular section of radius r h , one has: (x(t), y(t)) = (r h cos t, r h sin t), t ∈ [0; 2π]. It can be shown that the cross-section cut by the plane Y 3 =0 is parametrized as follows [START_REF] Frikha | Mechanical modeling of helical structures accounting for translational invariance. part 1: Static behavior[END_REF][START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF][START_REF] Treyssède | Mode propagation in curved waveguides and scattering by inhomogeneities: application to the elastodynamics of helical structures[END_REF]:

Y 1 (t) = (τ -1 tan φ + x(t)) cos (τy(t) sin φ) -y(t) cos φ sin (τy(t) sin φ) Y 2 (t) = (τ -1 tan φ + x(t)) sin (τy(t) sin φ) + y(t) cos φ cos (τy(t) sin φ) (57)
The above parametrization yields the cross-section shape that must be meshed in the two-dimensional formulation. Note that the shape changes with the helix parameters and that, for φ = 0 (straight wire), it degenerates to a circle. Figure 4a gives the normalized bending stiffness K hom 33 /Eπr 4 c computed from the numerical model as a function of the lay angle ranging from 1 • to 80 • (R h /r h = 10). As it can be observed, the numerical model matches the analytical solution (55).

Figures 4b andc depict the real and imaginary parts of the normalized stress field σ 1 33 /E computed for r h = 1, R h = 10, φ = 80 • and E C 1 = 1. As already shown, the imaginary part is equal to the solution for E C 2 = 1. It can be observed that the stress distributions both satisfy the symmetry of the cross-section about the y 1 axis. The FE mesh (3645 degrees of freedom) is superimposed in the figures. This example has a strong lay angle, which clearly illustrates that the cross-section shape, properly parametrized by Eq. (57), is not circular.

Seven-wire strands

The analysis of seven-wire strands under bending loads is now conducted. As shown by the previous results, the two-dimensional formulation expressed in the twisting system allows the consideration of both straight circular wires and helical wires. The twisting formulation is hence also suitable for seven-wire strands.

Typical strands are characterized by clearances between peripheral wires: their diameter is usually smaller than that of the central wire so that they do not contact each other. This is a widespread design criterion to minimize frictional effects and secondary tensile stresses [START_REF] Nawrocki | A finite element model for simple straight wire rope strands[END_REF][START_REF] Foti | Analytical and finite element modelling of the elastic-plastic behaviour of metallic strands under axial-torsional loads[END_REF]. In this paper, we set r h = 0.967r c (R h = 1.967r c ), where the subscripts c and h are used for the central and peripheral wires respectively. Clearance between peripheral wires is geometrically satisfied under the following condition [START_REF] Foti | Analytical and finite element modelling of the elastic-plastic behaviour of metallic strands under axial-torsional loads[END_REF]:

φ < φ max = arccos 3 (1 + r c /r h ) 2 -1 (with r c > r h ) (58)
yielding φ max = 12.1 • . The FE mesh of peripheral wires is generated by taking care of their cross-section shape owing to Eq. (57). The Poisson ratio is set to ν = 0.28.

As far as the interaction between the central wire and the peripheral ones is concerned, two extreme cases are commonly considered [START_REF] Zhang | Finite element solutions to the bending stiffness of a single-layered helically wound cable with internal friction[END_REF]: full slip-state (no friction) and full stick-state (no slip).

Full slip-state

In this case, the wires can freely slip relative to each other. This configuration occurs when the bending curvature is sufficiently high to neglect interwire friction. The bending stiffness reaches its minimum [START_REF] Papailiou | On the bending stiffness of transmission line conductors[END_REF]. The strand can be treated as an assemblage of independent wires, which from Eq. ( 53) and (55) leads to the following analytical solution [START_REF] Costello | Theory of Wire Rope[END_REF]:

(K hom 33 ) min = EI c + 6EI h 2 cos φ 2 + ν sin 2 φ (59)
In the numerical model, the continuity of the displacement in the contact region is enforced only in the radial direction. For simplicity, the contact region is restricted to a single node (this restriction will be removed later in Sec. 6.2.3).

Figure 5a shows the homogenized stiffness K hom 33 computed with the numerical model as a function of the lay angle of peripheral wires, ranging from φ = 1 • to φ = 15 • . Although the peripheral wires are geometrically touching each other from φ = 12.1 • (see Eq. ( 58)), peripheral contact is neglected in the whole angle range for simplicity. The numerical solution is in very close agreement with the analytical solution (59).

Figures 5b andc show the real and imaginary parts of the normalized stress field σ 1 33 /E computed for r c = 1, φ = 7.9 • and E C 1 = 1, as well as the FE mesh used for the computation (12, 369 degrees of freedom). As expected, each wire behaves freely (the axial stress varies linearly inside each wire). Note that the imaginary part clearly corresponds to the bending solution about the y 2 axis (E C 2 = 1), as expected from Eq. ( 41).

Full stick-state

In this case, the entire cable behaves like a single body: the peripheral wires are fully bonded to the central wire. Such a configuration typically occurs at the beginning of bending when the radial contact forces are high enough to prevent slipping. According to Ref. [START_REF] Papailiou | On the bending stiffness of transmission line conductors[END_REF], the bending stiffness is maximum and can be approximated by the following analytical formula:

(K hom 33 ) max = EI c + 6EI h cos φ + 3EA h R 2 h cos 3 φ (60)
where A h = πr 2 h is the area of a peripheral wire cross-section. Compared to Eq. ( 59), the last term in the right hand side can be viewed an additional term to the bending stiffness, when all wires are bonded together as a solid body.

In the numerical model, the continuity of the displacement in the contact region is enforced in the three directions. As previously, the contact region is restricted to a single node.

Figure 6a shows the homogenized stiffness K hom 33 computed with the numerical model as a function of the lay angle of peripheral wires. The numerical solution is in good agreement with the analytical approximation (60) for φ = 0 • but deviates as the lay angle increases. Numerical tests have been conducted by refining the FE mesh and by setting the Poisson ratio ν to zero (as implicitly supposed in the formula (60), showing that the numerical bending stiffness found in Fig. 6a remains nearly unchanged (results not shown for conciseness). These discrepancies with the analytical solution will indeed be explained later in Sec. 6.2.4.

Figures 6b andc give the real and imaginary parts of σ 1 33 /E computed for r c = 1, φ = 7.9

• and E C 1 = 1 (same FE mesh as before, involving 12, 369 dofs). As opposed to the full slip-state, the bending motion of the strand appears of global type and resembles that of a single beam (the stress varies continuously in the whole cross-section).

Effect of contact area due to an applied extension

In the previous simulations, the strand was initially unloaded, that is, no extension was applied (E E = 0). The interwire contact region was restricted to a single node, hence corresponding to the limiting case of a zero contact width. Such a configuration is somehow purely theoretical. A single contact point can indeed be viewed as a singular point. In contact mechanics, it is well known that a minimum number of elements is required to discretize a given contact width (see Refs. [START_REF] Chan | A finite element method for contact problems of solid bodies -Part I. Theory and validation[END_REF][START_REF] Chandrasekaran | Finite element analysis of Hertz contact problem with friction[END_REF] for instance). Therefore, it could be argued that numerical convergence can hardly be achieved for the theoretically unloaded case. In this section, the increase of the contact width due to an applied extension E E 0 is accounted for. The loading process consists in two steps.

First, the strand is gradually loaded up to E E = 0.02 (E T = E C α = 0). A tensile strain tends to increase the interwire contact width because of the helical geometry of peripheral wires, yielding a radial compression of the central wire. These contact phenomena are neglected when the linear static problem given by Eq. ( 32) is solved in a unique step. Instead, an iterative procedure must be used in order to account for the contact evolution (nonlinear by nature). Our modeling approach is based on a node-to-node contact procedure using a direct elimination method [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. A matching mesh is used inside the interfacial zone. For the two-dimensional cross-section FE model used in this paper, the computation starts with a single point contact. As the external load is incremented, contact pairs of nodes are successively formed: once the gap between a pair of nodes is closed, the continuity of displacement is enforced at these nodes by the direct elimination method. The continuity is enforced along the three directions, hence corresponding to the assumption of stick contact conditions (full stick-state), or along the radial direction only for the full slip-state. The reader may refer to [START_REF] Treyssède | Dispersion curve veering of longitudinal guided waves propagating inside prestressed seven-wire strands[END_REF] for additional details on this iterative procedure. In a second step, the bending load is finally applied (E C α 0). Note that the geometrical non-linear effects are not accounted for in this paper (a linear strain-displacement relation has been assumed).

The FE mesh has been designed such that the contact width for E E = 0.02 is discretized with six three-node line elements, resulting in thirteen contact nodes and six load increments. Figures 7a andb depict the evolution of the homogenized bending stiffness computed for φ = 7.9 • and 15 • as the contact half-width, denoted by a, increases up to the extension E E = 0.02. It can be observed that that the bending stiffness remains constant in the full slip-state, and equal to the analytical solution given by Eq. (59) (in the figures, gray solid lines are perfectly superimposed onto gray dashed lines). However in the full stick-state, the bending stiffness tends to increase with the contact width: the contact width indeed tends to stiffen the bending behavior of the strand. This stiffening effect remains moderate for φ = 7.9 • but becomes more significant for 15 • (i.e. as the lay angle increases). Despite this increase, note that the bending stiffness yet remains far from the analytical formula given by Eq. (60).

In order to check our static modeling of contact, Fig. 7c shows the evolution of the contact half-width with the normal contact force N for a given contact interface. These numerical results are compared to the solution calculated from Hertz theory for parallel cylinders, given by the following analytical solution [START_REF] Johnson | Contact Mechanics[END_REF]:

a = 8 π N 1 -ν 2 E 1 1 r c + 1 r h (61)
The numerical results are in agreement with Hertz solution.

Full stick-state with lateral contact

In order to explain the discrepancies with the analytical solution observed in Fig. 6, Fig. 8a shows the numerical results obtained by artificially connecting together the peripheral wires with their nearest neighbouring nodes (case of lateral contact). These artificial contacting nodes are indicated in Fig. 8b (by black bullets). The numerical solution is now in good agreement with the analytical approximation (60) for any value of φ. The difference between Fig. 6a and 8a can be explained by the rotation of peripheral wires, free when they are not connected to each other, but blocked when connected. In the former case, the plane section hypothesis is not fully valid, leading to a reduced bending stiffness. In the latter case, the whole cross-section behaves like a single solid body in accordance with the beam approximation of Ref. [START_REF] Papailiou | On the bending stiffness of transmission line conductors[END_REF] used to derive Eq. (60). Note that the rotational motion of peripheral wires under bending occurs due to their helical nature. Conversely, such a rotational motion does not occur for straight wires, which explains why the analytical and the numerical solutions agree as the lay angle tends to zero in Fig. 6a. These numerical results are consistent with those found in Ref. [START_REF] Zhang | Finite element solutions to the bending stiffness of a single-layered helically wound cable with internal friction[END_REF] and obtained from a full three-dimensional FE model.

Consequently, approximate beam-like solutions of plane section type can significantly overestimate the actual bending stiffness of strands in the full stick-state, and hereby the critical curvature at which the wires can start to slide relative to each other [START_REF] Foti | An analytical approach to model the hysteretic bending behavior of spiral strands[END_REF]. Figure 6b andc give the real and imaginary parts of σ 1 33 /E computed for r c = 1, φ = 7.9

• and E C 1 = 1. It can be observed that the whole cross-section behaves like a solid body, as assumed in beam theory. The stress values are slightly greater with peripheral contact than without, which is consistent with the increase of bending stiffness.

Conclusion

A two-dimensional formulation has been proposed for modeling the bending behavior of helical beam-like structures. Within the framework of homogenization theory and using the asymptotic expansion method, the initial threedimensional problem has been transformed to a microscopic problem reduced to the two-dimensional cross-section of the structure, thus offering a tremendous reduction of computational cost compared to full three-dimensional FE models. This two-dimensional reduction exploits the continuous helical symmetry of the problem by rewriting the equilibrium equations in a helical curvilinear coordinate system (twisting system). Such a coordinate system supports a key property: the separation of solution fields in terms of the axial variable. It has been shown that the treatment of bending loads in the twisting system leads to a specific analytical separation of variables, of complex type, along the beam axis. Once the microscopic problem is solved, the effective elastic properties of the beam can then be readily post-processed from the expression of the total strain energy. Numerical results, including the distribution of microscopic stresses over the cross-section, have been presented for cylinders, springs, and seven-wire strands. The formulation has been validated by comparison with several analytical solutions.

More particularly, the two-dimensional formulation allows a closer examination of mechanical contact effects on the homogenized properties of seven-wire strands. It has been shown that the increase of contact width due to extensional loads has a negligible effect in the full slip-state but tends to increase the homogenized bending stiffness in the full stick-state. It has also been shown that approximate solutions of rigid-body type (relying on plane section hypothesis) can significantly overestimate the actual bending stiffness in the full stick-state configuration by neglecting the rotational motion of peripheral wires. This rotational motion, inexistent for straight peripheral wires, tends to gradually increase with the lay angle and leads to a bending behavior softer than predicted with analytical solutions of plane section type.

A. Strain-displacement operator in the twisting frame

Introducing the slow helical variables (x 1 , x 2 , x 3 ), let us consider the following change of variables: (X 1 , X 2 , X 3 ) -→ (x 1 , x 2 , x 3 ) -→ (x 1 , x 2 , x 3 , y 1 , y 2 , y 3 ) (A.1)

We have: (y 1 , y 2 ,

y 3 ) = x 1 ε , x 2 ε , x 3 ε (A.2)
From Eqs. (1), (A.2) and ( 15), it can be deduced that the torsion operate as τ/ε for the slow variables. This is consistent with the fact that τ operates at the microscale, as already noticed in Sec. 3. Defined from ∇ s x u = L x u, the strain-displacement operator L x expressed in twisting coordinates can be found in [START_REF] Onipede | Propagating waves and end modes in pretwisted beams[END_REF][START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF][START_REF] Treyssède | Mode propagation in curved waveguides and scattering by inhomogeneities: application to the elastodynamics of helical structures[END_REF] and is:

L x =                             ∂ ∂x 1 0 0 0 ∂ ∂x 2 0 0 0 τ ε x 2 ∂ ∂x 1 -τ ε x 1 ∂ ∂x 2 + ∂ ∂x 3 ∂ ∂x 2 ∂ ∂x 1 0 τ ε x 2 ∂ ∂x 1 -τ ε x 1 ∂ ∂x 2 + ∂ ∂x 3 -τ ε ∂ ∂x 1 τ ε τ ε x 2 ∂ ∂x 1 -τ ε x 1 ∂ ∂x 2 + ∂ ∂x 3 ∂ ∂x 2                             (A.3)
From the chain rule:

∂ ∂x i = ∂ ∂x i + 1 ε ∂ ∂y i = ∂X α ∂x i ∂ ∂X α + ∂X 3 ∂x i ∂ ∂X 3 + 1 ε ∂ ∂y i (i = 1, 2, 3 and α = 1, 2) (A.4)
and noticing ∂ ∂X α = 0 (α = 1, 2) for a beam-like homogenization process, together with X 3 = x 3 , we have:

∂ ∂x i = ∂ ∂X 3 δ i3 + 1 ε ∂ ∂y i (i = 1, 2 , 3) (A.5)
where δ i3 denotes the Kronecker's symbol. Using Eqs. (A.5) into (A.3) yields:

L x = L X 3 + 1 ε L y (A.6)
where L X 3 and L y are given by Eqs. [START_REF] Synge | Tensor Calculus[END_REF].

B. Expression of the displacement field in the twisting frame

From Eqs. ( 5), ( 6), ( 15) and [START_REF] Karathanasopoulos | Finite element modeling of the elastoplastic axial-torsional response of helical constructions to traction loads[END_REF], one gets the following expression of the displacement field in the twisting coordinate system, at the zeroth and first orders: Equation ( 30) follows from Eq. (B.2) written using the notation of Eq. ( 26). As a side remark, it could be checked from Eq. ( 28), (B.1) and (B.2) that the zeroth and first order displacement fields satisfy the following equations:

1 ε L y u 0 = 0, L X 3 ∂u 0 ∂X 3 + L y u 1 = 0 (B.3)
which means that the strain at the order 1/ε and ε 0 vanishes. The displacement up to the first order, given in twisting coordinates by Eq.s (B.1) and (B.2), hence asymptotically satisfies the equilibrium equation, which is consistent with the homogenization theory of beam-like structures [START_REF] Buannic | Higher-order effective modelling of periodic heterogeneous beams -Part 1: Asymptotic expansion method[END_REF].

Figure 1 :

 1 Figure 1: Twisting coordinate system (y 1 , y 2 , y 3 ) for the example of a seven-wire strand. As one travels along the straight axis Y 3 = y 3 , the cross-section axes y 1 and y 2 rotate by following the helical geometry of peripheral wires.

5 Table 1 :Figure 2 :

 512 Figure 2: Circular beam test case with: r c = 1, τ = 1, E C 1 = 1. Normalized stress field σ 1 33 /E: (a) real part, (b) imaginary part.

Figure 3 :

 3 Figure 3: Circular beam test case with: r c = 1, τ = 1, E C 1 = 1. Microscopic displacement in the y 2 direction: (a) real part of the numerical solution, (b) analytical solution. (c) Imaginary part of the microscopic displacement in the y 1 direction (numerical).

Figure 4 :

 4 Figure 4: Spring test case. (a) Homogenized bending stiffness as a function of the lay angle computed by the two-dimensional FE formulation (× dashed line) and analytical solution (solid line). Normalized stress field σ 1 33 /E computed for r h = 1, R h = 10, φ = 80 • and E C 1 = 1: (b) real part, (c) imaginary part.

Figure 5 :

 5 Figure 5: Strand test case, full slip-state. (a) Homogenized bending stiffness as a function of the lay angle computed by the two-dimensional FE formulation (× dashed line) and analytical solution (solid line). Normalized stress field σ 1 33 /E computed for r c = 1, φ = 7.9 • and E C 1 = 1: (b) real part, (c) imaginary part.

Figure 6 :

 6 Figure 6: Strand test case, full stick-state. (a) Homogenized bending stiffness as a function of the lay angle computed by the two-dimensional FE formulation (× dashed line) and analytical solution (solid line). Normalized stress field σ 1 33 /E computed for r c = 1, φ = 7.9 • and E C 1 = 1: (b) real part, (c) imaginary part.

Figure 7 :

 7 Figure 7: Strand test case. Homogenized bending stiffness as a function of the contact half-width a/r c computed by the two-dimensional FE formulation for: (a) φ = 7.9 • , (b) φ = 15 • (solid lines: numerical solutions, dashed lines: analytical solutions, gray lines: full slip-state, black lines: full stick-state). (c) Contact half-width as a function of the normal contact force N/Eπr c computed for φ = 7.9 • (×) and 15 • (•) (dashed line: Hertz analytical solution).

Figure 8 :

 8 Figure 8: Strand test case, full stick-state, artificial contact between peripheral wires. (a) Homogenized bending stiffness as a function of the lay angle computed by the two-dimensional FE formulation (× dashed line) and analytical solution (solid line). Normalized stress field σ 1 33 /E computed for r c = 1, φ = 7.9 • and E C 1 = 1: (b) real part, (c) imaginary part. The artificial contact nodes are indicated by the six black bullets in Fig. (b).

1 + y 1 Φ 1 e 2 +

 112 u 0 (X 3 , y) = U 0 1 cos τy 3 + U 0 2 sin τy 3 e 1 + -U 0 1 sin τy 3 + U 0 2 cos τy 3 e 2 (B.1) and:u 1 (X 3 , y) = -y 2 Φ 1 e