Romain Lagrange

Yann Fraigneau

New estimations of the added mass and damping of two cylinders vibrating in a viscous uid, from theoretical and numerical approaches

Keywords: Vibration, Fluid-structure interaction, Fluid forces, Coupling coecients, Added mass, Added damping, Viscosity eect, Stokes number, Penalization method

This paper deals with the small oscillations of two circular cylinders immersed in a viscous stagnant uid. A new theoretical approach based on an Helmholtz expansion and a bipolar coordinate system is presented to estimate the uid forces acting on the two bodies. We show that these forces are linear combinations of the cylinders accelerations and velocities, through viscous uid added coecients. The self-added mass and damping coecients are shown to decrease with both the Stokes number and the separation distance. The cross-added mass and damping coecients tend to increase with the Stokes number and the separation distance. Compared to the inviscid results, the eect of viscosity is to add a correction term which scales as Sk -1/2 . When the separation distance is suciently large, the two cylinders behave as if they were independent and the Stokes predictions for an isolated cylinder are recovered. Compared to previous works, the present theory oers a simple and exible alternative for an easy determination of the uid forces and related added coecients. To our knowledge, this is also the rst time that a numerical approach based on a penalization method is presented in the context of uid-structure interactions for relatively small Stokes numbers, and successfully compared to theoretical predictions.

Introduction

The determination of the uid force acting on an immersed body has been the topic of considerable experimental and theoretical studies, covering a full range of applications, from turbomachinery [START_REF] Furber | Is the weis-fogh principle exploitable in turbomachinery?[END_REF], heat exchangers tube banks [START_REF] Chen | Vibration of nuclear fuel bundles[END_REF][START_REF] Chen | Dynamics of heat exchanger tube banks[END_REF] to biomechanics of plants [START_REF] Langre | Eects of wind on plants[END_REF] or energy harvesting of exible structures [START_REF] Doare | Piezoelectric coupling in energy-harvesting uttering exible plates: linear stability analysis and conversion eciency[END_REF][START_REF] Singh | Energy harvesting from axial uid-elastic instabilities of a cylinder[END_REF][START_REF] Michelin | Energy harvesting eciency of piezoelectric ags in axial ows[END_REF][START_REF] Virot | Coupling between a ag and a spring-mass oscillator[END_REF][START_REF] Eloy | Aeroelastic instability of cantilevered exible plates in uniform ow[END_REF]. Early researches were stimulated by the need of understanding the eect of the inertia of a surrounding uid on the frequency of an oscillating pendulum [10]. Assuming an inviscid uid, [START_REF] Poisson | Sur les mouvements simultanés d'un pendule et de l'air environnant[END_REF][START_REF] Green | Researches on the vibration of pendulums in uid media[END_REF][START_REF] Stokes | On some cases of Fluid Motion[END_REF] showed that the uid make the mass of the pendulum to increase by a factor that depends on the uid density and the geometry of the pendulum. Since these pioneer works, this apparent increase of mass has commonly been referred as the added mass concept. It has been investigated in various experiments [START_REF] Clough | Eects of eathquakes on underwater structures[END_REF][START_REF] Chandrasekaran | Virtual mass of submerged structures[END_REF][START_REF] Stelson | Virtual mass and acceleration in uids[END_REF][START_REF] Keulegan | Forces on cylinders and plates in an oscillating uid[END_REF][START_REF] Garrison | Hydrodynamic loads induced by earthquakes[END_REF][START_REF] Sarpkaya | Separated ow about lifting bodies and impulsive ow about cylinders[END_REF][START_REF] Chandrasekaran | Vibration of submerged structures[END_REF][START_REF] Skop | Add Mass and Damping Forces on Circular Cylinders[END_REF][START_REF] Sarpkaya | Forces on cylinders and spheres in a sinusoidally oscillating uid[END_REF] in which a single body is accelerated in a uid initially at rest. The acceleration of the body induces a uid motion which in returns induces an inertia eect from which an added mass coecient is computed.

The concept of added mass also applies to multiple immersed bodies, although its formulation is more complex as it involves "self-added" and "cross-added" mass coecients. The self-added mass coecient characterizes the force on a body due to its own motion. The cross-added coecient characterizes the uid-coupling force on a stationary body due to the motion of an other body. Considering multiple arrangements, many experimental rigs have been built [START_REF] Yamamoto | Hydrodynamic forces on groups of cylinders[END_REF][START_REF] Yamamoto | Hydrodynamic forces on multiple circular cylinders[END_REF][START_REF] Dalton | Potential ow past a group of circular cylinders[END_REF][START_REF] Chen | Design guide for calculating hydrodynamic mass. part 1: Circular cylindrical structures[END_REF][START_REF] Chen | Added mass and damping of a vibrating rod in conned viscous uids[END_REF][START_REF] Moretti | Hydrodynamic inertia coecients for a tube surrounded by rigid tubes[END_REF][START_REF] Shin | Flow-induced vibration in lmfbr steam generators: A state-of-the-art review[END_REF][START_REF] Chen | Vibration of nuclear fuel bundles[END_REF][START_REF] Chen | Vibration of a group of circular cylinders in a conned uid[END_REF][START_REF] Chen | Vibrations of a group of circular cylindrical structures in a liquid[END_REF][START_REF] Wu | Dynamic analysis of multibody system immersed in a uid medium[END_REF] to obtain precise measurements of these coecients. From a theoretical stand point, the added coecients should be computed from the Navier-Stokes equations. However, in many practical situations, the eects of uid viscosity and compressibility are neglected and a potential theory is carried out. A method of images [START_REF] Hicks | On the motion of two cylinders in a uid[END_REF][START_REF] Greenhill | Functional images in cartesians[END_REF][START_REF] Basset | A Treatise on Hydrodynamics[END_REF][START_REF] Carpenter | On the motion of two cylinders in an ideal uid[END_REF][START_REF] Birkho | Hydrodynamics[END_REF][START_REF] Gibert | Vibration of structures in a static uid medium[END_REF][START_REF] Landweber | Added masses and forces on two bodies approaching central impact in an inviscid uid[END_REF] or a complex analyis based on conformal transformations [START_REF] Wang | Interaction of two circular cylinders in inviscid uid[END_REF][START_REF] Burton | Hydrodynamic forces on two moving discs[END_REF][START_REF] Tchieu | Fluid-structure interaction of two bodies in an inviscid uid[END_REF][START_REF] Scolan | On the use of conformal mapping for the computation of hydrodynamic forces acting on bodies of arbitrary shape in viscous ow. part 2: multi-body conguration[END_REF][START_REF] Crowdy | Analytical solutions for uniform potential ow past multiple cylinders[END_REF][START_REF] Crowdy | A new calculus for two-dimensional vortex dynamics[END_REF] are usually derived to solve the boundary value problem governing the uid potential function. For small amplitude motions not entailing ow separation, the potential theory will accurately give the added mass coecients, and tabulated results are available in the literature for a wide variety of immersed geometries [START_REF]Tables of hydrodynamic mass factors for translational motion[END_REF]. All of the above-mentioned studies have dealt with an ideal uid, whereas the viscous eects may be important for some applications such as bodies relatively close to each other. Considering the small oscillations of a single body in a viscous uid, Stokes [START_REF] Stokes | On the eect of the internal friction of uids on pendulums[END_REF] solved the linearized Navier-Stokes equations and showed that the uid force is a linear combination of two components related to the acceleration of the body and its velocity. The coecients of this linear combination are commonly referred as the viscous added mass and the viscous added damping, respectively. Stokes found that the eect of viscosity is to add to the ideal uid added mass coecient a correction term which depends on the uid mass density and viscosity, the frequency of oscillation, and a characteristic length scale. All of these eects can be regrouped in a single dimensionless number, the Stokes number.

The extension of this work to the case of multiple bodies remains a challenging theoretical problem, mainly due to the viscous boundary conditions to account for. One approach developped in [START_REF] Chen | Vibration of nuclear fuel bundles[END_REF] is to associate to each body a uid potential and a stream function, governed by a Laplace and an Helmholtz equation, respectively. Introducing a polar coordinate system attached to each body, a method of separation of variables is used to expand the uid functions as an innite trigonometric series with unknown coecients. Applying the viscous boundary conditions into each local coordinate system yields a set of linear equations for these unknowns. The number of equations depends on the number of bodies and the number of terms used in the series expansions. In the end, the set of linear equations has to be solved numerically. The two cylinders problem could be solved in this framework, but even for such a restricted number of bodies, the method of [START_REF] Chen | Vibration of nuclear fuel bundles[END_REF] is hardly tractable.

In this paper, we build on our previous work which dealt with ideal uids [START_REF] Lagrange | A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant uid[END_REF] to introduce a exible theoretical method and obtain an estimation of the viscous added coecients. In addition to this theoretical work, we perform some numerical simulations where the immersed boundary conditions are considered with a penalization method. The choice of this approach relies on its eectiveness and simplicity of implementation in CFD codes, without deep modication of the algorithmic structure. The basic idea is to add a forcing term in the Navier-Stokes equation set over the area of the immersed body in order to locally impose the velocity of the body [START_REF] Peskin | The immersed boundary method[END_REF]. The method does not require any mesh update related to the motion of the body, any complex geometrical considerations on the position of the wall in regard to the computational grid or any high order interpolations as done with some other approaches (e.g. ALE methods [START_REF] Loubère | Reale: A reconnectionbased arbitrary-lagrangianâeulerian method[END_REF], cut-cell methods [START_REF] Cheny | Ls-stag method: A new immersed boundary/level-set method for the computation of incompressible viscous ows in complex moving geometries with good conservation properties[END_REF], immersed body methods [START_REF] Gronski | A simple and ecient direct forcing immersed boundary method combined with a high order compact scheme for simulating ows with moving rigid boundaries[END_REF]). In the present work, we actually use a variant method initially proposed by [START_REF] Pasquetti | A pseudo-penalization method for high reynolds number unsteady ows[END_REF], called the pseudo-penalization method, in which disappears the stiness nature of the Navier-Stokes equations due to the forcing term. The penalization and pseudo-penalization methods are particularly ecient in uid problems with moderate or high Reynolds numbers (see e.g. [START_REF] Mittal | Immersed boundary methods[END_REF][START_REF] Kadoch | A volume penalization method for incompressible ows and scalar advection diusion with moving obstacles[END_REF][START_REF] Kolomenskiy | Two-and threedimensional numerical simulations of the clap-ing-sweep of hovering insects[END_REF][START_REF] Schneider | Immersed boundary methods for numerical simulation of conned uid and plasma turbulence in complex geometries : a review[END_REF][START_REF] Minguez | High-order large-eddy simulation of ow over the Ahmed body car model[END_REF][START_REF] Nore | Numerical simulation of the von karman sodium dynamo experiment[END_REF]) but has never been tested in problems with low Reynolds numbers, as considered in the present work. This paper is organized as follows. Section 2 presents the problem and the governing equations for two circular cylinders immersed in a viscous uid at rest. In section 3, we propose a theoretical approach based on an Helmholtz decomposition and a bipolar coordinate system to obtain an approximate solution of the uid problem. We derive expressions for the uid potential and stream functions, from which we compute the uid forces on the cylinders. In section 4 we describe the numerical simulations that we have performed to solve the uid problem. The results of our investigation are presented in section 5. Throughout, we directly compare the theoretical predictions to the numerical simulations. We start with comparing the time evolutions of the uid forces acting on the cylinders, when one is stationary while the other is imposed a sinusoidal vibration. We then analyze the dependance of the uid added coecients with the Stokes number and the separation distance. Some scaling laws are derived in the limit of large Stokes numbers. Finally, section 6 summarizes our ndings. We consider the simple harmonic motions of two rigid circular cylinders C j , (j = 1, 2), with centers O j , radii R j , boundaries ∂C j , immersed in an innite 2D viscous uid domain, as illustrated in gure 1. The angular frequency of the cylinders is Ω and their displacement vectors are U j . The uid is Newtonian, homogeneous, of volume mass density ρ and kinematic viscosity ν. The Navier-Stokes equations and the boundary conditions for the incompressible uid ow

(V, P ) write ∇ • V = 0, (1a) ∂V ∂T + ∇V • V + 1 ρ ∇P -ν∆V = 0, (1b)
V -

dU j dT = 0 on ∂C j , j = {1, 2}. (1c)
The third equation expresses the continuity of velocities at the cylinder boundaries. The uid force acting on C j is the sum of a pressure and a viscous term, and writes

F j = - ∂C j P n j dL j + ρν ∂C j ∇V + (∇V) T • n j dL j . (2)
In this equation, n j is the outward normal unit vector to ∂C j , (∇V) T the transposate tensor of ∇V and dL j an innitesimal line element of integration.

Dimensionless equations

In what follows, we use R 2 and Ω -1 as a characteristic length and time. Introducing t = T Ω, we dene the dimensionless cylinders displacements u * j , uid ow (v * , p * ) and uid force f * j as

U j = U u * j , V = U Ω v * , P = ρU R 2 Ω 2 p * , F j = ρU (R 2 Ω) 2 f * j , (3)
with

U = max (|U 1 |, |U 2 |).
To reduce the number of parameters of the problem we also introduce the rescaled quantities

r = R 1 R 2 , ε = E R 2 , KC = U R 2 , Sk = R 2 2 Ω ν , (4)
as the radius ratio, separation distance, Keulegan-Carpenter number and Stokes number (i.e. vibration Reynolds number), respectively.

Introducing (3) in (1), the dimensionless Navier-Stokes equations write

∇ • v * = 0, ( 5a
) ∂v * ∂t + KC∇v * • v * + ∇p * - 1 Sk ∆v * = 0, ( 5b
) v * - du * j dt = 0 on ∂C j , j = {1, 2}. ( 5c
)
Introducing (3) in (2), the dimensionless uid force acting on C j write

f * j = - ∂C j p * n j dl j + 1 Sk ∂C j ∇v * + (∇v * ) T • n j dl j , (6)
with dl j = dL j /R 2 .

Theoretical approach

In the limit of small oscillations, i.e. KC = o(1), the nonlinear convective term in the Navier-Stokes equations is negligible. Introducing u * j = {e it u j }, v * = {e it v}, p * = {e it p}, the equations (5) rewrite

∇ • v = 0, ( 7a
) iv + ∇p - 1 Sk ∆v = 0, (7b)
v -iu j = 0 on ∂C j , j = {1, 2}, (7c)
with the real part operator.

Helmholtz decomposition

We seek a solution of (7) as a superposition of an irrotational and a divergencefree ow (Helmholtz decomposition)

v = ∇ϕ + ∇ × A, (8)
with ϕ and A = Ae z some unknown potential and stream functions. Introducing this decomposition in [START_REF] Michelin | Energy harvesting eciency of piezoelectric ags in axial ows[END_REF] yields

∆ϕ = 0, (9a)
∇ × (∆A -iSkA) -Sk∇ (iϕ + p) = 0, ( 9b
) ∇ϕ + ∇ × A -iu j = 0 on ∂C j , j = {1, 2}. (9c)
Taking the divergence and the curl of (9b) yields two equations p = -iϕ and ∆A + β 2 A = 0 with β = √ -iSk, (10) from which the pressure and the stream functions can be determined. Let z = x + iy be the complex number whose real and imaginary parts are the cartesian coordinates x and y, measured from the midpoint O of the two cylinders centers, O 1 and O 2 , see gure 2.

Bipolar coordinates

x ε 2 2r O y z ( ) h z ζ = σ τ 2 τ 1 τ ζ 1 O 2 O 1 C 2 C 2π
Let h(z) be the conformal mapping dened as

ζ = σ + iτ = h (z) = i ln z -x B + a z -x B -a , (11)
with x B = (r 2 -1) / (2d) and

a = d 2 -(1 + r) 2 d 2 -(1 -r) 2 2d , d = r + ε + 1. (12)
In [START_REF] Poisson | Sur les mouvements simultanés d'un pendule et de l'air environnant[END_REF], 0 < σ ≤ 2π and τ ∈ R are the real and imaginary parts of ζ, respectively. They are also the bipolar coordinates of a point in the plane (x, y). The images of C 1 and C 2 are the straight lines with ordinates τ 1 and τ 2 given by τ 1 = -sinh -1 (a/r) < 0 and τ 2 = sinh -1 (a) > 0.

The Laplace operator and the uid velocity vector in bipolar coordinates are

∆ϕ = 1 κ στ 2 ∂ 2 ϕ ∂σ 2 + ∂ 2 ϕ ∂τ 2 , ( 14a
) v = 1 κ στ ∂ϕ ∂σ + ∂A ∂τ e σ + ∂ϕ ∂τ - ∂A ∂σ e τ , (14b)
with κ στ = a/[cosh (τ ) -cos (σ)] the Lamé coecient and

the physical basis vectors. The uid equations [START_REF] Eloy | Aeroelastic instability of cantilevered exible plates in uniform ow[END_REF] in the bipolar coordinates system write

∂ 2 ϕ ∂σ 2 + ∂ 2 ϕ ∂τ 2 = 0, (16a)
∂ 2 A ∂σ 2 + ∂ 2 A ∂τ 2 + β 2 κ στ 2 A = 0, (16b)
∂ϕ ∂σ + ∂A ∂τ = (iu jx )g jy -(iu jy )g jx on τ = τ j , j = {1, 2}, (16c) ∂ϕ ∂τ - ∂A ∂σ = (iu jx )g jx + (iu jy )g jy on τ = τ j , j = {1, 2}, (16d)
with g jx = κ στ j e x • e τ j , g jy = κ στ j e y • e τ j . These are 2π periodic functions of σ given by

g jx (σ) = -a cos (σ) cosh (τ j ) -1 (cosh (τ j ) -cos (σ)) 2 = ∞ n=1 g jn cos (nσ), (17a)
g jy (σ) = -a sin (σ) sinh (τ j ) (cosh (τ j ) -cos (σ)) 2 = ∞ n=1 g jn sgn (τ j ) sin (nσ), (17b)
with g jn = -2nae -n|τ j | .

Ad-hoc problem, uid forces and added coecients

Since the problem is linear in u jx and u jy , the uid functions are linear combinations of the form

ϕ = (u 1x ϕ 1x + u 2x ϕ 2x ) + (u 1y ϕ 1y + u 2y ϕ 2y ) , (18a) A = (u 1x A 1x + u 2x A 2x ) + (u 1y A 1y + u 2y A 2y ) . (18b)
The diculty in nding ϕ jα and A jα arises from the fact that the Helmolhz equation (16b) has a variable coecient, κ στ . Instead, we consider the adhoc problem in which κ στ is replaced by some unknown constant k, that will be determined later on. A method of separation of variables is then used to nd the ad-hoc functions ϕ jα and A jα . The boundary conditions (16c), (16d) along with [START_REF] Keulegan | Forces on cylinders and plates in an oscillating uid[END_REF] indicate that ϕ jα and A jα are linear combinations of cos (nσ) and sin (nσ). Introducing these linear combinations in the Laplace and the Helmholtz equations, we also obtain that ϕ jα (resp. A jα ) is a linear combination of cosh (nτ ) and sinh (nτ ) (resp. cosh (lτ ) and sinh (lτ

) with l = n 2 -(βk) 2 .
All in all, the ad-hoc functions write

ϕ = (iu 1x ) ∞ n=1 cos (nσ) ϕ (1) n (τ 1 , τ 2 , l) cosh (nτ ) + ϕ (2) n (τ 1 , τ 2 , l) sinh (nτ ) -(iu 1y ) ∞ n=1 sin (nσ) ϕ (1) n (τ 1 , τ 2 , l) cosh (nτ ) + ϕ (2) n (τ 1 , τ 2 , l) sinh (nτ ) + (iu 2x ) ∞ n=1 cos (nσ) ϕ (1) n (τ 2 , τ 1 , l) cosh (nτ ) + ϕ (2) n (τ 2 , τ 1 , l) sinh (nτ ) + (iu 2y ) ∞ n=1 sin (nσ) ϕ (1) n (τ 2 , τ 1 , l) cosh (nτ ) + ϕ (2) n (τ 2 , τ 1 , l) sinh (nτ ) , (19a)
A = (iu 1x ) ∞ n=1 sin (nσ) A (1) n (τ 1 , τ 2 , l) cosh (lτ ) + A (2) n (τ 1 , τ 2 , l) sinh (lτ ) + (iu 1y ) ∞ n=1 cos (nσ) A (1) n (τ 1 , τ 2 , l) cosh (lτ ) + A (2) n (τ 1 , τ 2 , l) sinh (lτ ) + (iu 2x ) ∞ n=1 sin (nσ) A (1) n (τ 2 , τ 1 , l) cosh (lτ ) + A (2) n (τ 2 , τ 1 , l) sinh (lτ ) -(iu 2y ) ∞ n=1 cos (nσ) A (1) n (τ 2 , τ 1 , l) cosh (lτ ) + A (2) n (τ 2 , τ 1 , l) sinh (lτ ) , (19b)
with ϕ (j) n and A (j) n given in Appendix B. Plugging the Helmholtz decomposition v = ∇ ϕ + ∇ × A and the pressure equation p = -i ϕ in (6) yields the ad-hoc uid forces f

* j = {e it f j }           f 1x f 1y f 2x f 2y           = π ([M ] -i [C])           u 1x u 1y u 2x u 2y           , (20)
with [M ] and [C] the added mass and damping matrices

[M ] =           m (1) self 0 m cross 0 0 m (1) self 0 -m cross m cross 0 m (2) self 0 0 -m cross 0 m (2) self           , [C] =           c (1) self 0 c cross 0 0 c (1) self 0 -c cross c cross 0 c (2) self 0 0 -c cross 0 c (2) self           . ( 21
)
The self-added mass m (j) self and damping c (j) self relate the uid force on C j to its own motion. The cross-added mass m cross and damping c cross relate the uid force on C m to the motion of C j , j = m.

All the uid added coecients in ( 21) are functions of the radius ratio r, the dimensionless separation distance ε and the Stokes number Sk. A general closed-form expression for these coecients is not tractable, but some simplications are possible in particular cases. For example, as Sk → ∞ (inviscid uid), the ow is purely potential, i.e. ( A, c (j) self , c cross ) → (0, 0, 0), and the added mass coecients simplify to

m (1) self → m (1)P OT self = ∞ n=1 4na 2 e 2nτ 1 tanh [n (τ 2 -τ 1 )] as Sk → ∞, (22a)
m (2) self → m (2)P OT self = ∞ n=1 4na 2 e -2nτ 2 tanh [n (τ 2 -τ 1 )] as Sk → ∞, (22b)
m cross → m P OT cross = ∞ n=1 -4na 2 e -n(τ 2 -τ 1 ) sinh [n (τ 2 -τ 1 )] as Sk → ∞. (22c)
For the sake of clarity, we have reported the study of the variations of m (j)P OT self and m P OT cross in appendix A.

Determination of the ad-hoc constant k

In the previous section, we have obtained solutions of an ad-hoc problem in which the Lamé coecient κ στ has been replaced by some constant k. As a result, the ad-hoc functions A, ϕ and p do not satisfy the Navier-Stokes equation (9b), leading to a non zero local residual

W = u 1x W 1x + u 2x W 2x + u 1y W 1y + u 2y W 2y , (23)
with

W jα = ∇ × ∆ A jα -iSk A jα and A jα = A jα e z .
The constant k is determined from the condition that the weigthed residual

W = 2π 0 τ 2 τ 1 (|W 1x | w 1x + |W 2x | w 2x + |W 1y | w 1y + |W 2y | w 2y )κ στ 2 dτ dσ, (24)
must vanish for some given weight functions w jα . In this study, we consider two families of weight functions, which yields two sets of ad-hoc functions. In the least squares method, the weight functions are chosen in the form

w jα = d dk |W jα | , (25)
such that the residual W vanishes when

χ (k) = 2π 0 τ 2 τ 1 |W 1x | 2 + |W 2x | 2 + |W 1y | 2 + |W 2y | 2 κ στ 2 dτ dσ, (26)
is minimum. We call χ LS this minimum, reached for k = k LS .

In the collocation method, the residual W is forced to vanish on the cylinders boundaries. The weight functions are chosen to be the Dirac functions δ

w jα = d dk |W jα | (δ (τ -τ 1 ) + δ (τ -τ 2 )) , (27)
such that the residual W vanishes when

χ (k) = 2π 0 |W 1x | 2 + |W 2x | 2 + |W 1y | 2 + |W 2y | 2 κ στ 2 (σ, τ 1 ) dσ + 2π 0 |W 1x | 2 + |W 2x | 2 + |W 1y | 2 + |W 2y | 2 κ στ 2 (σ, τ 2 ) dσ, (28)
is minimum. We call χ COL this minimum, reached for k = k COL .

The evolutions of k LS , k COL , χ LS and χ COL , versus the Stokes number Sk are shown in Fig. 3, for equal size cylinders (r = 1) and three dimensionless separation distances ε = {0.5, 1, 2}. We nd that both k LS and k COL decrease with Sk, increase with ε, but remain close to 1. This can be explained from the fact that the bipolar coordinates (σ, τ ) are conformally equivalent to the cartesian coordinates (x, y), in which the Helmholtz equation is similar to (16b) under the change (σ, τ, κ στ ) → (x, y, 1). The numerical method to solve the Navier-Stokes equations ( 5) is based on the projection method of [START_REF] Guermond | An overview of projection methods for incompressible ows[END_REF] and the delta formulation of [START_REF] Goda | A multistep technique with implicit dierence schemes for calculating two or three dimensional cavity ows[END_REF]. The equations are discretized following a nite volume approach on a staggered structured grid (MAC procedure) with a second order approximation in time and space.

A dierentiation forumula (BDF2) is used for the time discretization of (5b), leading to

3v * (n+1) 2δt + -4v * (n) + v * (n-1) 2δt + NL (n+1) + ∇p * (n) - 1 Sk ∆v * (n+1) = 0, (29)
with NL = KC ∇v * • v * and n the subscript for the time step. The convective term at time (n+1)δt is computed from a linear extrapolation of the estimated values at time nδt and (n -1)δt, i.e. 1) . The space discretization of the convective and viscous terms are approximated with a second order centered-scheme. An implicit discretization is applied to the viscous term in order to increase the numerical stability. The pressure gradient is explicitly dened, as suggested in the projection method.

NL (n+1) = 2NL (n) -NL (n-
Introducing δv * i (n+1) = v * i (n+1) -v * i (n)
as the time increment of the i-th component of the velocity vector v * , the equation ( 29) reduces to a Helmholtz equation

δv * i (n+1) - 2 3 δt Sk ∆ δv * i (n+1) = S (n,n-1) i , (30)
where S (n,n-1) i contains all the explicit terms of [START_REF] Shin | Flow-induced vibration in lmfbr steam generators: A state-of-the-art review[END_REF]. The equation ( 30) is solved by means of an Alternating Direction Implicit method, see [START_REF] Peaceman | The numerical solution of parabolic and elliptic dierential equations[END_REF].

The Helmholtz decomposition of v * (n+1) with a potential function Φ yields the two equations

∆Φ = ∇ • v * (n+1) δt and Φ = p * (n+1) -p * (n) - 1 Sk ∇ • v * (n+1) . (31)
The Poisson's equation is solved using a direct method based on the partial diagonalization of the Laplace operator. Having determined Φ, the pressure at time (n + 1)δt is computed from the second equation of [START_REF] Chen | Vibrations of a group of circular cylindrical structures in a liquid[END_REF]. Finally, the velocity eld v * (n+1) is corrected in order to satisfy the divergence-free condition

v * (n+1) := v * (n+1) - 3 2 δt∇Φ. (32)

The pseudo penalization method

The pseudo penalization method is based on the standard volume penalty method, see [START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Mittal | Immersed boundary methods[END_REF][START_REF] Kadoch | A volume penalization method for incompressible ows and scalar advection diusion with moving obstacles[END_REF], and has shown to be eective in solving uid-structure interaction problems involving moving bodies, see [START_REF] Pasquetti | A pseudo-penalization method for high reynolds number unsteady ows[END_REF][START_REF] Nore | Numerical simulation of the von karman sodium dynamo experiment[END_REF]. The principle is to solve some penalized Navier-Stokes equations over a single domain, instead of considering two separate domains (uid and solid) interacting through a set of boundary conditions. The original contribution of [START_REF] Pasquetti | A pseudo-penalization method for high reynolds number unsteady ows[END_REF] relies on the removal of specic terms in the Navier-Stokes equations in order to turn them into steady penalized Stokes equations in the solid domains, where the penalty term is directly provided by the time-discretization scheme.

The penalization of (29) writes

3v * (n+1) 2δt + (1 -χ) -4v * (n) + v * (n-1) 2δt + NL (n+1) + ∇p * (n) - 1 Sk ∆v * (n+1) = 0, (33)
with χ a penalty function dened as χ = 1 in the solid domains and χ = 0 in the uid domain. In [START_REF] Hicks | On the motion of two cylinders in a uid[END_REF], 3v * (n+1) /(2δt) can be seen as a forcing term that makes v * to tend to zero in the solid domains. Although v * does not strictly vanishes in the solid domains, the consistency of the method scales as δt/Sk. Since the forcing term is provided by the time step, 3/(2δt), it does not aect the stiness of the equations, preventing spurious eects or stability constraints, unlike the standard penalization methods.

For a body moving with a velocity v 0 * , (33) can be reformulated as

3v * (n+1) 2δt + (1 -χ) -4v * (n) + v * (n-1) 2δt + NL (n+1) + ∇p * (n) - 1 Sk ∆v * (n+1) = χ 3v 0 * 2δt , (34)
and solved with the numerical method mentionned in section 4.1.

Results

We now present the results of our predictions, considering the case in which C 1 is stationary while C 2 is imposed a sinusoidal displacement in the x -direction.

For the geometric parameters, we have investigated the case of two equal size cylinders, corresponding to a radius ratio r = 1. Three representative values were chosen for the dimensionless separation distance (depicted in the insets of gures 4, 5 and 6): a small gap, ε = 0.5; a gap with size one radius, ε = 1; and a large gap, ε = 2. In the presentation of our results, we rst consider the eect of the Stokes number 100 ≤ Sk ≤ 900 and the dimensionless separation distance on the time evolution of the uid forces. We then analyze the evolution of the magnitude h j and phase φ j of the forces, including the case ε → ∞ for which Stokes [START_REF] Stokes | On the eect of the internal friction of uids on pendulums[END_REF] obtained

f 2x = π m ISO self -ic ISO self u 2x = π   1 + 4 √ iSk K1 √ iSk K0 √ iSk   u 2x , (35)
with K 0 and K 1 the modied Bessel functions of second kind. We nally study the evolution of the uid added coecients and derive some scaling laws for large Stokes numbers. Throughout the study, we perform some numerical simulations to corroborate the theoretical predictions, also providing a discussion on the limitations of both approaches.

Theoretical predictions

Since the problem is symetric about the axis τ = 0, we have

τ 1 = -τ 2 , m (1)
self = m (2) self = m self , c (1)
self = c (2)
self = c self and m (1)P OT self = m

(2)P OT self = m P OT self . The dimensionless ad-hoc uid forces are computed from [START_REF] Chandrasekaran | Vibration of submerged structures[END_REF], with u 1 = 0, u 2x = -i and u 2y = 0, leading to

f 1x = π (m cross -ic cross ) u 2x = h 1 e iφ 1 u 2x , (36a)
f 2x = π (m self -ic self ) u 2x = h 2 e iφ 2 u 2x . (36b)

Numerical setup

Concerning the numerical simulations, the computational domain size S is considered suciently large to minimize the end eects. For the small and medium separation distances (ε = 0.5 and ε = 1), we set S = 20 × 17. For ε = 2, we set S = 22 × 17 so that the distance between the cylinders and the domain ends is similar to the cases ε = 0.5 and ε = 1. For all the simulations, the Keulegan-Carpenter number is set to KC = 10 -2 . The cartesian grid is built with a regular distribution over the cylinder domains, including the displacement zone. The dimensionless cell size is 2 × 10 -3 in both the x and y directions. It follows that the smallest spatial scale of our problem, i.e. the cylinder displacement, is discretized over ten square cells, which yields a satisfying spatial resolution. The cell-size distribution outside the cylinder domain is performed with a hyperbolic tangent function and vary from 2×10 -3 to 3.25×10 -2 , with a maximum size ratio of 1.42%. The meshsize is 3060 × 1850 for ε = 0.5 and ε = 1, and 3300 × 1850 for ε = 2. The time step is set to δt = 2 × 10 -3 for Sk = 100 and δt = 5 × 10 -3 for Sk > 100. Regarding the boundary conditions at the domain ends, the normal velocity is set to zero to ensure a null ow rate far from the cylinders and the normal derivative of the tangential component is imposed to zero. The normal component of the pressure gradient is also set to zero, which is the usual boundary condition for the pressure eld when the ow rate is imposed.

Fluid forces

The time evolutions of the uid forces are depicted in gures 4, 5 and 6. The theoretical predictions show that the forces are sinusoidal functions whose amplitude and phase depend on Sk (viscous eects) and ε (connement eects).

To study this sensitivity, we plot in gure 7 the evolutions of the magnitude The solid lines refer to the least squares approximation and the dotted lines refer to the collocation approximation. The dimensionless separation distance is ε = 0.5 (black color), ε = 1 (red color), ε = 2 (blue color) and ε → ∞ (green color). On a) and b), the horizontal dashed lines are the asymptotic limits as Sk → ∞. Closed circles correpond to numerical predictions. On b) and c) the black, blue and red circles are indistinguishable because superimposed.

ε = 1 C 2 C 1 C 1 C 1 C 1 C 1 C 2 C 2 C 2 C 2 C 2 C
0.5 ε = 1 ε = 2 ε = ε → ∞ 0.5 ε = 1 ε = 2 ε = ε → ∞ -1 -0.8 -0.6 -0.4 -0.2 0
h j and the phase φ j . We observe that h j is maximum for the moving cylinder, diverges to innity when Sk → 0 and decreases to h 1 → π|m P OT cross | and

h 2 → π|m P OT self | as Sk → ∞ (inviscid uid).
The magnitude is also shown to be maximum for the small values of ε (strong connement) and to decrease to h 1 → 0 and h 2 → π|m ISO self | as ε → ∞ (isolated cylinders). Thus, as one would expect, the uid forces are all the more intense as both the viscous and connement eects are important. The gure 7 b) shows that the forces are in phase opposition, i.e. φ 1 = φ 2 + π, with φ 1 increasing from φ 1 → π/2 as Sk → 0 to φ 1 → π as Sk → ∞. We note that the connement has a very weak eect on the phase, leading to a slight increase of φ j with ε. The variations of φ j imply that the direction of the uid forces depends on Sk and, to a lesser extent on ε. From ( 36), the uid forces vanish and reverse their direction when {e it f jx } = {e i(t+φ j ) h j u jx } = h j sin (t + φ j ) = 0, i.e. t = -φ j + kπ, k ∈ Z. At that time, the dimensionless displacement u 2 = sin(t) of the moving cylinder equals u * 2 = ± sin (φ j ). In gure 7 c), we show that the uid forces cause the cylinders to attract (resp. repel) each other when

-1 ≤ u 2 < -|u * 2 | (resp. |u * 2 | ≤ u 2 < 1). In the nar- row range -|u * 2 | ≤ u 2 < |u * 2 |
, the cylinders are attracted (resp. repelled) to each other if the velocity of the moving cylinder is positive (resp. negative). An estimation of u * 2 is made possible from the observation that it is weakly sensitive to ε (at least for ε ≥ 0.5) and thus can be approximated by its limit as ε → ∞. From [START_REF] Basset | A Treatise on Hydrodynamics[END_REF]

and u * 2 = ± sin(φ 2 ) = ± sin(arg(f 2x /u 2x )), it comes that u * 2 ≈= ± sin arctan c ISO self m ISO self ≈ ± c ISO self m ISO self 2 + c ISO self 2 , (37)
which is the equation of the green line (ε → ∞) shown in gure 7 c). An asymptotic expansion of the modied Bessel functions K j entering in the definition of m ISO self and c ISO self , see [START_REF] Basset | A Treatise on Hydrodynamics[END_REF], yields that u 2 * = O Sk -1/2 as Sk → ∞. Finally, we note that the theoretical predictions for h j and φ j are successfully corroborated by the numerical simulations, in the sense that similar trends are clearly recovered. Still, we note that the numerical simulations are poorly sensitive to ε and slightly understimate the magnitude h 2 of the uid force acting on the moving cylinder, especially in the range of low Stokes numbers. A detailed discussion on the dierences between the theoretical and numerical approaches is reported in section 5.5.

Fluid added coecients

We now proceed with analyzing the evolutions of the uid added coecients m self , c self , m cross and c cross entering in the computation of the uid forces. The solid lines refer to the least squares approximation and the dotted lines refer to the collocation approximation. The dimensionless separation distance is ε = 0.5 (black color), ε = 1 (red color), ε = 2 (blue color) and ε → ∞ (green color). The horizontal dashed lines are the asymptotic limits [START_REF] Sarpkaya | Forces on cylinders and spheres in a sinusoidally oscillating uid[END_REF] as Sk → ∞. Closed circles correpond to numerical predictions.

1 -1/2 1 -1/2 1 1/2 1 1/2 ε → ∞ 0.5 ε = 1 ε = 2 ε = ε → ∞ 0.5 ε = 1 ε = 2 ε = ε → ∞ 0.5 ε = 1 ε = 2 ε = 0.5 ε = 1 ε = 2 ε = ε → ∞
The evolutions of m self and c self are depicted in gures 8 a) and b). We observe that m self and c self diverge to innity as Sk → 0 and decrease to m self → m P OT self and c self → 0 as Sk → ∞ (inviscid uid). The log-log plots in the insets of gures 8 a) and b) indicate that [START_REF] Gibert | Vibration of structures in a static uid medium[END_REF] In addition to the dependence on the Stokes number, m self and c self are also sensitive to the connement. The two coecients are maximum for the small values of ε (strong connement) and decrease to m self → m ISO self and

m self = m P ot self + O Sk -1/2 and c self = O Sk -1/2 as Sk → ∞.
c self → c ISO self as ε → ∞ (isolated cylinders).
As both Sk and ε tend to inn-ity, we recover the classical results for an isolated cylinder in a perfect uid, m self → 1 and c self → 0.

The evolutions of m cross and c cross are depicted in gures 8 c) and d). We observe that m cross is negative and converges to m cross → 0 as Sk → 0. As Sk increases, m cross rst decreases, then hits a minimum, and nally increases to

m cross → m P OT cross as Sk → ∞.
We hypothesize that the non-monotic variations of m cross are related to an antagonist competition between the viscous and the connement eects. The term c cross is also negative, diverges to c cross → -∞ as Sk → 0 and increases to c cross → 0 as Sk → ∞. The log-log plots in the insets of gures 8 c) and d) indicate that [START_REF] Landweber | Added masses and forces on two bodies approaching central impact in an inviscid uid[END_REF] The coecients m cross and c cross are also sensitive to the connement: they are minimum for the small values of ε (strong connement) and increase to m cross → 0 and c cross → 0 as ε → ∞ (isolated cylinders). In such a case, and as expected, there is no uid force acting on the stationary cylinder.

m cross = m P ot cross + O Sk -1/2 and c cross = O Sk -1/2 as Sk → ∞.
Here again, the theoretical predictions for the uid added coecients are successfully corroborated by the numerical simulations, in the sense that similar variations are recovered. However, we note that both approaches do not exactly exhibit the same sensitivity to the connement eect, leading to some deviations in the predictions, in particular concerning the self added coecients at low Stokes numbers. We discuss the possible origins of these deviations in the next section.

Discussion on numerics versus theory

The gure 8 shows that the simulations tend to underestimate m self and c self , and surestimate m cross and c cross . To quantify this deviation, we introduce the quantity ι, dened as the relative distance between the numerical and the theoretical predictions of some quantity Q : ι = |Q num. -Q th. |/|Q num. |. The gure 9 and the tables in appendix C show that ι is maximum for the small values of Sk and ε. We attribute this deviation to the fact that the theoretical approach is based on an approximation (least squares or collocation method) which loses its accurary when Sk and ε become small, as shown in the study of the residuals in gure 3 b). Also, the numerical simulation, which is based on a penalization method, hardly makes the dierence between the solid and the uid domains for the low values of Sk. In any case, the relative deviation for m self (resp. m cross ) is always smaller than ι ≤ 10% (resp. ι < 20%). The deviation for the damping coecients c self and c cross is more pronounced, with ι ≤ 50% and ι ≤ 35%, respectively. Note that the maximum deviations are observed for Sk ∈ [0, 400] , ε < 1, and are less important when using the least squares method. Even if the approximations of the theoretical and numerical approaches can be invoked, the slope steepness of the damping coecients also contributes to the enhancement of the relative deviation in such a range of Sk and ε. It follows that both approaches yield similar trends, bringing out the same behavior of the uid coecients, despite some deviations in the particular case of a very viscous uid (low Sk) in a conned environnement (ε < 1).

Conclusion

We have considered the problem of the small oscillations of two cylinders immersed in a viscous uid initially at rest. A theoretical approach based on an Helmholtz decomposition of the uid velocity vector and a bipolar coordinate system has been carried out to estimate the uid forces acting on the two cylinders. In addition to this new theoretical work, we also have developed a numerical approach based on a pseudo-penalization method. Such a numerical method has been shown particularly ecient in solving uid-structure interaction problems, in particular for moderate or high Stokes numbers.

We studied the case in which one cylinder is stationnary while the other one is imposed an harmonic motion. We show that the amplitude, the phase and the direction of the uid forces are sensitive to the Stokes number and the separation distance between the cylinders. The two forces are in phase opposition and their amplitude decreases to the inviscid limits as Sk increases. The eect of viscosity is to add to the ideal uid added coecients a correction term which scales as Sk -1/2 . When the separation distance increases, the uid coecients converge to the limits of an isolated cylinder derived by Stokes [START_REF] Stokes | On the eect of the internal friction of uids on pendulums[END_REF]. The theoretical predictions are successfully corroborated by the numerical simulations, in the sense that similar trends are recovered, despite some deviations for low Sk and ε.

As an improvement to our previous work on ideal uids [START_REF] Lagrange | A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant uid[END_REF], the new theoretical approach carried out in the present article is able to capture the eects of viscosity on the uid forces. It oers a simple and exible alternative to the fastidious and hardly tractable approach developed by [START_REF] Chen | Vibration of nuclear fuel bundles[END_REF]. To our knowledge, this is also the rst time that the pseudo-penalization method is presented in the context of relatively small Stokes numbers. As such, the present work should foster further developements of this easy to implement numerical method, to tackle complex uid-structure interaction problems.

A Evolutions of m (j)P OT self and m P OT cross

In this appendix, we study the variations of the uid added coecients m (j)P OT self and m P OT cross , given by ( 22). We show in gure A.1 that m (j)P OT self (resp. m P OT cross ) increases (resp. decreases) with the radius ratio r while it decreases (resp. increases) with the dimensionless separation distance ε. When r → 0, the cylinder C 1 transforms to a point and the system is equivalent to an isolated cylinder C 2 , leading to the classical result m (2)P OT self → 1. On the other hand, when r → ∞, the cylinder C 1 transforms to an innite plane and the system is equivalent to a cylinder C 2 near a wall. In such a case, we obtain

m (1)P OT self → ∞, (A.1a) m (2)P OT self → m W ALL self = -4 ∞ n=1 nε (2 + ε) 4 n + 16 n 2ε + 2 ε (2 + ε) + 2 -2 n -2ε + 2 ε (2 + ε) + 2 2n + 4 n , (A.1b) m P OT cross → m W ALL cross = 8 ∞ n=1 4 n nε (2 + ε) 4 n -2 ε (2 + ε) + 2ε + 2 2n .
(A.1c) Values of m W ALL self are presented in gure A.2, showing a perfect agreement with the predictions of [START_REF] Mazur | Motion of a circular cylinder near a vertical wall[END_REF] and [START_REF] Chen | Flow-Induced Vibration of Circular Cylindrical Structures[END_REF].

B Functions ϕ (j)

n and A (j) n

The functions ϕ (j) n and A (j) n appearing in [START_REF] Sarpkaya | Separated ow about lifting bodies and impulsive ow about cylinders[END_REF] are determined from the boundary conditions (16c), (16d). It yields a linear system of equations, whose solution is

              ϕ (1) n (ξ 1 , ξ 2 , l) ϕ (2) n (ξ 1 , ξ 2 , l) A (1) n (ξ 1 , ξ 2 , l) A (2) n (ξ 1 , ξ 2 , l)               = [M n (ξ 1 , ξ 2 , l)] -1               -2 nae -n|ξ 1 | sgn (ξ 1 ) -2 nae -n|ξ 1 | 0 0               , (B.1) with [M n (ξ 1 , ξ 2 , l)] =              
-n cosh (nξ 1 ) -n sinh (nξ 1 ) sinh (lξ 1 ) l cosh (lξ 1 ) l n sinh (nξ 1 ) n cosh (nξ 1 ) -n cosh (lξ 1 ) -n sinh (lξ 1 )

-n cosh (nξ 2 ) -n sinh (nξ 2 ) sinh (lξ 2 ) l cosh (lξ 2 ) l n sinh (nξ 2 ) n cosh (nξ 2 ) -n cosh (lξ 2 ) -n sinh (lξ 2 )               . (B.2)

C Tables of comparison numerics versus theory

In this appendix, we report the theoretical and numerical values of the uid added coecients m self , c self , m cross and c cross , for ε = 0.5 (

2Figure 1 .

 1 Figure 1. Schematic diagram of the system: two oscillating cylinders C j with radii R j , centers O j , displacement vectors U j (T ), are immersed in a uid of kinematic viscosity ν. The small oscillations of C j generate an incompressible uid ow. The midpoint of O 1 and O 2 is O and the separation distance is E.

Figure 2 .

 2 Figure 2. Sketch of the conformal mapping ζ = h(z) dened by equation (11). The conformal function ζ = h(z) maps C 1 and C 2 into two parallel lines of equation ζ = iτ 1 and ζ = iτ 2 .

Figure 3 .

 3 Figure 3. Evolutions of k LS , k COL , χ LS , χ COL , versus the Stokes number Sk. The solid lines refer to the least squares approximation and the dotted lines refer to the collocation approximation. The dimensionless separation distance is ε = 0.5 (black color), ε = 1 (red color) and ε = 2 (blue color). The radius ratio is r = 1.

2 C 2 C 2 C 2 C 2 C 1 C 1 C 1 C 1 C 1 CFigure 4 .

 22222111114 Figure 4. Dimensionless uid force {e it f jx } as a function of the dimensionless time t, for various Stokes numbers Sk. The dimensionless separation distance is ε = 0.5.

1

 1

Figure 5 .

 5 Figure 5. Dimensionless uid force {e it f jx } as a function of the dimensionless time t, for various Stokes numbers Sk. The dimensionless separation distance is ε = 1.

2 C 2 C 2 C 2 C 2 C 1 C 1 C 1 C 1 C 1 CFigure 6 .

 22222111116 Figure 6. Dimensionless uid force {e it f jx } as a function of the dimensionless time t, for various Stokes numbers Sk. The dimensionless separation distance is ε = 2.

Figure 7 .

 7 Figure 7. Evolutions of a) the magnitude h j = |f jx /u 2x | and b) the phase φ j = angle (f jx /u 2x ) versus the Stokes number Sk. On c), the arrows show the direction of the uid forces depending on Sk, ε and the dimensionless displacement of C 2 . The solid lines refer to the least squares approximation and the dotted lines refer to the collocation approximation. The dimensionless separation distance is ε = 0.5 (black color), ε = 1 (red color), ε = 2 (blue color) and ε → ∞ (green color). On a) and b), the horizontal dashed lines are the asymptotic limits as Sk → ∞. Closed circles correpond to numerical predictions. On b) and c) the black, blue and red circles are indistinguishable because superimposed.

Figure 8 .

 8 Figure 8. Evolutions of the uid added coecients versus the Stokes number, Sk.The solid lines refer to the least squares approximation and the dotted lines refer to the collocation approximation. The dimensionless separation distance is ε = 0.5 (black color), ε = 1 (red color), ε = 2 (blue color) and ε → ∞ (green color). The horizontal dashed lines are the asymptotic limits[START_REF] Sarpkaya | Forces on cylinders and spheres in a sinusoidally oscillating uid[END_REF] as Sk → ∞. Closed circles correpond to numerical predictions.

Figure 9 .

 9 Figure 9. Evolutions of the relative deviation, ι, versus the Stokes number, Sk. The solid lines refer to the least squares approximation while the dotted lines refer to the collocation approximation. The dimensionless separation distance is ε = 0.5 (black color), ε = 1 (red color), ε = 2 (blue color).

Figure A. 1 .

 1 Figure A.1. Fluid added coecients, given by (22), versus the radius ratio r. Evolution of a) m (1)P OT self , b) m (2)P OT self and c) m P OT cross . The horizontal dotted lines on b) and c) show the limits as r → ∞, see eq. (A.1). The dimensionless separation distance is ε = 0.5 (black color), ε = 1 (red color), ε = 2 (blue color).

ε 2 C 1 CFigure A. 2 .

 212 Figure A.2. Added mass coecient m W ALL self , given by (A.1) for a cylinder vibrating near a wall.

  0 10 -3 (regular over cylinder areas) -hmax= 3.25 10 -2 Time step : t= 2.0 10 -3 (SK=100) , t= 5.0 10 -3 (SK=300), t= 1.0 10 -2 (SK>300)

  0 10 -3 (regular over cylinder areas) -hmax= 3.25 10 -2 Time step : t= 2.0 10 -3 (SK=100) , t= 5.0 10 -3 (SK=300), t= 1.0 10 -2 (SK>300)

  The evolutions of χ LS and χ COL indicate that the theory becomes less accurate as the Stokes number and the dimensionless separation distance decrease (i.e. as the viscous and the connement eects becomes preponderant).

	4 Numerical simulation
	4.1 Solving the Navier-Stokes equations

Table C .

 C 1 Table of the uid added coecients and the relative deviation, ι. The notations LS and COL refer to the Least Squares and Collocation methods. The dimensionless separation distance is ε = 0.5.

  table C.1), ε = 1 (table C.2) and ε = 2 (table C.3). The numerical values correspond to the closed symbols shown in gure 8. The relative deviation ι is also reported in the tables.

	Comparison theory / numerics
	Sk		self m				c	self			cross m			cross c
		LS	COL	LS	COL	LS	COL	LS	COL
			Numerics		Numerics		Numerics		Numerics
			LS		COL		LS		COL		LS		COL		LS		COL
	100	1.35	1.36	0.390	0.413	-0.293	-0.296 -0.0910 -0.0990
			1.30				0.326			-0.283			-0.0754
		3.8%	4.6%	20%	27%	3.5%	4.6%	21%	31%
	300	1.22	1.22	0.210	0.219	-0.266	-0.268 -0.0475 -0.0508
			1.19				0.189			-0.240			-0.0419
		2.5%	2.5%	11%	16%	11%	12%	13%	21%
	500	1.18	1.18	0.159	0.165	-0.258	-0.259 -0.0357 -0.0379
			1.14				0.156			-0.231			-0.0348
		3.5%	3.5%	1.9%	5.8%	12%	12%	2.6%	8.9%
	700	1.15	1.16	0.133	0.138	-0.253	-0.254 -0.0297 -0.0314
			1.13				0.135			-0.228			-0.0300
		1.8%	2.6%	1.5%	2.2%	11%	11%	1.0%	4.7%
	900	1.14	1.14	0.116	0.120	-0.250	-0.251 -0.0258 -0.0273
			1.12				0.122			-0.226			-0.0270
		1.8%	1.8%	4.9%	1.6%	11%	11%	4.4%	1.1%
	Relative deviation :
	Theory Numerics
	Numerics

Table C .

 C 2 Table of the uid added coecients and the relative deviation, ι. The notations LS and COL refer to the Least Squares and Collocation methods. The dimensionless separation distance is ε = 1.

	Comparison theory / numerics
	Eps= 2.0
	Sk		self m				c	self			cross m				cross c
		LS	COL	LS	COL	LS	COL	LS	COL
			Numerics		Numerics		Numerics		Numerics
			LS		COL		LS		COL		LS		COL		LS		COL
	100	1.31	1.32	0.329	0.338	-0.160	-0.162 -0.0465 -0.0498
			1.27			0.310			-0.135			-0.0372
		3.1%	3.9%	6.1%	9.0%	19%	20%	25%	34%
	300	1.19	1.19	0.185	0.188	-0.146	-0.147 -0.0242 -0.0254
			1.16			0.181			-0.123			-0.0208
		2.6%	2.6%	2.2%	3.9%	19%	20%	16%	22%
	500	1.15	1.15	0.142	0.144	-0.142	-0.142 -0.0182 -0.0189
			1.12			0.149			-0.119			-0.0174
		2.7%	2.7%	4.7%	3.4%	19%	19%	4.6%	8.6%
	700	1.13	1.13	0.120	0.121	-0.139	-0.140 -0.0152 -0.0157
			1.11			0.130			-0.117			-0.0150
		1.8%	1.8%	7.7%	6.9%	19%	20%	1.3%	4.7%
	900	1.11	1.11	0.105	0.106	-0.138	-0.138 -0.0132 -0.0136
			1.10			0.117			-0.116			-0.0135
		0.91%	0.91%	10%	9.4%	19%	19%	2.2%	0.74%
	Relative deviation :
	Theory Numerics
	Numerics

Table C .

 C 3 Table of the uid added coecients and the relative deviation, ι. The notations LS and COL refer to the Least Squares and Collocation methods. The dimensionless separation distance is ε = 2.

Email addresses: romain.lagrange@cea.fr, romain.g.lagrange@gmail.com (Romain Lagrange).


R Lagrange

D Averbuch
email: daniel.averbuch@ifpen.fr

Solution methods for the growth of a repeating imperfection in the line of a strut on a nonlinear foundation

Keywords: Buckling, Nonlinear elastic foundation, Imperfection, Stability and bifurcation

This paper is a theoretical and numerical study of the uniform growth of a repeating sinusoidal imperfection in the line of a strut on a nonlinear elastic Winkler type foundation. The imperfection is introduced by considering an initially deformed shape which is a sine function with an half wavelength. The restoring force is either a bi-linear or an exponential prole. Periodic solutions of the equilibrium problem are found using three dierent approaches: a semi-analytical method, an explicit solution of a Galerkin method and a direct numerical resolution. These methods are found in very good agreement and show the existence of a maximum imperfection size which leads to a limit point in the equilibrium curve of the system. The existence of this limit point is very important since it governs the appearance of localization phenomena.

Using the Galerkin method, we then establish an exact formula for this maximum imperfection size and we show that it does not depend on the choice of the restoring force. We also show that this method provides a better estimate with respect to previous publications. The decrease of the maximum compressive force supported by the beam as a function of the imperfection magnitude is also determined. We show that the leading term

Introduction

The subject of beam buckling can be found in several situations in industrial applications. Among the most studied are the thermal track buckling and the buckling of subsea pipelines under the eect of temperature and/or pressure. In the latter case, buckling can appear and in both the vertical and horizontal planes, according to the existing restraints imposed by the environment and the backll. Numerous authors have studied these two applications due to their high practical importance, and have proposed solution methods to determine buckling loads and post-buckling situations. Along time, the techniques used have progressed, based rstly on analytical analyses and latter on numerical methods mostly derived from nite elements models. Thus, based on some early work by [START_REF] Kerr | On the stability of the railroad track in the vertical plane[END_REF][START_REF] Kerr | Analysis of thermal track buckling in the lateral plane[END_REF] studying the stability of railway tracks subjected to thermal buckling, several authors such as [START_REF] Bournazel | Vertical buckling of buried pipes[END_REF]; [START_REF] Hobbs | Pipeline buckling caused by axial loads[END_REF][START_REF] Hobbs | In-service buckling of heated pipelines[END_REF] have proposed solution methods where the equilibrium equations were solved in post-buckling congurations to establish relevant buckling loads. In these works, the soil was supposed rigid, while the external forces acting on the beam was assumed as constant as a dead weight or constant friction force. One of the key features of these theories is the fact that the loss of contact (or movement) induces a loss of global stiness of the structure which leads to subcritical buckling and innite buckling loads if no imperfection is assumed. Using slightly dierent arguments, other models were proposed in [START_REF] Croll | A simplied model of upheaval thermal buckling of subsea pipelines[END_REF] and Maltby and Calladine (1995b). In the latter work, equilibrium equations were obtained by assuming sine deections in the post-buckling situations and using an approximate Galerkin solution method. This method was compared with numerical solutions in Maltby and Calladine (1995a) and against experimental results in Maltby and Calladine (1995b). Though using an approximate solution method, the approach showed good results with respect to the numerical simulation. In order to improve the earlier methods, numerical models were developed in [START_REF] Ju | Thermal buckling of oshore pipelines[END_REF]; [START_REF] Klever | A dedicated nite-element model for analyzing buckling response of submarine pipelines[END_REF]; [START_REF] Leroy | Behavior of buried exible pipelines[END_REF]; [START_REF] Yun | Model for beam-mode buckling of buried pipelines[END_REF] to incorporate for instance additional non linear effects in the models, such as non linear geometric and material models. One of the key aspects of work related to the study of the upheaval buckling is the study of the localization phenomenon, which was suspected early in the pioneering work of Tveergard andNeedleman (1980, 1981). This aspect was analyzed through numerical simulations in these rst papers and in [START_REF] Hunt | Comparative lagrangian formulations for localized buckling[END_REF]; [START_REF] Hunt | Principles of localized buckling for a strut on an elastoplastic foundation[END_REF], and through analytical approaches based on a double scale expansion of the equilibrium equations in [START_REF] Potier-Ferry | Foundations of elastic postbuckling theory[END_REF], for a beam resting on an elastic non-linear foundation. Based on these results, the study of the localization phenomenon, was continued by using Galerkin techniques (see [START_REF] Wadee | Eects of periodic and localized imperfections on struts on nonlinear foundations and compression sandwich panels[END_REF][START_REF] Whiting | A galerkin procedure for localized buckling of a strut on a nonlinear elastic foundation[END_REF] using the displacement envelopes obtained through the double-scale expansion of [START_REF] Potier-Ferry | Foundations of elastic postbuckling theory[END_REF]. A lot of attention has also been paid to the estimation of the mechanical restraint induced by the soil friction and to the eect of backll on the pipeline behavior, since the corresponding forces were found to highly inuence the mechanical behavior of the pipe. In order to feed the corresponding models, experiments were performed (see [START_REF] Palmer | Uplift resistance of buried submarine pipelines: comparison between centrifuge modeling and full-scale tests[END_REF][START_REF] Schaminee | Soil response for pipeline upheaval buckling analyses: full-scale laboratory tests and modelling[END_REF]Trautmann et al., 1985a,b) either through centrifuge testing of small-scale pipeline models or through direct testing of buried full scale pipe sections.

The present paper is an attempt to provide additional solutions for the study of the growth of a repeating imperfection in the line of a strut on a nonlinear foundation. In this work, the foundation is supposed to act through an either bi-linear or exponential regularized friction model relating the interaction line force to the transverse displacement (see section 2). These two kinds of models are indeed found in the above mentioned papers describing the soilpipe interaction models. In the former case, a solution method (piecewise solution) in section 3 is proposed by explicitly solving the equilibrium equation in the regions where the foundation acts linearly and where the friction force is constant, and by connecting the two solutions by adequate boundary solutions. Alternatively, a Galerkin approach of the same problem is developed in section 3 and leads to an explicit solution of the problem, which is developed for the two regularization models. The piecewise solution and the Galerkin approach are consequently compared together and with numerical solutions of the problem in section 4. The post-buckling problem is then studied through the Galerkin approach which provides precise analytical solutions, focusing on the characteristics of a limit point (see section 5) in the equilibrium curve depending on the magnitude of the initial imperfection. The imperfection is W 0 , the compressive load is P and the buckling displacement is W . Springs have a nonlinear elastic force-displacement relationship P (W ).

Formulation of the dierential equation

This section formulates the dierential equation for the growth of an imperfection in the line of a strut on a nonlinear Winkler-type foundation, see gure 1. The imperfection is introduced by supposing an initially deformed shape W 0 whose form is

W 0 = A 0 sin ( π L X ) , ( 1
)
with A 0 the amplitude of the imperfection, L its length and X the longitudinal coordinate. The compressive load is P and the restoring force per unit length is P . These two forces are assumed to be conservative. The dierential equation governing the deection W may be derived either: directly by equilibration of forces; or from the Principal of Virtual Work; or using an energy formulation. The latter approach is adopted here; the total potential energy at rst order being

V = L ∫ 0   1 2 EIW ′′ 2 -P ( 1 2 W ′ 2 + W 0 ′ W ′ ) - W ∫ 0 P (t) dt   dX, (2)
where a prime indicates dierentiation with respect to X. The rst term is the strain energy of bending (EI is the bending stiness of the strut), the second is the work done by the load P and the remainder is the energy stored in the elastic foundation. Equilibrium is given by stationary values of V . In what follows the strut is assumed to be simply supported, such that the conditions at the boundaries are W (0) = W (L) = 0. The calculus of variations on (2) gives, for a virtual displacement δW such that δW (0) = δW (L) = 0

δV = EI [ δW ′ W ′′ ] 0 L + L ∫ 0 [ EIW ′′′′ + P ( W ′′ + W ′′ 0 ) -P (W ) ] δW dX, (3)
so that the Euler-Lagrange equation and the conditions at the boundaries for a simply supported strut are

EIW ′′′′ + P W ′′ -P (W ) = -P W ′′ 0 , ( 4a
)
W (0) = 0, (4b)
W ′′ (0) = 0, (4c)
W (L) = 0, (4d)
W ′′ (L) = 0. (4e)
In equation (4a), P is the compressive load before buckling. The compressive load after buckling, considering a strut of section S, should be written as N = P -ES 2L L ∫ 0 (W ,X ) 2 dX, last term of this expression being a geometric shortening which allows for the additional length introduced by the lateral movement. Therefore, N should be used in the equation for equilibrium. However, [START_REF] Tveergard | On localized thermal track buckling[END_REF] have shown that the buckle will only become unstable if N (y) has a maximum is correct for an isolated halfwave but is not correct for a long strut which contains a sequence of half-waves end-to-end. In such a case the key point is that a localization of the buckling, in which one particular half-wave grows at the expenses of its neighbours, can occur whenever the curve P (y) has a maximum. Under this consideration, and as Maltby and Calladine (1995b) did, we use P instead of N as the load parameter.

The restoring force

The restoring force per unit length is assumed to be nonlinear and two particular P functions are considered. The rst one is referred as the bi-linear function and is dened by

P (W ) =    -KW if |W | < ∆, -K∆ if W > ∆, K∆ if W < -∆, ( 5
)
where K is the linear stiness and ∆ the mobilization. These two constants are positive. The second P function considered in this paper is referred as the exponential prole and is dened by

P (W ) =    -K∆ ( 1 -e -W ∆ ) if W > 0, K∆ ( 1 -e W ∆ ) if W < 0. (6)
The two P functions share the same initial slope and limiting force (see gure 2(a)).

Nondimensionalization

Let's introduce a characteristic length

L char = ( EI K )1 4 and nondimensional quantities l = L L char , x = X L char , w = W ∆ , w 0 = W 0 ∆ , λ = P KL char 2 , p = P K∆ . ( 7
)
Hence, from (4) the deection w is solution of the dierential problem

w ′′′′ + λw ′′ -p (w) = -λw ′′ 0 , ( 8a
)
w (0) = 0, (8b)
w ′′ (0) = 0, (8c)

w (l) = 0, (8d)
w ′′ (l) = 0, (8e)
the dimensionless imperfection being

w 0 = a 0 sin ( π l x ) , (9)
with a 0 = A 0 ∆ . After nondimensionalization, the bi-linear function rewrites

p (w) =    -w if |w| < 1, -1 if w > 1, 1 if w < -1, ( 10
)
such that the dimensionless mobilization and limiting force equal 1. The nondimensional exponential prole is given by

p (w) = { -(1 -e -w ) if w > 0, 1 -e w if w < 0. (11)
Figure 2(b) shows the evolution of these two dimensionless restoring forces for w > 0.

Theoretical resolution

In this paper, imperfections with dimensionless lengths l < √ 2π are considered. In order to give a practical meaning to this inequation, let's consider the case of an imperfect railway. The most common rails in France (rail type 50E6) have a weight per meter of about m = 500 N m -1 , a length from 10 to 400 m (length between two joints connection) and a bending stiness of about EI = 4 × 10 6 N m 2 . Considering a coecient of friction φ between steel and concrete of 0.4 and a mobilization ∆ from 0.1 to 1 mm it comes a linear stiness K = mgφ ∆ (g being the gravity) from 2 × 10 5 to 2 × 10 6 N m -2 and a characteristic length L char from 1 to 2 m. Therefore l < √ 2π would correspond to a repeating sinusoidal imperfection in the railway whose half-wavelength is no more than 10 m. In a more general way, l < √ 2π deals with imperfections whose length does not exceed some meters. Therefore the specic calculations described in this paper have been made in the context of a small-scale experimental setup.

For l < √ 2π, the rst buckling mode predicted by the linear analysis is excited when

λ = λ c = ( π l ) 2 + ( π l ) -2 , ( 12
)
and it has the same shape as the imperfection. In what follows we will also introduce the Euler load λ e

λ e = ( π l ) 2 , ( 13
)
which is the buckling load of the rst mode when the restoring force equals 0.

Two theories are developed to solve the equilibrium problem. The rst one, named piecewise solution theory is an exact resolution of the equilibrium problem when the bi-linear restoring force is considered. The second theory is based on a Galerkin method: it leads to an approximate resolution of the equilibrium problem by considering equally the bi-linear or the exponential restoring force. To initiate this method, the deection shape is assumed to be a sinusoid, as the imperfection. Explicit solutions of the Galerkin equation are obtained without any assumptions.

Piecewise solution theory

The principle of the piecewise solution theory is to solve (8a) on each piece of the bi-linear function. Then, the solutions are connected thanks to the boundary conditions and assuming the continuity of w, w ′ , w ′′ and w ′′′ at two connecting points x 1 and x 2 . Substituting p (w) by -w and by -1 in (8a) yields respectively

w ′′′′ 1 + λw ′′ 1 + w 1 = -λw ′′ 0 , ( 14a
)
w ′′′′ 2 + λw ′′ 2 = -1 -λw ′′ 0 . ( 14b
)
The solutions w 1 and w 2 belong to two ane spaces of dimension 4 given by

w 1 = w 1h + w 1part , ( 15a
)
w 2 = w 2h + w 2part , ( 15b
)
where w 1h and w 2h satisfy the homogeneous equations

w 1h ′′′′ + λw 1h ′′ + w 1h = 0, (16a)
w 2h ′′′′ + λw 2h ′′ = 0. ( 16b
)
Inserting w 1h = Ae r 1 x and w 2h = Be r 2 x in ( 16) yields two algebraic equations

r 4 1 + λr 1 2 + 1 = 0, (17a)
r 4 2 + λr 2 2 = 0, (17b)
whose solutions are

r 1 (1) = 1 2 [ -2λ + 2 ( λ 2 -4 ) 1/2 ] 1/2 , (18a) r 1 (2) = - 1 2 [ -2λ + 2 ( λ 2 -4 ) 1/2 ] 1/2 , (18b) r 1 (3) = 1 2 [ -2λ -2 ( λ 2 -4 ) 1/2 ] 1/2 , (18c) r 1 (4) = - 1 2 [ -2λ -2 ( λ 2 -4 ) 1/2 ] 1/2 , ( 18d
)
and

r 2 = 0, ± √ λi. ( 19
)
Let's introduce α i and ω i the real and imaginary parts of r 1 (i) . For λ < 2, the roots r 1 are complex numbers and the function w 1h is

w 1h = A 1 e α 1 x cos (ω 1 x) + A 2 e α 2 x cos (ω 2 x) + A 3 e α 1 x sin (ω 1 x) +A 4 e α 2 x sin (ω 2 x) , ( 20
)
with A 1 , A 2 , A 3 and A 4 four real constants. For λ = 2, the roots r 1 are double imaginary numbers (r

1 (1) = r 1 (3) and r 1 (2) = r 1 (4) ). The function w 1h is w 1h = (A 1 x + A 2 ) cos (x) + (A 3 x + A 4 ) sin (x) . ( 21
)
For λ > 2, the roots r 1 are imaginary numbers. The solution w 1h is

w 1h = A 1 cos (ω 1 x) + A 2 cos (ω 3 x) + A 3 sin (ω 1 x) + A 4 sin (ω 3 x) . ( 22
)
For λ ̸ = 0 the function w 2h is given by

w 2h = B 1 + B 2 x + B 3 cos ( √ λx ) + B 4 sin ( √ λx ) , ( 23
)
with B 1 , B 2 , B 3 and B 4 four real constants.

Since w 1h and w 2h depend on 4 undetermined constants, these functions will be noted

w 1h (x, A 1 , A 2 , A 3 , A 4 ) and w 2h (x, B 1 , B 2 , B 3 , B 4 ).
The functions w 1part and w 2part appearing in (15a) and (15b) are particular solutions of (14a) and (14b), with w 0 given by (9). Searching w 1part and w 2part as w 1part = Aw 0 and w 2part = Cx 2 + Aw 0 , yields

w 1part = λ λ c -λ w 0 , ( 24a
)
w 2part = - 1 λ x 2 2 + λ λ e -λ w 0 . (24b)
Finally, the functions w 1 and w 2 , dened by ( 15a) and (15b), are

w 1 = w 1h (x, A 1 , A 2 , A 3 , A 4 ) + λ λ c -λ w 0 , ( 25a
)
w 2 = w 2h (x, B 1 , B 2 , B 3 , B 4 ) - 1 λ x 2 2 + λ λ e -λ w 0 . ( 25b
)
In what follows, we will introduce the function w 3 dened as

w 3 = w 1h (x, C 1 , C 2 , C 3 , C 4 ) + λ λ c -λ w 0 . ( 26
)
This function is a solution of (14a), as w 1 , but with 4 dierent undetermined constants C 1 , C 2 , C 3 and C 4 .

The aim of the piecewise solution theory is to search for a deection w solution of ( 8) by connecting the functions w 1 , w 2 and w 3 at two unknown points x 1 and x 2 such that

w =            w 1 if x ∈ [0, x 1 [ , 1 if x = x 1 , w 2 if x ∈ ]x 1 , x 2 [ , 1 if x = x 2 , w 3 if x ∈ ]x 2 , l] . ( 27
)
With these notations, the boundary conditions (8b), (8c), ( 8d) and ( 8e) rewrites

w 1 (0) = 0, (28a)
w ′′ 1 (0) = 0, (28b)
w 3 (l) = 0, ( 28c
)
w ′′ 3 (l) = 0. ( 28d
)
The continuity of the displacement, the tangent, the curvature and the shear at x 1 yields

w 1 (x 1 ) -w 2 (x 1 ) = 0, (29a)
w ′ 1 (x 1 ) -w ′ 2 (x 1 ) = 0, (29b)
w ′′ 1 (x 1 ) -w ′′ 2 (x 1 ) = 0, (29c)
w ′′′ 1 (x 1 ) -w ′′′ 2 (x 1 ) = 0, (29d)
and at x 2

w 2 (x 2 ) -w 3 (x 2 ) = 0, (30a) w ′ 2 (x 2 ) -w ′ 3 (x 2 ) = 0, ( 30b
)
w ′′ 2 (x 2 ) -w ′′ 3 (x 2 ) = 0, (30c) w ′′′ 2 (x 2 ) -w ′′′ 3 (x 2 ) = 0. (30d)
Equations ( 28), ( 29) and ( 30) lead to a linear system with 12 equations and 12 unknowns (i.e. the amplitudes A i , B i and C i ). A matrix representation of this system is

G (x 1 , x 2 ) a = b (w 0 , x 1 , x 2 ) , ( 31
)
with G a 12 by 12 real matrix and a the vector of the unknown amplitudes.

The vector b contains the particular solutions w 1part and w 2part which depend on w 0 , x 1 and

x 2 . If det (G) ̸ = 0 then a (w 0 , x 1 , x 2 ) = G -1 (x 1 , x 2 ) b (w 0 , x 1 , x 2 ) . ( 32
)
Equation ( 32) express the amplitudes as functions of x 1 and x 2 . These two connecting points are obtained by solving the nonlinear system

f 1 (w 0 , x 1 , x 2 ) = w 1 (x 1 ) -1 = 0, (33a)
f 2 (w 0 , x 1 , x 2 ) = w 2 (x 2 ) -1 = 0. (33b)
The numerical resolution of this system has been carried out with Matlab, using the fzero function. This function tries to nd a zero of (33) near X 0 , X 0 being a vector of length two. Depending on λ (33) has 0 or several solutions. For small λ, the restoring force p remains linear, such that w = w 1 for any x ∈ [0, l]. Thus, (33) has no solution. For high λ, the restoring force p is nonlinear: the existence and the uniqueness of a solution for (33) is not trivial. In this paper, we only search for a solution which satises

x 1 ∈ [0, l/2],
x 2 ∈ [l/2, l] and x 2 = l -x 1 : the connecting points are symmetric relative to x = l/2. This condition is specied adjusting the X 0 vector used by the fzero function. Typically, X 0 = [l/2; l/2] is a good candidate to easily nd the symmetric connecting points. Once the connecting points are calculated, we determine the amplitudes A i , B i and C i thanks to (32), the functions w 1 , w 2 and w 3 thanks to (25a), (25b), ( 26) and nally the deection w via (27).

Galerkin method

The Galerkin procedure (see [START_REF] Fox | An introduction to the Calculus of Variations[END_REF]) may be seen as being derived from (3) by assuming that the modes which go to make up w are given by

w = n ∑ i=1 y i ϕ i , ( 34
)
where each y i is an undetermined amplitude of each shape function ϕ i . Depending on the form of w in (34), we can perform periodic or localized buckling analysis. For a very long imperfection in the line of a strut on a cubic foundation, [START_REF] Whiting | A galerkin procedure for localized buckling of a strut on a nonlinear elastic foundation[END_REF] performed the latter, using the functions predicted by the asymptotical analysis (see [START_REF] Wadee | Asymptotic and rayleigh-ritz routes to localized buckling solutions in an elastic instability problem[END_REF] as test functions.

The amplitudes of each shape function are determined numerically using a variable-order variable-step Adams method. In this paper we do not search for localized solutions: we use a unique test function which has the same shape as the imperfection. It means that the deection w is searched as

w = y sin ( π l x ) , ( 35
)
with y > 0.

Inserting δw in the dimensionless form of (3) gives

l ∫ 0 sin ( π l x ) [ w ′′′′ + λ ( w ′′ + w ′′ 0 ) -p (w) ] dx = 0. (36)
Writing the restoring force as p (w) = -w -N (w) yields a relation between the load λ and the amplitude y

λ = 1 a 0 + y [ λ c y + Q (y) λ e ] , (37)
with

Q (y) = 2 l ∫ l 0 sin ( π l x ) N ( y sin ( π l x )) dx. (38)
The assumption N ( y sin

( π l x )) = N (y) sin ( π l x
) introduced by Maltby and Calladine (1995b) yields N = Q. Such an assumption will not be introduced in this paper. However, we will see how this assumption aects the result.

Note that since the integrand function in ( 38) is l/2-periodic, the Q function does not change when the deection is searched as w = y sin ( nπ l x

) . Then, when the imperfection is w 0 = a 0 sin ( nπ l x

) , the equilibrium paths are still given by (37) with λ e = ( nπ l ) 2 and λ c = λ e + λ -1 e .

In practice, the equilibrium paths predicted by the Galerkin method are plotted using (37), varying y and evaluating λ.

Bi-linear restoring force

Considering the bi-linear restoring force, the Q function rewrites

Q (y) = 2H (y -1) π { y [ arcsin ( 1 y ) - π 2 ] + ( y 2 -1 y 2 )1 2 } , ( 39
)
where H is the Heaviside function dened by H (y -1) = 0 if y < 1 and H (y -1) = 1 if y ≥ 1. The proof of this result is reported in Appendix A.

Exponential restoring force

Considering the exponential restoring force, the Q function rewrites

Q (y) = 2 [ I1 (y) -L1 (y)] -y, ( 40
)
where I 1 and L 1 are respectively the modied Bessel and Struve functions of parameter 1. The proof of this result is reported in Appendix B.

Theoretical and numerical results

In this section, the equilibrium paths predicted by the piecewise solution theory, the Galerkin procedure and a numerical resolution of (8) are compared. We also determine the inuence of the restoring force (bi-linear or exponential) on the shape of the equilibrium paths.

Piecewise solution theory

The equilibrium paths predicted by the piecewise solution theory are depicted in gure 3(a) for a 0 from 0 to 1.19. For a small imperfection size, the equilibrium path shows the load increasing at rst but then hits a maximum (limit point, or saddle-node bifurcation point) that is below λ c (or equals for a 0 = 0), and the rest of the path asymptotically decreases to the Euler load when max (w) → ∞. Thus, the system is subcritical. Moreover, the greater the imperfection size, the greater the reduction in the maximum load. Thus the system is imperfection sensitive. For high imperfection sizes, the equilibrium path increases monotically and λ → λ e < λ c when max (w) → ∞.

From these observations we infer the existence of a critical amplitude a 0c such that Function w 1 (solid line), function w 2 (dashed line), function w 3 (dash-dotted line). The functions w 1 and w 2 are connected at x = x 1 ≈ 0.618 (rst vertical dotted line). The functions w 2 and w 3 are connected at x = x 2 ≈ 2.382 (second vertical dotted line). The horizontal dotted line represents the mobilization. l = 3, a 0 = 0, λ = 1.75.

• if a 0 > a 0c then the equilibrium paths do not have a limit point and the equilibrium states are stable,

• if a 0 < a 0c then the equilibrium paths have a limit point (y m , λ m ).

For y < y m (resp. y > y m ) the equilibrium states are stable (resp. unstable). An unstable equilibrium state is represented in gure 4 by connecting the functions w 1 , w 2 , and w 3 .

The determination of a 0c is reported in section 5.

Galerkin method and numerical results

Bi-linear restoring force

For the bi-linear function, the equilibrium paths predicted by the Galerkin method and those obtained via a numerical resolution of (8), using the ODE45 solver from Matlab (this routine uses a variable step Runge-Kutta method), are represented in gures 3(b) and 3(c). The equilibrium paths are identical to those predicted by the piecewise solution theory (the relative error between the two theories and the numerical resolution being less than 0.1%). Since the Galerkin test function and the imperfection have the same shape, in the case of a bi-linear restoring force the deection is an amplication of the imperfection.

Exponential restoring force

For the exponential prole, the equilibrium paths predicted by the Galerkin method and those obtained via a numerical resolution of ( 8), are represented in gures 5 and 6. Once again they are identical, thus the deection w is an amplication of the imperfection. The equilibrium paths predicted by Maltby and Calladine (1995b) have also been reported in gure 5. The assumption introduced by Maltby and Calladine (1995b) (see section 3.2) leads to an underestimation of the load. In section 5 we will show that the existence of a limit point is also aected by this assumption.

In gure 5 the equilibrium paths predicted for a bi-linear restoring force have been reported in order to discuss the inuence of the restoring force. It appears that the restoring force (bi-linear or exponential) has no inuence on the shape of the equilibrium paths (the variations are preserved, it always exists a limit point for small imperfection sizes, the same asymptotes are recovered). Nevertheless, the equilibrium paths for an exponential prole are below the equilibrium paths for a bi-linear prole. Therefore, the choice of the restoring force has a non negligible inuence on the maximal load acceptable by the system. This point is the purpose of the next section.

Limit point

A limit point corresponds to a maximum of λ. Therefore, if it exists, this point (y m , λ m ) satises dλ dy (y m ) = 0, that is to say using (37)

(a 0 + y m ) dQ dy (y m ) -Q (y m ) = -a 0 λ c λ e . ( 41
)
This equation yields

Q(ym) λe = a 0 λ c -(a 0 + y m ) f (a 0 )
with f a function of a 0 . Inserting this last expression in (37) yields

λ m = λ (y m ) = λ c -f (a 0 ) . ( 42
)
In practice, the coordinates of a limit point are obtained by varying y m and looking for an amplitude a 0 which satises (41). that y m is an increasing function of a 0 and diverges for the same critical amplitude. This amplitude is determined by enforcing y m → ∞ in (41). It yields

a 0c = 4 π λ -2 e , ( 43
)
whose dimensional equivalent form is

A 0c = 4 π L 4 π 4 EI K∆. ( 44
)
This result is of great importance since it governs the appearance of localized solutions, as demonstrated in [START_REF] Tveergard | On localized thermal track buckling[END_REF]. In that paper, they indeed showed the close link between the existence of a limit point in the λ (y) curve and localization, the bifurcation into a localized form occurring a little way beyond this limit point.

Considering the bi-linear restoring force, when a 0 → 0, equation ( 41) has an innity of solutions y m in ]0, 1]. This observation is coherent with the equilibrium paths depicted in gures 3 since the points lying on the plateau λ = λ c are a maximum.

Considering the exponential restoring force, y m → 0 when a 0 → 0. This result is coherent with the equilibrium paths depicted in gures 5(a) and 6(a) since (0, λ c ) is a maximum.

Substituting Q (y m ) by N (y m ) = -y m + 1 -e -ym in (41) (assumption made by Maltby and Calladine (1995b)) and enforcing y m → ∞ yields

a 0c = λ e , ( 45
)
whose dimensional equivalent form is

A 0c = L 4 π 4 EI K∆, (46)
which is 4 π times smaller than the exact prediction (43). This result means that a stable equilibrium (i.e. an equilibrium path without a limit point) predicted with the assumption introduced by Maltby and Calladine (1995b), could in fact be unstable (see gure 9).

Figure 7(b) represents the evolution of the limit point λ m as a function of a 0 . Considering the two restoring forces, the same asymptotes are observed: λ m → λ e when a 0 → a 0c and λ m → λ c when a 0 → 0. We also observe that the assumption introduced by Maltby and Calladine (1995b) (see section 3.2) leads to an underestimation of the limit load. These observations are coherent with the equilibrium paths depicted in gures 3, 5 and 6.

Figure 8 shows the evolution of λ c -λ m as a function of a 0 , in a logarithmic scale. This picture allows to determine a scaling for the function f appearing in (42). Considering the bi-linear restoring force, for small amplitudes a 0 , λ m -λ c scales as

λ m -λ c ∼ -a 0 . ( 47
)
When the exponential restoring force is considered this scaling becomes

λ m -λ c ∼ -a 1/2 0 . ( 48
)
We observe in gure 8 that the assumption introduced by Maltby and Calladine (1995b) has no inuence on this last scaling.

Conclusion

In this paper, the growth of a repeating sinusoidal imperfection in the line of a strut on a nonlinear elastic Winkler type foundation is considered. The imperfection is introduced by considering an initially sinusoidal deformed shape with an half wavelength. The imperfection length is chosen such that the buckle mode predicted by the linear theory has the same shape as the imperfection (rst buckle mode). The nonlinearities are only due to the restoring force provided by the foundation. This restoring force is expressed as a force-displacement relationship which is either a bi-linear or an exponential function. The equilibrium problem is solved using three dierent methods. The rst one, named piecewise solution theory, is dedicated to the bi-linear prole and leads to an exact resolution of the equilibrium problem. The second one is available whatever the restoring force and is based on a Galerkin procedure. This procedure is initiated with a test function which has the same shape as the imperfection. It yields an explicit relation between the compressive load and the amplitude of the test function. This expression is an exact solution of the Galerkin equation and gives an approximate solution of the equilibrium problem. The last method is a numerical resolution of the equilibrium problem, using the ODE45 solver from Matlab. These three solving methods yield the same results: whatever the restoring force (bi-linear or exponential), the bifurcation is subcritical, the system is imperfection sensitive and the deformed shape is an amplication of the default. Moreover, it exists a critical imperfection size a 0c = 4 π λ -2 e (λ e being the Euler load) which does not depend on the restoring force and such that • if a 0 > a 0c , then the equilibrium path shows the load increasing monotically and remains asymptotic to the Euler load.

• if a 0 < a 0c , then the equilibrium path shows the load increasing at rst but then hits a limit point and the rest of the path is asymptotic to the Euler load.

This paper provides a better estimate of a 0c with respect to previous publications.

For each restoring fore, an approximate mathematical rule is derived relating the imperfection size a 0 to the corresponding limit load λ m . Considering the bi-linear prole (resp. the exponential prole) the limit point scales as λ m -λ c ∼ -a 0 (resp. λ m -λ c ∼ -a 1/2 0 ), where λ c is the critical load issued from the classical linear analysis. Therefore, the scaling of the limit point depends on the regularization method.

In this paper, the restoring force and the compressive load are independant. Nevertheless, in some industrial applications (such as in drilling problems) the restoring force slightly depends on the axial compressive load. Therefore, we are currently carrying out a study with a bi-linear restoring force proportional to the axial load. First results issued from the Galerkin approach indicate that it is necessary to redene the dimensionless parameters, leading to new scalings for the critical imperfection size and the limit load.

where Sgn is the sign function. Since p = -w -N (w) it comes N (w) = (Sgn (w) -w) H (|w| -1) .

(A.2)

For an imperfection with an half-wavelength, it can be assumed that w > 0.

Then, the Q function is (see (38))

Q (y) = - 2 l ∫ l 0 [ y sin ( π l x ) -1 ] sin ( π l x

Figure 1 :

 1 Figure1: Half-wavelength imperfection in the line of a strut resting on an elastic foundation. The imperfection is W 0 , the compressive load is P and the buckling displacement is W . Springs have a nonlinear elastic force-displacement relationship P (W ).

Figure 2 :

 2 Figure 2: (a) Restoring force. (b) Dimensionless restoring force. Bi-linear prole (solid line), exponential prole (dashed line). Dotted lines: limiting plateau (horizontal line), mobilization (vertical line).

Figure 3 :

 3 Figure3: Equilibrium paths predicted by (a) the piecewise solution theory, (b) the Galerkin method, (c) the numerical resolution of (8), case of a bi-linear restoring force. On each graph, the equilibrium paths are plotted (from top to bottom) for a 0 = 0, a 0 = 0.0238, a 0 = 0.595, a 0 = 1.19 and l = 3. Dotted lines: critical load (upper line), Euler load (lower line), mobilization (vertical line).

Figure 4 :

 4 Figure 4: Deection predicted by the piecewise solution theory. (b) is a zoom of (a). Function w 1 (solid line), function w 2 (dashed line), function w 3 (dash-dotted line). The functions w 1 and w 2 are connected at x = x 1 ≈ 0.618 (rst vertical dotted line). The functions w 2 and w 3 are connected at x = x 2 ≈ 2.382 (second vertical dotted line). The horizontal dotted line represents the mobilization. l = 3, a 0 = 0, λ = 1.75.

Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :

 56789 Figure5: Equilibrium paths predicted by the Galerkin method. Bi-linear restoring force (solid line), exponential restoring force (dashed lines). Galerkin method (upper dashed line), Galerkin method with the Maltby and Calladine (1995b) assumption (lower dashed line). Dotted lines: critical load (upper line), Euler load (lower line), mobilization (vertical line). l = 3, a 0 = 0 (a), a 0 = 0.0238 (b), a 0 = 0.595 (c), a 0 = 1.19 (d).

Q

  0. Therefore, the function Q can be written asQ (y) = -2H (y -(t) -1] sin (t) H (y sin (t) -1) dt. (A.5)The argument of the Heaviside function under the integral sign equals 0 whensin (t) = 1 y , that is to say for t = t 1 = arcsin (1 y) and t = π -t 1 . It comes that the function Q is non zero between t 1 and π -t 1 from (39) is recovered.The terms W n are classical to calculate. For n = 2p and n = 2p + 1 it yields appearing in (B.5) into odd and even indices gives Q (y) = Σ 1 (y) + Σ 2 (y) with Σ 1 (y) =and I1 the modied Struve and Bessel functions of pa- rameter 1. FinallyQ (y) = Σ 1 (y) + Σ 2 (y) = 2 [ I1 (y) -L1 (y)] -y,(B.11) and the result from (40) is recovered.

Appendix A. Function Q, bi-linear restoring force

The aim of this Appendix is to calculate the function Q appearing in (39) when the bi-linear restoring force is considered.

Introduce H the Heaviside function dened as

restoring force

The aim of this Appendix is to calculate the function Q appearing in (40) when the exponential restoring force is considered. For this calculus, we recall that the modied Struve and Bessel functions of parameter 1 can be expended as power series

Substituting the exponential function in (11) by its power series yields