Yongqi Wang
email: wang@fdy.tu-darmstadt.de

Wei Wu

Unsteady Flow of a Fourth Grade Fluid Due to an Oscillating Plate

Keywords: Fourth-grade fluid, non-Newtonian fluid, magnetohydrodynamic fluid, time-dependent flow

The governing non-linear high-order, sixth-order in space and third-order in time, differential equation is constructed for the unsteady flow of an incompressible conducting fourth-grade fluid in a semi-infinite domain. The unsteady flow is induced by a periodically oscillating two-dimensional infinite porous plate with suction/blowing, located in a uniform magnetic field. It is shown that by augmenting additional boundary conditions at infinity based on asymptotic structures and transforming the semi-infinite physical space to a bounded computational domain by means of a coordinate transformation, it is possible to obtain numerical solutions of the nonlinear magnetohydrodynamic equation. In particular, due to the unsymmetry of the boundary conditions, in numerical simulations non-central difference schemes are constructed and employed to approximate the emerging higher-order spatial derivatives. Effects of material parameters, uniform suction or blowing past the porous plate, exerted magnetic field and oscillation frequency of the plate on the time-dependent flow, especially on the boundary layer structure near the plate, are numerically analysed and discussed. The flow behaviour of the fourth-grade non-Newtonian fluid is also compared with those of the Newtonian fluid.

Introduction

The inadequacy of the classical Navier-Stokes theory to describe rheologically complex fluids such as polymer solutions, blood, paints, certain oils and greases, has led to the

A c c e p t e d m a n u s c r i p t

development of several theories of non-Newtonian fluids. In recent years, amongst the many models which have been used to describe the non-Newtonian behaviours exhibited by certain fluids, the fluids of differential type [START_REF] Truesdell | The Non-Linear Field Theories of Mechanics[END_REF] have received special attention, as well as much controversy, see e.g. [START_REF] Dunn | Fluids of differential type: critical review and thermodynamic analysis[END_REF] for a complete discussion of the relevant issues. The fluids of second and third grade, which form a subclass of the fluids of the differential type, have been studied successfully in various types of flow situations. We mention here some of the studies such as Akyildiz [START_REF] Akyildiz | A note on the flow of a non-Newtonian fluid film[END_REF], Benharbit and Siddiqui [START_REF] Benharbit | Certain solutions of the planar motion of a second grade fluid for steady and unsteady cases[END_REF], Erdogan [START_REF] Erdogan | Plane surface suddenly set in motion in a non-Newtonian fluid[END_REF], Gupta and Massoudi [START_REF] Gupta | Flow of a generalized second grade fluid between heated plates[END_REF], Hayat et. al. [START_REF] Hayat | Fluctuating flow of a third grade fluid on a porous plate in a rotating medium[END_REF][START_REF] Hayat | Peristaltic transport of a third-order fluid in a circular cylindrical tube[END_REF] and Rajagopal [START_REF]On the creeping flow of second order fluid[END_REF].

Although the second-grade model is found to predict the normal stress differences, it does not properly respond to shear thinning or thickening due to its constant apparent shear viscosity. For this reason, some experiments may be well described by the fluids of grade three or four [START_REF] Beavers | The rotating rod viscometer[END_REF]4,[START_REF] Joseph | The free surface on a liquid between cylinders rotating at different speeds[END_REF]. Keeping this fact in view, the aim of the present analysis is to venture further into the regime of fourth-grade fluids. Only fairly scarce literatures are available [START_REF] Hayat | Flow of fourth-order fluid[END_REF][START_REF] Kaloni | A note on the flow of a viscoelastic fluid between eccentric disks[END_REF][START_REF] Wang | Hydromagnetic rotating flow of a fourth-order fluid past a porous plate[END_REF]. No attempt has been made to model unsteady flows of fourth-grade fluids in the context of magnetohydrodynamics. In the present analysis, such an attempt is made. Literature survey reveals no previous attempts at studying such a high-order nonlinear partial differential equation even in the absence of plate porosity and magnetic field. The modelled partial differential equation is in a generalized form and is a significant contribution to understand the behaviour of fourthgrade fluids both from physical and mathematical standpoints. We examine numerically the unsteady flow in a semi-infinite space caused by the periodic motion of an infinite plate. In solving this problem the additional boundary conditions that are imposed at infinity are of particular interest. Furthermore, by a coordinate transformation, the flow problem needs to be solved only in a bounded computational domain instead of the investigated semi-infinite physical space. Due to the inequality of the numbers of the boundary conditions at both boundaries, i.e., much more boundary conditions at infinity than at the plate, in numerical simulations non-central difference schemes are constructed and employed to approximate the emerging higher-order spatial derivatives. These difference schemes are inclined toward the boundary at infinity and but still of second-order accuracy. The effects of the material parameters of the fourth-grade fluid, magnetic field and suction through the plate on the velocity spatial distribution and its time series are investigated and compared with the Newtonian fluid.

Equation of Motion

The basic field equations governing the flow of an incompressible conducting fluid are

div v = 0, (1)
ρ dv dt = div T + J × B, ( 2
)
where ρ is the density, v the velocity, T the Cauchy stress tensor, J the current density and B the total magnetic field so that B = B 0 + b, b is the induced magnetic field. d/dt denotes the material time derivative.

The above system of equations will be closed by a constitutive equation for the stress tensor T . We examine the flow of an incompressible fourth-grade fluid described by the

with n = 4, in which

S 1 = µA 1 , (4)
S 2 = α 1 A 2 + α 2 A 2 1 , (5)
S 3 = β 1 A 3 + β 2 (A 1 A 2 + A 2 A 1 ) + β 3 (tr A 2 1 )A 1 , (6)
S 4 = γ 1 A 4 + γ 2 (A 3 A 1 + A 1 A 3 ) + γ 3 A 2 2 + γ 4 (A 2 A 2 1 + A 2 1 A 2 ) +γ 5 (tr A 2 )A 2 + γ 6 (tr A 2 )A 2 1 + [γ 7 tr A 3 + γ 8 tr (A 2 A 1 )] A 1 , (7)
where µ is the coefficient of shear viscosity, and

α i (i = 1, 2), β i (i = 1, 2, 3), γ i (i = 1, 2, . . . , 8
) are material parameters. The Rivlin-Ericksen tensors A n are defined by the recursion relation:

A n = dA n-1 dt + A n-1 (grad v) + (grad v) T A n-1 , n > 1, (8)
A 1 = (grad v) + (grad v) T .
We note that when γ i = 0 (i = 1, 2, . . . , 8) the fourth-grade model reduces to the third-grade model, i.e., n = 3 in (3), while when β i = 0 (i = 1, 2, 3) and γ i = 0 (i = 1, 2, . . . , 8) the model reduces to a second-grade fluid, corresponding to the case of n = 2 in [START_REF] Benharbit | Certain solutions of the planar motion of a second grade fluid for steady and unsteady cases[END_REF]. If all of the high-order material parameters vanish, i.e., α i = 0 (i = 1, 2), β i = 0 (i = 1, 2, 3) and γ i = 0 (i = 1, 2, . . . , 8), the model reduces to the classical Navier-Stokes fluid.

We consider the unsteady flow of the incompressible, conducting fourth-grade fluid, excited by the longitudinal oscillation of an infinite porous plate with suction/blowing. We employ a spatial coordinate system with the origin on the plate, x-axis along the plate in the direction of the plate movement, and y-axis perpendicular to the plate directed toward the fluid, which is in contact with the plate and occupies the whole semi-infinite region y ≥ 0. A uniform magnetic field B 0 = B 0 k (k is a unit vector normal to the plate in the y-direction) is applied to the fluid system. In the lowmagnetic-Reynolds-number approximation, in which the induced magnetic field b can be ignored, the magnetic body force becomes (see e.g. [START_REF] Wang | On nonlinear magnetohydrodynamic problems of an Oldroyd 6-constant fluid[END_REF])

J × B = -σB 2 0 v, (9)
in which σ is the electrical conductivity of the fluid. Under these conditions, the flow is independent of x. Thus, it follows from the equation of continuity (1) that ∂v ∂y = 0, [START_REF] Hayat | Flow of fourth-order fluid[END_REF] hence the velocity component in y-direction, v, is a function of time only. We take

v = const = -v 0 , (11)

A c c e p t e d m a n u s c r i p t

where v 0 is the suction (v 0 > 0) and blowing (v 0 < 0) velocity normal to the plate. We shall seek a velocity field of the form v = (u, v, w) = (u(y, t), -v 0 , 0) ,

where u, v and w are the x, y and z components of the velocity field. Substituting ( 4)-( 8), [START_REF] Joseph | The free surface on a liquid between cylinders rotating at different speeds[END_REF] and ( 12) into (3), then the resulting constitutive relations together with (9) again into the momentum equation ( 2), and subsequently eliminating the emerging pressure gradients ∂p/∂x and ∂p/∂y by cross differentiation, after some intricate computations, we can obtain a differential equation

ρ ∂ 2 u ∂y∂t -v 0 ∂ 2 u ∂y 2 = µ ∂ 3 u ∂y 3 -σB 2 0 ∂u ∂y + α 1 ∂ 4 u ∂y 3 ∂t -v 0 ∂ 4 u ∂y 4 + β 1 ∂ 5 u ∂y 3 ∂t 2 -2v 0 ∂ 5 u ∂y 4 ∂t + v 2 0 ∂ 5 u ∂y 5 +6(β 2 + β 3 ) ∂ ∂y ∂u ∂y 2 ∂ 2 u ∂y 2 + γ 1 ∂ 6 u ∂y 3 ∂t 3 -3v 0 ∂ 6 u ∂y 4 ∂t 2 + 3v 2 0 ∂ 6 u ∂t∂y 5 -v 3 0 ∂ 6 u ∂y 6 +(6γ 2 + 2γ 3 + 2γ 4 + 2γ 5 + 6γ 7 + 2γ 8 ) ∂ 2 ∂y 2 ∂u ∂y 2 ∂ 2 u ∂y∂t -v 0 ∂ 2 u ∂y 2 . ( 13
)
Its dimensionless form can be written as

∂ ∂ t ∂ ū ∂ ȳ -ᾱ1 ∂ 3 ū ∂ ȳ3 + 2 β1 v0 ∂ 4 ū ∂ ȳ4 -3γ 1 v2 0 ∂ 5 ū ∂ ȳ5 -2γ ∂ 2 ū ∂ ȳ2 2 + ∂ ū ∂ ȳ ∂ 3 ū ∂ ȳ3 ∂ ∂ t ∂ ū ∂ ȳ -4γ ∂ ū ∂ ȳ ∂ 2 ū ∂ ȳ2 ∂ ∂ t ∂ 2 ū ∂ ȳ2 -γ ∂ ū ∂ ȳ 2 ∂ ∂ t ∂ 3 ū ∂ ȳ3 + ∂ 2 ∂ t2 3γ 1 v0 ∂ 4 ū ∂ ȳ4 -β1 ∂ 3 ū ∂ ȳ3 -γ1 ∂ 3 ∂ t3 ∂ 3 ū ∂ ȳ3 = - B ∂ ū ∂ ȳ + v0 ∂ 2 ū ∂ ȳ2 + ∂ 3 ū ∂ ȳ3 -ᾱ1 v0 ∂ 4 ū ∂ ȳ4 + β1 v2 0 ∂ 5 ū ∂ ȳ5 + β 2 ∂ ū ∂ ȳ ∂ 2 ū ∂ ȳ2 2 + ∂ ū ∂ ȳ 2 ∂ 3 ū ∂ ȳ3 -γ 1 v3 0 ∂ 6 ū ∂ ȳ6 -γv 0 2 ∂ 2 ū ∂ ȳ2 3 + 6 ∂ ū ∂ ȳ ∂ 2 ū ∂ ȳ2 ∂ 3 ū ∂ ȳ3 + ∂ ū ∂ ȳ 2 ∂ 4 ū ∂ ȳ4 , (14)
where the emerging dimensionless variables and parameters are defined by

ū = u U , ȳ = Uy ν , v0 = v 0 U , t = ρU 2 µ t, B = νσB 2 0 ρU 2 , ᾱ1 = α 1 U 2 ρν 2 , β1 = β 1 U 4 ρν 3 , γ1 = γ 1 U 6 ρν 4 , (15) β
= 6(β 2 + β 3 )U 4 ρν 3 , γ = 2 (3γ 2 + γ 3 + γ 4 + γ 5 + 3γ 7 + γ 8 ) U 6 ρν 4 ,

A c c e p t e d m a n u s c r i p t

in which U is a reference velocity and ν = µ/ρ is the fluid kinematic viscosity. We note that ( 14) is a high-order partial differential equation (sixth-order in space and third-order in time) and thus requires additional boundary conditions except for the usual boundary conditions for a Newtonian fluid [START_REF] Rajagopal | An exact solution for the flow of a non-Newtonian fluid past an infinite porous plate[END_REF]. As we are solving the problem in an unbounded domain, it is possible to augment the boundary conditions by enforcing additional asymptotic structures at infinity. Here, we augment the additional boundary conditions by requiring that the several high-order derivatives in space vanish at infinity. Similar additional initial conditions are also postulated. Their different choices may have influence only on the numerical results for the initial time steps. We will solve [START_REF]On the creeping flow of second order fluid[END_REF] in the domain ȳ ∈ [0, ∞) with the following boundary conditions

ū = sin(ω t) at ȳ = 0, ū = 0, ∂ n ū ∂ ȳn = 0 (n = 2, 3, 4, 5) at ȳ → ∞, (16)
and the initial conditions

ū = 0, ∂ n ū ∂ tn = 0, (n = 2, 3) at t = 0. ( 17
)
For the boundary condition at the plate, ( 16) 1 , we assume a non-slip condition and the plate oscillates in a sinusoidal form, in which ω is the corresponding dimensionless frequency. It should be pointed out that the additional boundary conditions at ȳ → ∞ in ( 16) 2 are specified for n = 2, 3, 4, 5 instead of those for n = 1, 2, 3, 4. This choice is for no other reason than the consideration of numerical simulations, because the physical value at the first interior grid point near the boundary ȳ → ∞ would be determined directly by the boundary conditions instead of the governing differential equation, if the other choice with n = 1, 2, 3, 4 would be employed. It is also the similar case for the additional initial conditions (17) 2 .

For simplicity, in the following we will drop the bars of the dimensionless variables.

Numerical Method

Approximate solutions can be obtained by truncating the semi-infinite domain. In this approximation, the semi-infinite domain y ∈ [0, ∞) is replaced by a finite domain y ∈ [0, H] and the boundary conditions at infinity y → ∞ are enforced at y = H. The length H is chosen in such a manner that any further increase in its value does not significantly alter the solution. It should be emphasized here that proper care should be exercised in the determination of the solution because it is well known that even in the case of a Navier-Stokes fluid spurious solutions are possible due to the truncation of an infinite domain to a bounded one [START_REF] Mcleod | On the uniqueness of flow of the Navier-Stokes fluid due to a stretching boundary[END_REF]. In order to avoid such possible spurious solutions, in this paper the coordinate transformation η = 1/(y + 1) is applied to transform the semi-infinite physical domain y ∈ [0, ∞) to a finite calculation domain η ∈ [0, 1], by the following transformation relations

y = 1 η -1, ∂ ∂y = -η 2 ∂ ∂η , ∂ 2 ∂y 2 = η 4 ∂ 2 ∂η 2 + 2η 3 ∂ ∂η ,

A c c e p t e d m a n u s c r i p t

∂ 3 ∂y 3 = -η 6 ∂ 3 ∂η 3 -6η 5 ∂ 2 ∂η 2 -6η 4 ∂ ∂η , ∂ 4 ∂y 4 = η 8 ∂ 4 ∂η 4 + 12η 7 ∂ 3 ∂η 3 + 36η 6 ∂ 2 ∂η 2 + 24η 5 ∂ ∂η , (18)
∂ 5 ∂y 5 = -η 10 ∂ 5 ∂η 5 -20η 9 ∂ 4 ∂η 4 -120η 8 ∂ 3 ∂η 3 -240η 7 ∂ 2 ∂η 2 -120η 6 ∂ ∂η , ∂ 6 ∂y 6 = η 12 ∂ 6 ∂η 6 + 30η 11 ∂ 5 ∂η 5 + 300η 10 ∂ 4 ∂η 4 + 1200η 9 ∂ 3 ∂η 3 + 1800η 8 ∂ 2 ∂η 2 + 720η 7 ∂ ∂η .
With these transformations, the governing differential equation ( 14) can be rewritten in the η coordinate. Here, we avoid writing its explicit form because of its complexity.

The boundary conditions ( 16) in terms of η can be rewritten as

u = sin(ωt), at η = 1, u = 0, ∂ n u ∂η n = 0 (n = 2, 3, 4, 5), at η = 0. ( 19
)
In the subsequent simulations for most cases the dimensionless oscillation frequency of the plate, ω, is simply taken as ω = 1, but we will also discuss the effect of various values of ω on the time series of velocity and the structure of the boundary layer near plate.

The differential equation ( 14) is a nonlinear high-order differential equation. We can solve this initial-boundary value problem ( 14), ( 17) and ( 19) by finite difference method. We can discretize [START_REF]On the creeping flow of second order fluid[END_REF] (in fact its form in terms of η) for M uniformly distributed discrete points η = (η 1 , η 2 , . . . , η M ) ∈ (0, 1) with a space grid size of ∆η = 1/(M + 1) and time steps t = (t (1) , t (2) , ..., t (n) , ...) = (∆t, 2∆t, ..., n∆t, ...), where ∆t is the time step size. Due to the inequality of the numbers of the boundary conditions at both boundaries, i.e., five conditions at η = 0, but only one at η = 1, we cannot use central differences to approximate the high-order derivatives emerging in ( 14) and (19) (if the order > 2). The following finite difference schemes are constructed and employed for different order derivatives. Except for the first and second derivatives they are not central differences, but inclined towards the boundary η = 0, where more boundary conditions are specified; however they are still of second-order accuracy. Here, we omit the laborious derivation and only list the formulas for the derivatives at point η = η j :

∂u ∂η | η=η j = u j+1 -u j-1 2(∆η) + o(∆η 2 ), ∂ 2 u ∂η 2 | η=η j = u j+1 -2u j + u j-1 (∆η) 2 + o(∆η 2 ), ∂ 3 u ∂η 3 | η=η j = 3u j+1 -10u j + 12u j-1 -6u j-2 + u j-3 2(∆η) 3 + o(∆η 2 ), (20) ∂ 4 u ∂η 4 | η=η j = 2u j+1 -9u j + 16u j-1 -14u j-2 + 6u j-3 -u j-4 (∆η) 4 + o(∆η 2 ),
A c c e p t e d m a n u s c r i p t

∂ 5 u ∂η 5 | η=η j = 5u j+1 -28u j + 65u j-1 -80u j-2 + 55u j-3 -20u j-4 + 3u j-5 2(∆η) 5 + o(∆η 2 ), ∂ 6 u ∂η 6 | η=η j = 3u j+1 -20u j + 57u j-1 -90u j-2 + 85u j-3 -48u j-4 + 15u j-5 -2u j-6 (∆η) 6
+o(∆η 2 ), in which o(∆η 2 ) denotes quantities which are higher order small than ∆η 2 . For the fifth-order derivative emerging in the boundary condition at infinity η = 0 (j = 0) and the sixth-order derivative emerging in the differential equation [START_REF]On the creeping flow of second order fluid[END_REF] for the first interior grid point near the boundary of infinity η = 0 (j = 1) only difference approximations with first-order accuracy are used to avoid additional unknown imaginary boundary values, i.e.,

∂ 5 u ∂η 5 | η=η j = u j+1 -5u j + 10u j-1 -10u j-2 + 5u j-3 -u j-4 (∆η) 5 + o(∆η), (j = 0) (21) ∂ 6 u ∂η 6 | η=η j = u j+1 -6u j + 15u j-1 -20u j-2 + 15u j-3 -6u j-4 + u j-5
(∆η) 6 +o(∆η) (j = 1).

Similarly, a temporally afterwards form for the third-order time derivative is used

∂ 3 u ∂t 3 | t=t (n) = 3u (n+1) -10u (n) + 12u (n-1) -6u (n-2) + u (n-3) 2(∆t) 3 + o(∆t 2 ), (22)
whilst for the first-and second-order time derivative still traditional central difference schemes are employed.

In so doing, for each time step a linear algebraic equation system emerges. Its coefficient matrix has a bandwidth of eight elements. This equation system can be solved e.g. by Gaussian elimination. The values of the velocity field u at the boundary points η 0 = 0 (j = 0), η M +1 = 1 (j = M + 1) and the imaginary points outside the domain (j = -1, -2, -3, -4) can be obtained by using the boundary conditions (19).

Numerical Results and Discussions

We compare the velocity profiles for two kinds of fluids: the Newtonian fluid, for which α 1 = 0, β = β 1 = 0 and γ = γ 1 = 0, and the full fourth-grade non-Newtonian fluid, in which for most cases we choose α 1 = β 1 = β = γ 1 = γ = 1 without physical reasons. In addition the influences of other different values of the material parameters are also investigated. In order to clearly observe the structure of the boundary layer near the plate, all numerical results are depicted in the y-coordinate only for a truncated physical domain y ∈ [0, 9] (corresponding to η ∈ [0.1, 1]) near the plate, after finishing the simulations in the whole computational domain η ∈ [0, 1], corresponding to the physical domain z ∈ [0, ∞). The time series of the flow velocity in several different positions (different distances from the plate) for the first seven periods ωt ∈ [0, 14π] are displayed in Fig. 1 for the fourth-grade non-Newtonian fluid and in Fig. 2 for the Newtonian fluid, respectively. Obviously, periodic flows can be reached almost immediately after the set-up of the flow conditions, at most after two periods. Due to the no-slip boundary condition at the solid boundary the fluid near the plate oscillates together with the plate in the same phase and the amplitude. The velocity amplitude decays rapidly with the increase of the distance from the plate. Whilst the flow of the fourth-grade non-Newtonian fluid in the whole flow domain oscillates approximately in phase with the driving plate movement (Fig. 1), a phase lag behind the excitation occurs for the Newtonian fluid which increases with the increase of the distance away from the plate (Fig. 2). Actually, these phase difference

A c c e p t e d m a n u s c r i p t

depends on the ratio of viscous to inertial forces for the investigated fluids. If the apparent viscosity of the fluid is sufficiently large, i.e., the viscous forces dominate, the flow can simultaneously adapt to the current driving force and oscillates approximately with the same phase in the whole flow domain. This is the case for the fourth-order non-Newtonian fluid. For the Newtonian fluid displayed in Fig. 2 the inertial forces dominate and the flow slowly reacts to the changes of the plate movement, hence in the flow velocity a phase lag behind the excitation occurs. This time delay increases with the increase of the distance from the plate. From Figs. 3 and4 it can be clearly seen that the investigated fourth-grade non-Newtonian fluid shows a larger apparent viscosity than the corresponding Newtonian fluid. In both figures the velocity profiles at two different phases ωt = (2n + 1/2)π (wave peak) and ωt = (2n + 3/2)π (wave trough) within a period ωt ∈ [2nπ, 2(n + 1)π] with a sufficiently large value of n, for which the periodic flow has been reached (e.g. after several periods), are represented for various values of the second-order (α 1 ) in Fig. 3 and third-order material parameters (β, β 1 ) in Fig. 4 of the non-Newtonian fluid, respectively. For comparison, the corresponding velocity profiles of the Newtonian fluid are also displayed (solid lines in the figures). In general, the non-Newtonian fluid exhibits a much thicker boundary layer, indicating a larger apparent viscosity than that of the Newtonian fluid. Increasing these second-and third-order material parameters (α 1 in Fig. 3 andβ, β 1 in Fig. 4) causes further thickening of the boundary layer. Furthermore, increasing the fourth-order material parameters, γ and γ 1 , indicates a similar tendency, that is for briefness not shown here.

The numerical results obtained by using various suction velocities v 0 through the plate are illustrated in Fig. 5. The alteration tendencies of the velocity distributions by varying v 0 are practically converse for the Newtonian fluid and the fourth-grade fluid. For the Newtonian fluid (Fig. 5a), when the suction velocity v 0 increases, as one may expect, the boundary layer near the plate (y = 0) tends to become thinner. On the contrary, for the fourth-order fluid (Fig. 5b) the effect of increase in the suction velocity increases the boundary layer thickness. Test computations with other different highorder material parameters show similar tendency, except when all high-order material parameters are chosen to be very small, for which the fourth-grade fluid behaves like a Newtonian fluid. In general, the boundary layer for the Newtonian fluid is much thinner than that for the fourth-order fluid, especially for the case with a suction through the porous plate. For a blowing through the plate, the fluids behave qualitatively opposite to a suction velocity. In this case, with the increase of the blowing velocity, for the Newtonian fluid the shear boundary layer becomes thicker rapidly, whilst for the fourthgrade fluids the thickness of boundary layers is not so sensitive to variations in the blowing, and only slight thinning of the boundary layer occurs. To save space, the corresponding results with a blowing are not graphically displayed.

The dependence of the distributions of the flow velocity on the magnetic field is shown in Fig. 6. An applied magnetic field tends to restrict the shearing to a thinner boundary layer near the plate not only for the investigated fourth-order non-Newtonian fluid but also for the Newtonian fluid, although this thinning seems to be more obvious for the fourth-order fluid than that of the Newtonian fluid. The reason for this thinning is that the magnetic field provides a resistance to the flow and hence decrease the flow

α 1 = γ = γ 1 = 1.
Here n is chosen to be sufficiently large to assure that the periodic flow has been reached. The other parameters are chosen as B = 0, v 0 = 0 and ω = 1. velocity.

To our surprise, an applied magnetic field can put forward the oscillation phase in the domain away from the plate. This can be extracted from the time series of velocity displayed in Fig. 7 for a fourth-grade fluid with three different values of the magnetic field B = 0, 1, 2 at two different distances from the plate y = 1, 2, respectively. Even a phase lead before the driving plate movement occurs which increases with the increase of the distance away from the plate (by comparing the above and below panels) and the magnetic field strength B, as it can be easily identified in Fig. 7, in which for B = 0 (solid lines) the flow is almost in the same phase as the driving plate movement (as shown in Fig. 1). As expected, with the increase of B, the amplitude of the flow decreases because the magnetic force is a resistance to the flow. The influence of B on the velocity time series of a Newtonian fluid has a similar tendency and is not repeated here.

As we have pointed out before, to some extent increasing the high-order material parameters corresponds to the increase of the apparent viscosity of the fourth-grade fluid. This effect can also be seen from the time series of flow velocity at a fixed position y = 2 with different values of the second-order material parameter α 1 , displayed in Fig. 8. With the increase of α 1 , the amplitude of the flow velocity at a given distance from the plate increases, indicating the boundary layer becomes thickener. For cases in which a phase difference between the flow and the driving force occurs, increasing the high-order material parameters will cause the fluid to react to the driving force more rapidly and hence reduce the possible phase difference (lag or lead). This is the case for B = 2 displayed in Fig. 8 (below panel), in which the phase lead decreases with the increase of α 1 , whilst for the case of B = 0 (the above panel of Fig. 8) the flow velocities with the different values of α 1 are nearly in the same phase with the driving movement, hence such a phase change by varying α 1 is invisible. Numerical results by increasing the other high-order material parameters β, β 1 or γ, γ 1 show similar influences, but are not graphically presented here.

In the numerical results displayed before we took the dimensionless oscillation frequency of the plate as ω = 1. It is fairly interesting to observe the effect of various values of ω on the time series of the velocity and the structure of the boundary layer near the oscillating plate. The corresponding numerical results for three different values of ω, ω = 0.3, 1, 3, are graphically shown in Figs. 9 and 10, respectively. It is obvious that, for the investigated forth-order non-Newtonian fluid, the boundary layer near the oscillating plate becomes thicker when the dimensionless oscillation frequency ω deviates from its value of ω = 1, as seen in Fig. 9. The case with ω = 1 corresponds to the thinnest boundary layer. As shown in Fig. 1, for the investigated fourth-order fluid in the absence of the suction/blowing through the plate (v 0 = 0) and the magnetic field (B = 0), the flow for ω = 1 in the whole flow domain oscillates approximately in phase with the driving plate movement. It is also the case shown in Fig. 10 with ω = 1 (dashed line). For ω < 1 (solid line) a phase lag after the driving plate movement exists, whilst for ω > 1 even a phase lead occurs.

Final Remarks

The time-dependent flow field of a magnetohydrodynamic (MHD) fluid of grade four is investigated. The velocity field is governed by a non-linear higher-order, sixth-order in space and third-order in time, partial differential equation. The governing equation is sufficiently general to incorporate the effects of an applied magnetic field and allows the study of blowing/suction at the boundary. Numerical results for the unsteady MHD flow caused by a periodic motion of the infinite two-dimensional plate are presented,

A c c e p t e d m a n u s c r i p t

for which the surrounding fluid occupies a semi-infinite space. This offers difficulties due to the unbounded domain and the lack of sufficient physical boundary conditions, even in the case of the second-or third-grade fluids, more boundary conditions are necessary in addition to the usual ones for a Newtonian fluid. In the unbounded domain the additional boundary conditions are augmented by enforcing additional asymptotic structure at infinity, although these are relatively weak and physically plausible, see [START_REF] Rajagopal | An exact solution for the flow of a non-Newtonian fluid past an infinite porous plate[END_REF]. In order to circumvent the difficulty with the infinite flow domain, a coordinate transformation is employed, which causes some complicated derivations. Furthermore, to match the unsymmetrical boundary conditions, non-central difference schemes for the high-order spatial derivatives are constructed and used. The purpose of our analysis is to investigate flows of a fourth-grade fluid with a view towards understanding its response characteristics. The flow fields of the fourth-order non-Newtonian fluid are compared with those of a linear Newtonian fluid. Numerical results show that phase differences in flow fields may appear, which may be a phase lag or lead, depending on the chosen parameters. Increasing the strength of the applied magnetic field decreases the boundary-layer thickness, that is valid for both the Newtonian fluid and the fourth-order non-Newtonian fluid, while the effect of alteration in suction velocity through the porous plate shows obvious differences between the Newtonian and fourth-grade non-Newtonian fluids. In general, the fourth-order non-Newtonian fluid exhibits a thicker boundary layer than the Newtonian fluid, and increasing the highorder material parameters of the non-Newtonian fluid causes further thickening of the boundary layer.

It should be pointed out that the quantitative behaviour of high-order fluids is strongly dependent on the choice of the material parameters. At the present time there appears to be no experimental data determining these parameters, particularly for these material parameters in a theory of fourth-order fluids. In this paper, only the simplest prototype flow problem, containing the ingredients essential to understanding the phenomena, is presented. It is only a step towards understanding the behaviour of fourth-order fluids. A good deal of refinement and extension is required, especially the determination of physically reasonable high-order material parameters, before more quantitative answers can be expected.

Figure 1 :Figure 2 :

 12 Figure 1: Time series of the flow velocity at several different distances from the plate y = 0, 0.25, 1.0, 2.0 within the first seven periods ωt ∈ [0, 14π] for the fourth-grade non-Newtonian fluid with the chosen material parameters α 1= β = β 1 = γ = γ 1 = 1.The other parameters are chosen as B = 0, v 0 = 0 and ω = 1.

Figure 3 :Figure 4 :

 34 Figure 3: Profiles of the flow velocity u(y, t) at two different time slices ωt = (2n+1/2)π (left panel) and ωt = (2n + 3/2)π (right panel) for the Newtonian fluid (solid lines) and the fourth-grade non-Newtonian fluid with various values of the second-order material parameter, α 1 = 1, 2, 5, while the other high-order material parameters are fixed, β =β 1 = γ = γ 1 = 1.Here n is chosen to be sufficiently large to assure that the periodic flow has been reached. The other parameters are chosen as B = 0, v 0 = 0 and ω = 1.

Figure 5 :

 5 Figure 5: Profiles of the flow velocity u(y) at the time slice ωt = (2n + 1/2)π with a sufficiently large value of n for the Newtonian fluid (left panel) and the fourth-grade fluid α 1 = β = β 1 = γ = γ 1 = 1 (right panel) with different suction velocities in the absence of the magnetic field B = 0. The dimensionless oscillation frequency of the plate is fixed as ω = 1.

Figure 6 :

 6 Figure 6: Profiles of the flow velocity u(y) at the time slice ωt = (2n + 1/2)π with a sufficiently large value of n for the Newtonian fluid (left panel) and the fourth-grade fluid α 1 = β = β 1 = γ = γ 1 = 1 (right panel) with different values of the magnetic field B = 0, 1, 2 in the absence of suction/blowing through the plate v 0 = 0. The dimensionless oscillation frequency of the plate is fixed as ω = 1.

ωtFigure 7 :

 7 Figure 7: Time series of the flow velocity at two different distances from the plate, y = 1 (above panel) and y = 2 (below panel), within the first seven periods ωt ∈ [0, 14π] for a non-Newtonian fluid with the chosen material parameters α 1 = β = β 1 = γ = γ 1 = 1. For the magnetic field strength, three different values are taken: B = 0 (solid lines), B = 1 (dashed lines) and B = 2 (dotted lines). No suction/blowing past the plate exists v 0 = 0. The dimensionless oscillation frequency of the plate is fixed as ω = 1.

Figure 8 :

 8 Figure 8: Time series of the flow velocity at a fixed distance from the plate of y = 2 within the first seven periods ωt ∈ [0, 14π] for fourth-grade non-Newtonian fluids with different values of the second-order material parameter α 1 : α 1 = 1 (solid line), α 1 = 2 (dashed line) and α 1 = 5 (dotted line). The values of the other material parameters are fixed: β = β 1 = γ = γ 1 = 1. The magnetic field strength is chosen as B = 0 (above panel) and B = 2 (below panel). No suction/blowing past the plate exists: v 0 = 0. The dimensionless oscillation frequency of the plate is fixed as ω = 1.

Figure 9 : 3 ωtFigure 10 :

 9310 Figure 9: Profiles of the flow velocity u(y, t) at two different time slices ωt = (2n+1/2)π (left panel) and ωt = (2n + 3/2)π (right panel) for a fourth-grade non-Newtonian fluid with the chosen material parameters α 1= β = β 1 = γ = γ 1 = 1.Here n is chosen to be sufficiently large to assure that the periodic flow has been reached. Three different values of the dimensionless frequency ω are considered: ω = 0.3, 1, 3. The other parameters are chosen as B = 0, v 0 = 0.

A c c e p t e d m a n u s c r i p t


Stefano Lenci
email: lenci@univpm.it

Giuseppe Rega
email: giuseppe.rega@uniroma1.it

Dimension reduction of homoclinic orbits of buckled beams via the nonlinear normal modes technique

Keywords: Buckled beam, infinite dimensional system, homoclinic solutions, higher-order modes, reduced order models, nonlinear normal modes

The problem of detecting the homoclinic orbits of an initially straight buckled beam is addressed. Two families of b.c. are identified and investigated in detail. For the first family, the homoclinic orbits belong to a planar invariant manifold, and are easily computed in closed form. For the second family, the manifold is no longer planar, and is detected via the nonlinear normal modes technique by obtaining approximate expressions which are sufficient to highlight the effects of the non-flatness. A hierarchy of reduced order, single degree of freedom, models is determined. These are obtained by taking into account increasing degrees of nonlinearity in the potential energy, which allow for a more and more refined computation of the homoclinic solution. The various models are compared with each other and discussed in detail, and the nonplanarity of the manifold is illustrated through examples.

Introduction

All people working on nonlinear dynamics and chaos of engineering systems become aware of the importance of homoclinic orbits, sooner or later. These are particular solutions of the governing equations connecting, in an infinite time, a given unstable equilibrium point. In spite of the fact that they are special and isolated solutions, they play a fundamental role in the organization of the global dynamics of complex systems, as repeatedly verified in practice (e.g., [START_REF] Grebogi | Crises, Sudden Changes in Chaotic Attractors, and Transient Chaos[END_REF][START_REF] Moon | Evidence for Homoclinic Orbits as a Precursor to Chaos in a Magnetic Pendulum[END_REF][START_REF] Katz | From Single Well Chaos to Cross Well Chaos: A Detailed Explanation in Terms of Manifold Intersections[END_REF]). They also play a fundamental role in the study of spatial chaos [START_REF] Thompson | Homoclinic orbits, spatial chaos and localized buckling[END_REF].

In the past, these solutions have been studied on systems with few (usually one) mechanical degrees of freedom (d.o.f.), which are governed by ordinary differential equations (ODEs) [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields[END_REF].

However, real mechanical structures are infinite dimensional, namely they are governed by partial differential equations (PDEs), and call for a more refined analysis.

The study of the homoclinic solutions in infinite dimensional systems is the main objective of

A c c e p t e d m a n u s c r i p t

this paper, where it is addressed with reference to an initially straight, planar, undamped and unforced buckled beam. While in low dimensional systems the homoclinic solutions can usually be computed, so that their existence is guaranteed, in infinite dimension their existence is initially conjectured based on system mechanical properties. In this respect, the results of this paper could be used, at least in principle, to prove a posteriori that the homoclinic solutions actually exist, would one be able to compute all terms of the series used for their representation, and to prove its convergence. We will not address this issue in this mechanically oriented paper, being satisfied of the mechanical conjecture about their existence. For the same reason, only few terms of the mentioned series will be computed, by discussing how they can be sufficient for practical purposes.

We focus on the analytical detection of homoclinic orbits, which is possible because we consider the hilltop saddle of a conservative system. More general cases have to be dealt with numerically, even for low dimensional ODEs [START_REF] Hobson | An Efficient Method for Computing Invariant Manifolds of Planar Maps[END_REF]. Yet, we believe that the case considered in this work, which aims at enlightening how to address large dimensionality and obtaining information about its main effects, can help also in view of a purely numerical approach.

The main source of difficulty is that the homoclinic solutions are expected to activate all spatial linear normal modes of the beam, thus entailing a full nonlinear modal coupling. We assume that the homoclinic orbits belong to a two dimensional invariant manifold, which is not expected to be planar in the space spanned by the linear modes. Such a manifold is the set where the dynamics we are interested in lie. Its existence is conjectured based on similar results holding for finite dimensional systems [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields[END_REF], and again it would be proved a posteriori should one succeed in showing the convergence of the series used for its detection. In turn, if the manifold is two dimensional, the dynamics embedded on it are governed by two mathematical d.o.f., corresponding to a single mechanical d.o.f.. This suggests a 'structure' to occur in the associated homoclinic solutions, which can be described by a proper reduced order model, possibly different from those obtained with classical dimensional reduction techniques, such as the linear Galerkin method. Thus, we naturally fall in the realm of dimension reduction, which is a matter of intensive on-going research (see, e.g., [START_REF] Blekhman | Vibrational Mechanics[END_REF][START_REF] Steindl | Methods for dimension reduction and their application in nonlinear dynamics[END_REF][START_REF]Dimension Reduction of Dynamical Systems: Methods, Models, Applications[END_REF]) aimed at improving the accuracy of the reduced models.

Among various competing approaches, in this paper we use the nonlinear normal modes techniques to detect the manifold, the homoclinic orbits and the reduced order model lurking in the background, thus obtaining all information of interest within a comprehensive investigation. We closely follow the approach of [START_REF] Rand | Lecture Notes on Nonlinear Vibrations[END_REF], which is similar to the center manifold reduction [START_REF] Troger | Nonlinear Stability and Bifurcation Theory[END_REF], and permits the detection of searched quantities in either exact or approximate ways. In fact, the results will be given in terms of a series, and one can consider only few terms for practical purposes. Of course, the number of required terms depends on the searched accuracy.

A c c e p t e d m a n u s c r i p t

An appealing viewpoint on the nonlinear normal modes (NNMs) of continuous systems, which somehow enlightens the relevant nature, consists in regarding them as mathematical structures accounting for time change of the system spatial shape, which instead remains fixed in linear oscillations if initially corresponding to a linear normal mode.

The idea of NNMs was introduced by Rosenberg [START_REF] Rosenberg | The normal modes of nonlinear n-degrees-of-freedom systems[END_REF], and later theoretical and practical important contributions were given by Rand [START_REF] Rand | A direct method for non-linear normal modes[END_REF] and by Vakakis [START_REF] Vakakis | Analysis and identification of linear and nonlinear normal modes in vibrating systems[END_REF]. Shaw and Pierre [START_REF] Shaw | Non-linear normal modes and invariant manifolds[END_REF][START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF] were the first to consider the NNMs as an invariant manifold tangent to the space spanned by the linear modes in phase space, which is the viewpoint followed in this paper. When applied to infinite-dimensional systems, this approach requires a previous projection of the dynamics onto the linear normal modes, so that the PDE of motion is transformed into a sequence of coupled nonlinear ODEs, which is the natural framework for the geometric approach to NNMs. Then, the machinery of invariant manifolds applies [START_REF] Xie | Normal modes of a non-linear campled-clamped beam[END_REF].

A different approach to NNMs, not pursued in this work, has also been proposed [START_REF] Pak | Free vibrations of a thin elastica by normal modes[END_REF][START_REF] Nayfeh | On direct methods for constructing nonlinear normal modes of continuous systems[END_REF] and investigated at length. It is based on directly attacking the PDE via an asymptotic development method, with no a priori projection of the solution onto linear normal modes. The coupling among modes ensues a posteriori from the infinite-dimensional analysis.

NNMs were originally developed to deal with classical nonlinear oscillations, and synthesize the effects of the nonlinearity. In this context they have been investigated in depth and developed along various research paths, with reference to various mechanical and structural systems, which include shallow arches and cables [START_REF] Lacarbonara | Resonant nonlinear normal modes. Part I: Analytical treatment for onedimensional structural systems[END_REF][START_REF] Srinil | Two-to-one resonant multi-modal dynamics of horizontal/inclined cables. Part II: Internal resonance activation, reduced order models and nonlinear normal modes[END_REF], single beams [START_REF] Xie | Normal modes of a non-linear campled-clamped beam[END_REF], frames [START_REF] Mazzilli | Non-linear normal modes of a simply supported beam: continuos system and finite-element models[END_REF], rotor blades [START_REF] Pesheck | Accurate reduced-order models for a simple rotor blade model using nonlinear normal modes[END_REF] and MEMS [START_REF] Xie | Nonlinear dynamic analysis of MEMS switches by nonlinear modal analysis[END_REF]. In this paper, they are used to detect homoclinic orbits. The authors are not aware of previous works in this direction.

The paper is organized as follows. In Sect. 2 the mechanical model and the associated partial differential equation are set, paying attention to distinguish between two families of boundary conditions giving different behaviour in terms of homoclinic solutions. The linear modal analysis is performed in Sect 3, where the projected nonlinear equations are also obtained. The nonlinear normal modes are considered in Sect. [START_REF] Thompson | Homoclinic orbits, spatial chaos and localized buckling[END_REF], where the invariant manifold containing the homoclinic orbits is determined by solving the appropriate modal equations. The dynamics on the manifold, and the homoclinic solutions, are studied in Sect. 5, where a hierarchy of approximate reduced order models is established. The paper ends (Sect. 6) with some conclusions and suggestions for further developments.

The mechanical model

The dimensionless partial differential equation governing the unforced undamped planar nonlinear dynamics of an initially straight buckled beam is where w(z,t) is the time dependent transversal displacement of the point at abscissa z∈[0,1], dot and prime represent derivatives with respect to time and space, respectively, and k is the stiffness due to membrane effects. The parameter Γ>Γ cr represents the first-order axial load (positive=compression, Γ cr =buckling critical parameter) associated with a finite prescribed end-displacement of the beam in the axial direction [START_REF] Mettler | Dynamic Buckling[END_REF][START_REF] Villaggio | Mathematical Models for Elastic Structures[END_REF].

Equation (1) describes the finite amplitude transverse dynamics of the buckled beam under the classical Kirchhoff assumption [START_REF] Kauderer | Nichtlineare Mechanik[END_REF] of neglecting the longitudinal inertia, also known in the literature as statically condensed model. It permits feasible computations and analytical developments without losing the main mechanical phenomena, at least for moderate displacements, which are implicitly assumed without further remarks. System (1) is conservative, with the Hamilton and potential functions being, respectively

H( w ,w)= 2 1         ∫ 1 0 2 dz w +V(w), V(w)=- 2 Γ         ∫ 1 0 2 ' dz w + 2 1         ∫ 1 0 2 ' ' dz w + 4 k 2 1 0 2 '         ∫ dz w . (2)
When Γ>Γ cr and irrespective of the type of transverse boundary conditions, the undeformed equilibrium position of system (1) (which now is a hilltop saddle point) has one stable and one unstable direction, and we assume that they are connected by an homoclinic orbit. As the mechanical system is symmetric, namely, if w(z,t) is a solution also -w(z,t) is a solution, the homoclinic orbits are always in couple.

Boundary conditions

To make the problem (1) well-posed, initial and boundary conditions must be added. While initial conditions are not important, as we consider global system dynamics and not specific trajectories, boundary conditions (b.c.) play a significant role in successive developments. If the constraints are supposed to be not elastic, we can assume on each boundary one condition between w=0 and w′′′+Γw′=0, and another condition between w′′=0 and w′=0. Excluding repetitions and kinematically not admissible b.c., the six basically different b.c. are reported in Table 1, where they are grouped in two different families.

Looking for an exact solution of nonlinear equation (1) in the form w(z,t)=a(t)b(z) requires b′′′′∝b′′∝b, namely, b′′=κb for an arbitrary κ. If κ is positive, this equation has solution b(z)=c 1 sinh(z√κ)+ c 2 cosh(z√κ), while for κ negative the solution is b(z)=c 1 sin(z√(-κ))+c 2 cos(z√(κ)). As we have only the two unknowns c 1 and c 2 , in general there is no hope to satisfy all four b.c.

A c c e p t e d m a n u s c r i p t

associated to [START_REF] Grebogi | Crises, Sudden Changes in Chaotic Attractors, and Transient Chaos[END_REF]. This occurs only for the two special cases of Table 1, which constitute the first family. For the hinged-hinged beam (Fig. 1a) we have b(z)=w n (z)= 2 sin(nπz), n=1,2,3,…, while for the guided-hinged beam (Fig. 1b) we have b(z)=w n (z)= 2 cos[(nπ-π/2)z], n=1,2,3,….

In both cases we have w n =0 or w n ′=0 at the boundary, and

=-κ ∫ 1 0 ' ' dz w w m n n δ nm , (3)
where κ n =-(nπ) 2 for the hinged-hinged beam, and κ n =-(nπ-π/2) 2 for the guided-hinged beam, δ nm being the Kronecker delta. This is a very important orthogonality property of this family of b.c., which will be properly exploited in the following.

Linear normal modes and projected nonlinear equations of motion

The spatial linear operator appearing in ( 1) is L(w)=w′′′′+Γw′′, and their eigenvalues Λ n (Γ) and eigenfunctions w n (z;Γ) satisfy the equation w n ′′′′+Γw n ′′=Λ n w n , whose general solution can be written in the form w n (z)=c 1 cos(αz)+c 2 sin(αz)+c 3 cosh(βz)+c 4 sinh(βz),

α= 2 Λ + Γ + Γ , β= 2 4 2 n Λ + Γ + Γ - . (4)
The eigenvalues Λ n are determined by the boundary conditions. The mathematical properties of L(w) (basically because it is self-adjoint for all b.c.) guarantee that they are real, that =δ

∫ 1 0 dz w w m n nm , -Γ =Λ ∫ 1 0 ' ' ' ' dz w w m n ∫ 1 0 ' ' dz w w m n n δ nm , (5)
and that the eigenfunctions constitute a basis of the displacement field, so that the solution can be sought in the form

w (z,t)= ∑ ( ( a . ( 6
) ∞ =1 ) ) n n n z w t
The eigenfunctions are called linear normal modes. They span a planar invariant (with respect to the linearized equations) manifold in the space of admissible functions.

For the hinged-hinged and guided-hinged beams we have, respectively, the expressions

Λ n =(nπ) 2 [(nπ) 2 -Γ], w n (z)= 2 sin(nπz), (7)
Λ n =(nπ-π/2) 2 [(nπ-π/2) 2 -Γ], w n (z)= 2 cos[(nπ-π/2)z], (8)
which are somehow special because the two cases belong to the first family of b.c., as they satisfy w n ′′=κ n w n . The relevant eigenvalues are linear functions of Γ.

For the second family, the expressions of the eigenvalues and eigenfunctions are more involved, and are better illustrated by means of the diagrams in Fig. 3. The eigenvalues are decreasing functions of the axial parameter Γ. The unique value Γ cr,n where they vanish is the n-th buckling parameter. Positive and negative values of Λ n correspond to stable and unstable linear normal mode, respectively. In the former case, n Λ is the n-th natural frequency. The homoclinic orbits we are interested in may exist only when the beam is buckled, i.e., for Γ>Γ cr,1 and Λ 1 <0, otherwise there would be no unstable direction of the hilltop saddle.

A c c e p t e d m a n u s c r i p t

Consequently, we make this assumption in the following.

Inserting the expression ( 6) in (2), and using the properties (5) we get

H( a ,a n n )= 2 1 ∑ ∞ =1 2 ) ( n n a +V(a n ), V(a n )= 2 1 ∑ ∞ = Λ 1 2 ) ( n n n a + 4 k 2 1 1 ,       ∑∑ ∞ = ∞ = n m m n m n d a a , d m,n =d n,m = , (9)
+ n a n a V ∂ ∂ = a +Λ n n a n +k   =0. ( 10
)       ∑∑ ∞ = ∞ = 1 1 , i j j i j i d a a       ∑ ∞ =1 , m m n m d a
It is worth to note that the nonlinear coefficients d n,m depend on only the axial parameter Γ.

For the first family of b.c. we have d n,m =-κ n δ nm (see eq. ( 3)), so that the equations of motion become

+Λ n a n a n +kκ n a n       ∑ ∞ =1 2 i i i a κ =0, (11)
which exhibit only the minor coupling given by the single summation. From [START_REF] Troger | Nonlinear Stability and Bifurcation Theory[END_REF] In dynamical systems language [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields[END_REF], this means that the eigenspaces spanned by the linear modes, which are planar in the linearized case, remain planar also in the nonlinear one, or, equivalently, that the NNMs coincide with the linear ones. For the b.c. of the second family, on the other hand, we have a full nonlinear coupling among modes, and the invariant manifolds no longer remain planar when nonlinearities are considered.

Invariant manifold and homoclinic solutions (nonlinear normal modes technique)

The homoclinic solutions are special cases of the expression (6), w hom (z,t)= , whose time-dependent amplitudes a homoclinic solutions of the hilltop saddle can exist in principle, corresponding to stable and unstable directions of different order. However, we focus attention on only the principal one, namely on the homoclinic solution related to the lowest order, stable and unstable, directions. As we conjectured in the Introduction, they belong to a two dimensional invariant (with respect to the nonlinear equations) manifold Σ which in general is not planar and is difficult to be handled. To overcome this point, the NNM technique, which also allows us to obtain reliable reduced order models, is used [START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF][START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF][START_REF] Nayfeh | Nonlinear Interactions[END_REF]. Although it was originally developed to deal with nonlinear oscillations, herein it is applied to homoclinic orbits.

∑ ∞ =1 hom , ) ( ) (
Three issues are addressed in this section. (i) The presentation of the master-slave technique used to detect the invariant manifold, which includes obtaining the slave equations, (ii) its specification for the homoclinic orbits, and (iii) its application to the two families of b.c. As a matter of fact, the first family can be dealt with without resorting to the NNM technique, so that this section is really needed for only the second family.

Following the treatment of [START_REF] Rand | Lecture Notes on Nonlinear Vibrations[END_REF] (see also the introductory paper in [START_REF]Dimension Reduction of Dynamical Systems: Methods, Models, Applications[END_REF] for an overall framework), we identify a main, or master, variable (i.e., modal amplitude) x and consider the others as secondary, or slave, variables. Because in the buckled beam problem the unstable direction triggering the homoclinic orbit we are looking for is wherever along the first eigenfunction, apart from the range in Fig. 3a where Λ 2 <Λ 1 , it is natural to assume x=a 1 , while a i , i≥2, namely the higher order modal amplitudes, are the slave variables. Note that around the crossover points of Fig. 3a we need two master variables for a consistent analysis, but this is out of the scope of the present work.

The key idea of the method consists in assuming that the slave variables can be expressed as time-independent functions of the master one:

a i =a i (x), i≥2. ( 13
)
If we use only N terms in [START_REF] Hobson | An Efficient Method for Computing Invariant Manifolds of Planar Maps[END_REF] this relation provides an approximation of the searched manifold Σ (approximate manifold), while for N→∞ it exactly gives Σ (exact manifold).

Differentiating [START_REF] Rand | A direct method for non-linear normal modes[END_REF] with respect to time we get, by the chain rule,

= i a dt da i = dx da i dt dx =a i ′ , = a x i a i ′′ x 2 +a i ′ , i≥2. ( 14
) x Combining (14) 2 with (10) yields i a V ∂ ∂ =-a i ′′ x 2 +a i ′ x V ∂ ∂ , i≥2, (15) while
2[H(x,a i )-V(x,a i )]= x 2 [ ] ) ' ( 1 2 2 ∑ ∞ = + j j a . ( 16
)
Eliminating x 2 from ( 15) and ( 16), we obtain

2[V-H]a i ′′+[1+ ] ∑ ∞ =2 2 ) ' ( j j a       ∂ ∂ - ∂ ∂ i i a V x V a '   =0, i≥2, (17)
which are the slave equations (corresponding to the modal equations in Rand [START_REF] Rand | Lecture Notes on Nonlinear Vibrations[END_REF]).

To apply the previous general theory to homoclinic orbits we note that on these solutions it is H=0, so that the strongly nonlinear ODEs permitting to determine the unknown functions (13) read:

2 Va i ′′+[1+ ∑ ' ] ∞ =2 2 ) ( j j a      ∂ ∂ - ∂ i i a V x V a '    ∂ =0, i≥2, (18)
where we have still written a i instead of a i,hom , for the sake of simplicity. Actually, they are much more complex than the original system [START_REF] Rand | Lecture Notes on Nonlinear Vibrations[END_REF], and their introduction may appear useless.

Nonetheless, they contain only polynomial terms and, since we are interested in the finite -though non necessarily small -neighborhood of x=0 where the homoclinic solutions develop, it is natural to look for a relevant solution in the Taylor series form

a i = k 1 [a i,3 (x√k) 3 +a i,5 (x√k) 5 +…]=a i,3 kx 3 +a i,5 k 2 x 5 +…, i≥2. ( 19
)
and consider only a limited number of terms. The larger is this number, the larger is the accuracy of the solution. Note that eq. ( 10) also contains polynomial terms, but in that case Taylor series solutions would require infinite terms because the solution is needed over the whole time line, thus being practically useless. Note also that a solution in the form [START_REF] Nayfeh | On direct methods for constructing nonlinear normal modes of continuous systems[END_REF] resembles the classical machinery of center manifold theory [START_REF] Troger | Nonlinear Stability and Bifurcation Theory[END_REF].

The nonlinearity parameter k has been explicitly singled out so that the Taylor coefficients a i,l are functions of the axial parameter Γ only. The constant terms a i,0 vanish because the manifold passes through the saddle point x=a i =0. The linear terms a i,1 are equal to zero because the manifold is tangent to the linear manifold (which is described by a i =0, i≥2) at the saddle. Even powers a i,2n are not present because eq. ( 1) is symmetric, i.e., both w(z,t) and -w(z,t) are solutions.

To determine the remaining unknown numbers a i,l , we insert [START_REF] Nayfeh | On direct methods for constructing nonlinear normal modes of continuous systems[END_REF] in [START_REF] Pak | Free vibrations of a thin elastica by normal modes[END_REF] and expand in Taylor series. Accounting for (9) we obtain recursive sets of linear algebraic equations yielding the solution.

Equating the coefficients of x 3 to zero we get

a i,3 = i i d d Λ - Λ 1 , 1 1 , 1 9
, i≥2, [START_REF] Lacarbonara | Resonant nonlinear normal modes. Part I: Analytical treatment for onedimensional structural systems[END_REF] from which we see that the no internal resonance conditions 9Λ 1 ≠Λ i , i≥2 must be satisfied to avoid singularities due to interactions between modes. Note that these relations certainly hold if Λ 1 <Λ i , i≥2, as it always occurs, e.g., in the range Γ∈]Γ cr,1 ;Γ cr,2 [ (Fig. 3).

Equating the coefficients of x 5 to zero we get

a i,5 = i Λ - Λ 1 25 1   6 , i≥2, ( 21
)     + + - ∑ ∞ =2 , 1 , 1 , 1 , 1 3 , 3 , 2 1 , 1 ) 2 ( n n i n i n i d d d d a a d
where a i,3 are given by [START_REF] Lacarbonara | Resonant nonlinear normal modes. Part I: Analytical treatment for onedimensional structural systems[END_REF]. Equation ( 21) adds the no internal resonance conditions 25Λ 1 ≠Λ i , i≥2, which again are certainly satisfied if Λ 1 <Λ i , i≥2.

Obtaining higher order terms is much more difficult, and the coefficients a i,l , l=7,9,…, have very involved expressions. We do not pursue the relevant computations, but only infer from ( 20) and ( 21) that further no internal resonance conditions l 2 Λ 1 ≠Λ i , i≥2, are to be satisfied, which certainly

occurs if Λ 1 <Λ i , i≥2.
For the first family of b.c. we have d 1,i =0, i≥2, so that, from (20), a i,3 =0 and then, from (21), a i,5 =0. It is easy to infer that in this case a i,l =0, i≥2, ∀l. This confirms that the nonlinear invariant manifold coincides with the linear one, namely, it is planar. This case corresponds to what Rosenberg [START_REF] Rand | Lecture Notes on Nonlinear Vibrations[END_REF] called similar normal modes, which are solutions of (18) in the linear form a i (x)=c i x+e i , and practically entails solving

c i x V ∂ ∂ = i a V ∂ ∂ , i≥2. (22)
Actually, in the present case, being c i =e i =0, a special case of vanishing similar normal modes occurs. The previous points show a further property of the beam with the first family of b.c.: the relevant nonlinear dynamics at homoclinicity develop along similar normal modes, which in particular coincide with the linear normal modes.

For the second family of b.c. the coefficients a i,3 and a i,5 are not trivial. They are reported in Fig. 4, where four terms are considered in the sum for the computation of a i,5 (see eq. ( 21)). First, we note that the maximum values of both a i,3 and a i,5 for the fixed-free case (Fig. 4c) are O(1), while in the remaining cases they are at least one order of magnitude smaller. From an engineering point of view, this fact and the cubical vanishing of Σ at x=0 mean that the non-flatness of the manifold, while theoretically important, has a minor practical relevance, apart from the fixed-free case. This holds in a wide range of Γ values, and confirms a preliminary observation in [START_REF] Lenci | Optimal control of the homoclinic bifurcation in buckled beams: infinite-dimensional vs reduced order modeling[END_REF], which has shown Figure 4. The first four coefficients a i,3 and a i,5 for a) fixed-fixed (Γ cr,1 =4π 2 , Γ cr,2 =8.1830π 2 ; here a 2,3 =a 4,3 =a 2,5 =a 4,5 =0 because odd eigenfunctions satisfy w 2n-1 (1/2+z)=w 2n-1 (1/2-z) while even eigenfunctions satisfy w 2n (1/2+z)=-w 2n (1/2-z), so that d n,m =0 if n is odd and m is even), b) fixedhinged (Γ cr,1 =2.0457π 2 , Γ cr,2 = 6.0468π 2 ), c) fixed-free (Γ cr,1 =π 2 /4, Γ cr,2 =9π 2 /4) and d) fixed-guided (Γ cr,1 =π 2 , Γ cr,2 =4π 2 ) beams. The coefficients a i,5 are computed with 4 terms in the sum in [START_REF] Srinil | Two-to-one resonant multi-modal dynamics of horizontal/inclined cables. Part II: Internal resonance activation, reduced order models and nonlinear normal modes[END_REF].

A c c e p t e d m a n u s c r i p t

The second observation following from Fig. 4 is that, on average, a i,5 is one order of magnitude smaller than a i,3 , showing how the higher order nonlinearities (Taylor coefficients a i,l ) in the slave modal amplitudes are important only for some isolated cases. Furthermore, still on average, both a i,3

and a i,5 are decreasing functions of i in absolute value. This means that high order slave amplitudes are practically unessential. A remarkable exception is, again, the fixed-free b.c. (Fig. 4c).

Finally, we note that a i,3 and a i,5 can have different signs for fixed values of Γ, showing how the non-flat manifold Σ can have complicated shapes in phase space.

Another simplification occurs for symmetric b.c., as the fixed-fixed ones of Fig. 2a. In this case, the odd eigenfunctions are symmetric with respect to the point z=1/2, namely, w 2n-1 (1/2+z) = w 2n-1 (1/2-z), while the even eigenfunctions are anti-symmetric, w 2n (1/2+z)=-w 2n (1/2-z), so that they are mutually orthogonal and d n,m =0 if n is odd and m even (or viceversa). In other words, odd and even eigenfunctions span two orthogonal subspaces.

Reduced order models (ROMs)

By inserting the a i (x) obtained in the previous section in the expression (9) 2 of the potential V(x,a i (x)) we get, after some computations,

V(x)= k 2 1 {α 2 (
α 2 =Λ 1 <0, α 4 = 2 2 1 , 1 d >0, α 6 = = ∑ ∞ = + Λ 2 3 , 1 1 , 1 2 3 ] 2 ) ( [ n n n n n a d d a 2 1 , 1 d ∑ ∞ = Λ - Λ Λ - Λ 2 2 1 1 2 , 1 ) 9 ( 18 n n n n d , α 8 = , ( 23
) ∑ ∞ = + Λ - Λ 2 2 3 , 2 1 , 1 5 , 3 , 1 ] ) ( 6 ) 43 [( n n n n n a d a a
where again the parameter k has been singled out, and where a n,3 and a n,5 are given by ( 20) and [START_REF] Srinil | Two-to-one resonant multi-modal dynamics of horizontal/inclined cables. Part II: Internal resonance activation, reduced order models and nonlinear normal modes[END_REF].

In what follows we assume that Γ∈]Γ cr,1 ;Γ cr,2 [, which is the region of main practical interest. This implies that Λ 1 <0 and Λ n >0, n≥2, and thus easily proves that α 6 is negative. The sign of α 8 , on the other hand, is not known in advance (see Fig. 5

below).

The coefficients α i are at the base of the reduced order models (single d.o.f., indeed) obtained in this section, and α i , i=6,8,…, summarize the effects of non-flatness of the manifold Σ. They are The equation of motion is

α 2 =π 2 (π 2 -Γ), α 4 = 2 4 π , α i =0, i=6,8,…, (hinged-hinged beam), α 2 =(π/2) 2 [(π/2) 2 -Γ], α 4 = 32 4 π , α i =0, i=6,8,…, (guided-hinged beam), (24)
+ x dx dV = +α x 2 x+2α 4 kx 3 +3α 6 k 2 x 5 +4α 8 k 3 x 7 +…=0, (25)

A c c e p t e d m a n u s c r i p t

while the homoclinic orbit of the hilltop saddle, whose determination is the main object of this work, can be obtained by solving

dt dx = ) ( 2 x V - = ... 8 2 8 6 2 6 4 4 2 2 - - - - - x k x k kx x α α α α . ( 26
)
It is useful to perform the following change of variables

y=         -2 4 α α k x 2 , τ=2t 2 α -, (27)
so that ( 26) becomes

τ d dy = ... 1 3 3 2 2 + - + - y y y y β β , β 2 = 2 4 6 2 ) (α α α , β 3 = 3 4 8 2 2 ) (α α α , …. ( 28
)
which is k-independent. Note also that with this change of variables the two symmetric homoclinic loops coincide, so that in terms of y we have only one homoclinic solution.

Integrating [START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF] we obtain the expression τ=τ(y) of the homoclinic loop

τ=-∫ + - + - y y d ... 1 3 3 2 2 ξ β ξ β ξ ξ ξ , ( 29
)
where y =y hom (0) is the smallest positive solution of 1-y+β 2 y 2 -β 3 y 3 +….=0. It represents the maximum extension of the homoclinic loop, which in fact belongs to the interval y∈[0, y ]. Finally, y hom (τ) can be obtained by inverting [START_REF] Nayfeh | Nonlinear Interactions[END_REF].

Equations ( 23) and ( 25) provide a hierarchy of approximate reduced order, single d.o.f., models (ROMs), obtained by accounting for an increasing number l of higher order nonlinearities a i,l in the ith slave amplitude.

Unrefined ROM

If we neglect α 6 , α 8 ,…, i.e., consider a quartic potential, we obtain

+α x 2 x+2α 4 kx 3 =0, (30)
which is the well-known Duffing equation. It corresponds to the governing equation obtained by the classical (linear) Galerkin projection of the dynamics onto the first linear mode. In fact, it is easy to verify that if we assume w(z,t)=x(t)w 1 (z), insert this expression in (1), multiply by w 1 (z) and integrate from z=0 to z=1, we get [START_REF] Lenci | Optimal control of the homoclinic bifurcation in buckled beams: infinite-dimensional vs reduced order modeling[END_REF]. It is worth to remark that (30) provides an exact result for the first family of b.c., for which α i =0, i=6,8,…, whereas just an unrefined approximation for the second family. In the latter case, it can be considered as a minimal approximation, in the sense that

A c c e p t e d m a n u s c r i p t

it is the simplest ROM permitting to approximate the homoclinic orbits of interest.

The homoclinic loop of eq. ( 30) can be computed in closed form, by integrating and then inverting [START_REF] Nayfeh | Nonlinear Interactions[END_REF]. This yields

y 3,hom (τ)= ) cosh( 1 2 τ + , (31)
which, in these variables, is independent of k and Γ, namely of the b.c..

First refined ROM

For the second family of b.c. it is possible to improve the approximation by considering α 6 and neglecting higher order terms. The equation of motion becomes

+α x 2 x+2α 4 kx 3 +3α 6 k 2 x 5 =0. ( 32
)
To proceed with the computation of the homoclinic solution, we need y , which at the considered order is given by

y = 2 2 2 4 1 1 β β - - . (33)
We see that y exists if and only if β 2 <1/4. This suggests that this model must be considered with care. In fact, from [START_REF] Pesheck | Accurate reduced-order models for a simple rotor blade model using nonlinear normal modes[END_REF] we know that α 6 is always negative, see Fig. 5 (and then β 2 is positive), so that, strictly speaking, this approximation corresponds to a softening oscillator, which is in contrast with the original mechanical system, which is certainly hardening.

However, for α 6 sufficiently small (and thus for β 2 small), the softening "error" manifests itself for only large values of y, and far from the region of practical interest, which is that occupied by the homoclinic loop (see Fig. 12a below). In this case, y and the homoclinic solution certainly exist.

When |α 6 | (and β 2 ) increase, the escape region approaches the homoclinic orbit, and the approximation ( 32) is initially no longer reliable, and then fails to possess a homoclinic orbit (when point B in Fig. 12a collapses onto point A on the zero potential line). This occurs when β 2 =1/4. To illustrate how this threshold depends on the b.c., we have reported in Fig. 6 Above the thresholds β=0 highlighted in Fig. 6, the ROM (32) is meaningless, and one must retain higher order terms. To reproduce the hardening behaviour at infinity, one must have a positive highest order coefficient of the considered approximation. In several cases, but not always,

A c c e p t e d m a n u s c r i p t

this can be obtained by considering terms up to α 8 (see Fig. 5). Below the threshold, on the other hand, the present approximation can be practically used, although with a varying degree of reliability. The homoclinic loop of the saddle x=0 can again be computed in closed form, by integrating and then inverting (29):

y 5,hom (τ)= ) cosh( 1 2 τ β + . ( 34
)
The possibility to have an analytical expression for the homoclinic orbit makes model [START_REF] Fung | Buckling of Low Arches or Curved Beams of Small Curvature[END_REF] interesting in any case. By comparing (34) with [START_REF] Holmes | A Partial Differential Equation with Infinitely many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF] we can appreciate the effects of the nonplanarity of the manifold (up to its lowest order) on the homoclinic loop. Fig. 7 reports the time histories and the phase portraits, holding for all b.c., from which we note a marked difference, even for β not so far from 1, i.e., even for Γ not so larger than Γ cr,1 (Fig. 6). (36)

This model reproduces the hardening behaviour of the original system if α 8 >0, which is therefore assumed in the following (Fig. 5).

To determine y we have to solve 0=1-y+β 2 y 2 -β 3 y 3 . Depending on the values of β 2 and β 3 , this equation can have one or three real solutions. The regions of different behaviour are depicted in white and gray in Fig. 10, respectively. The former case corresponds to a two-well potential, while the latter corresponds to a four-well potential (Fig. 10), a situation which is not expected in practice.

In fact, it is also shown in Fig. 10 that the actual curve (β 2 (Γ),β 3 (Γ)) for the fixed-fixed b.c. does not cross the grey area, so that, as expected, in the real case there is only one solution y , which is depicted in Fig. 11 as a function of Γ. When compared with y (Γ)=1 of the unrefined ROM (see Sect. 5.1), this figure shows again that the non-flatness of the manifold Σ may have non-negligible effects. Figure 11 (as well as Fig. 7b) shows that the linear Galerkin method produces the maximum underestimation of y , which corresponds to the stiffest approximation of the real system ensuing from projection of its global dynamics onto the first mode (master amplitude) with no contributions from the higher order modes (slave amplitudes).

Similar situations occur with the other b.c. of the second family, at least for Γ far away from Γ cr,2 , which are not reported to limit the length of the paper. To extend the previous comparison to other values of the axial parameter, we have depicted in Fig. 13 the maximum extension x of the homoclinic orbit as a function of Γ for the three models.

Note that 5

x exists only up to (Γ-Γ cr,1 )/(Γ cr,2 -Γ cr,1 )≅0.76 (see Fig. 6a), and in the neighborhood of this limit value the first refined ROM is not accurate. Figure 13 shows that just above the static buckling load Γ cr,1 the effects of the non-planarity of the manifold Σ are negligible, while they become more and more important as the axial parameter increases.

) ( } )] ( [ )] ( [ { n n n n z w t x k a t x k a
In the unrefined model (which is nonetheless exact for the first family of b.c.), the un-corrected master amplitude and time-independent spatial shape do occur. In contrast, in the refined models, the approximation of the master modal amplitude is improved by accounting for (a theoretically infinite number of) slave modal amplitudes a n,l at different l-orders. Of course, for practical purposes, the n-summations in (37) can be limited to only few terms. Yet, irrespective of the considered number of terms, all of the refined models account, at different l-orders, for a twofold effect on the homoclinic orbit of the higher order modes described by the slave amplitudes. (i) They affect the solely time-dependent master amplitude of the first mode, which enters the homoclinic orbit expression at increasing power orders in the various refined models. (ii) They entail a time varying spatial shape of the homoclinic orbit ensuing from a nonlinear Galerkin projection of the solution onto the basis of linear normal modes. In both respects, apart from referring to homoclinic orbits and from being obtained based on a priori Galerkin projection of the dynamics in the spirit of the original (see the Introduction) NNM approach [START_REF] Rosenberg | The normal modes of nonlinear n-degrees-of-freedom systems[END_REF], the structure of the obtained solution closely resembles the refined structure of nonlinear displacement fields as provided a posteriori for a number of infinite-dimensional systems (see, e.g., [START_REF] Lacarbonara | Resonant nonlinear normal modes. Part I: Analytical treatment for onedimensional structural systems[END_REF]) by the alternative NNM approach [START_REF] Pak | Free vibrations of a thin elastica by normal modes[END_REF][START_REF] Nayfeh | On direct methods for constructing nonlinear normal modes of continuous systems[END_REF].

Finally, it would be worth: (i) To perform a convergence analysis aimed at evaluating the number of slave amplitudes to be considered for obtaining a reliable solution; (ii) to evaluate the relative importance of considering either one more higher-order nonlinearity within the chosen number of slave amplitude(s) or one more slave amplitude up to an overall lower-order level of nonlinearity; (iii) to compare the spatial shapes of the solution as provided by the various ROMs.

All of these issues of main practical interest are left for possible future research.

Conclusions and further developments

The homoclinic solutions of an initially straight, planar buckled beam have been investigated in an infinite-dimensional context. Upon conjecturing that the homoclinic solutions exist and belong to a two-dimensional invariant manifold Σ, we have considered all possible boundary conditions and shown how they can be grouped in two families. However, for the purposes of the present work, the most important properties are that (i) the invariant manifold Σ is flat (and coincides with the manifold spanned by the first eigenfunction), and that (ii) the homoclinic solutions can be easily calculated in analytical form. They coincide with those computed via the (linear) Galerkin projection onto the solely first mode.

These properties are no longer valid for the second family, which contains all of the remaining b.c.. In particular, no exact expressions of the homoclinic solutions are directly found for these b.c., so that we have attacked the problem by means of the NNM technique [START_REF] Rand | Lecture Notes on Nonlinear Vibrations[END_REF]. This provides series expressions of the manifold Σ and of the homoclinic orbits, and requires to solve an infinite sequence of linear algebraic problems. Only few terms are indeed sufficient for practical purposes, so that only approximate manifolds and orbits are considered.

As a matter of fact, there are two kinds of approximations. The first is concerned with the number of linear normal modes used to span the space of admissible functions. The second with the order of the highest nonlinear term considered in the Taylor expansion of the functions describing the non-planar manifold. As to the former, we have identified the amplitude of the first linear eigenfunction as the solely independent variable, with all of the other modes being considered as dependent variables. Thus, we have considered single d.o.f. models. As to the latter, we have obtained a hierarchy of reduced order models based on successive approximations of the potential.

The simplest one, corresponding to a quartic potential, is the classical Duffing oscillator obtained by the linear Galerkin projection. A first refined model is obtained by considering up to the sixth term in the potential, a second refined model up to the eighth term, and so on. The first three ROMs have been analyzed in detail and compared with each other.

It has been shown that, while being always important from a theoretical point of view, the nonplanarity of the manifold where the homoclinic solutions lie becomes increasingly important also from a practical point of view when the axial parameter grows much above the buckling value.

The analysis in this paper is a necessary step in the involved process of detecting chaotic dynamics in infinite dimensional systems. The next step consists of evaluating how damping and excitation modify the homoclinic orbits. Indeed, it is known that, depending on the relative amplitudes of these parameters, homoclinic orbits may survive or disappear, and that the boundary between the two different regimes is given by the homoclinic bifurcation. This global event can be detected by means of the Holmes and Marsden [START_REF] Holmes | A Partial Differential Equation with Infinitely many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF] theorem, which is an infinite-dimensional extension of the classical Melnikov method [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields[END_REF]. For fixed-fixed b.c. and for a given value of Γ close to Γ cr,1 , this has been done in [START_REF] Lenci | Optimal control of the homoclinic bifurcation in buckled beams: infinite-dimensional vs reduced order modeling[END_REF]. The results of the present work strongly call for a generalization

Figure 1 .A c c e p t e d m a n u s c r i p tFigure 2 .

 12 Figure 1. The first family of b.c.: a) hinged-hinged and b) guided-hinged.

Figure 3 .

 3 Figure 3. The first four eigenvalues for a) fixed-fixed, b) fixed-hinged, c) fixed-free and d) fixedguided beams.

  governing equations easily follow:

  ,hom (t) correspond to vanishing of the Hamiltonian (9) 1 . Many A c c e p t e d m a n u s c r i p t

A c c e p t e d m a n u s c r i p t

  [START_REF]Dimension Reduction of Dynamical Systems: Methods, Models, Applications[END_REF] 1 and (14) 1 give

A c c e p t e d m a n u s c r i p t

A c c e p t e d m a n u s c r i p tFigure 5 .

 5 Figure 5. The coefficients α 2 , α 4 , α 6 and α 8 for the second family of b.c.. The coefficients α 6 and α 8 are computed with 4 terms in the sums in (23). The abscissa ranges from Γ cr,1 to Γ cr,2 .

  33), which is central for this model (note that β=1 corresponds to the unrefined ROM). It is always imaginary for the case of fixed-free b.c., which again shows a different behaviour with respect to the other b.c. of the second family.

  Figure 6. The parameter β=

Figure 7 .A c c e p t e d m a n u s c r i p t 5 . 3 .

 753 Figure 7. Homoclinic loops y 3,hom (τ) (thick line, β=1) and y 5,hom (τ) (thin lines, β=0.75 and β=0.5). a) time histories and b) phase portraits.

  e d -f i x e d b . c .

Figure 10 .

 10 Figure 10. Behaviour chart for the solution of 0=1-y+β 2 y 2 -β 3 y 3 , and parametric curve (β 2 (Γ),β 3 (Γ)) for fixed-fixed b.c.

3 Figure 12

 312 Figure 12. a) The potentials kV 4 (x√k), kV 6 (x√k) and kV 8 (x√k) and b) the phase portrait of the homoclinic loops x 3,hom (t), x 5,hom (t) and x 7,hom (t). Fixed-fixed b.c., (Γ-Γ cr,1 )/(Γ cr,2 -Γ cr,1 )≅0.64.

3 Figure 13 .

 313 Figure 13. The maximum extension ( x √k) of the homoclinic orbit for the three ROMs. Fixedfixed b.c..

The

  first family includes the hinged-hinged and guided-hinged b.c.. Various relevant properties have been illustrated all along the paper. (i) An exact solution of the system PDE can be pursued in the form of separate variables; (ii) the eigenvalues of the linearized PDE depend linearly on the governing parameter Γ representing the axial load; (iii) the associated projected (onto the linear normal modes) equations of motion have a very weak coupling; and (iv) the relevant nonlinear A c c e p t e d m a n u s c r i p t dynamics at homoclinicity develop along similar normal modes (as defined by Rosenberg), which in particular coincide with the linear ones.

A c c e p t e d m a n u s c r i p t

  of the homoclinic bifurcation analysis to other values of the axial parameter and to different b.c.

Table 1 .

 1 The boundary conditions of the first and of the second family.

		z=0			z =1
	w w'''+Γw' w' w'' name w w'''+Γw' w' w'' name
	x		x hinged x		x hinged first family
	x		x hinged	x	x	guided first family
	x	x	fixed x		x	fixed second family
	x	x	fixed x		x hinged second family
	x	x	fixed	x	x	guided second family
	x	x	fixed	x	x	free second family

  x√k) 2 +α 4 (x√k) 4 +α 6 (x√k) 6 +α 8 (x√k) 8 +…}= 2 1 {α 2 x 2 +α 4 kx 4 +α 6 k 2 x 6 +α 8 k 3 x 8 +…},

2n

Γ Λ n Λ 1 Λ 2 Λ 3 Λ 4 Γ cr,1 Γ cr,2 Γ cr,3 Γ cr,4-80008000 350 0 Γ Λ n Λ 1 Λ 2 Λ 3 Λ 4 Γ cr,1 Γ cr,2 Γ cr,3 Γ cr,4 a) b) -4000 4000 250 0 Γ Λ n Λ 1 Λ 2 Λ 3 Λ 4 Γ cr,1 Γ cr,2 Γ cr,3 Γ cr,4-4000 4000 250 0Γ Λ 1 Λ 2 Λ 3 Λ 4 Γ cr,1Γ cr,2 Γ cr,3 Γ cr,4Λ n

this property to hold for a fixed parameter value and for fixed-fixed b.c.. a) 0 0.03 Γ a 3,3 Γ cr,1 Γ cr,2 -0.001 0 Γ a 3,5 Γ cr,1 Γ cr,2 b) 0 0.08 Γ a 2,3 Γ cr,1 Γ cr,2 a 3,3 a 4,3 0 0.005 Γ a 2,5 Γ cr,1 Γ cr,2 a 3,5 a 4,5 c) -0.3 0.5 Γ a 2,3 Γ cr,1 Γ cr,2 a 3,3 a 4,3 -0.15 0.35 Γ a 2,5 Γ cr,1 Γ cr,2 a 3,5 a 4,5 d) -0.01 0.15 Γ a 2,3 Γ cr,1 Γ cr,2 a 3,3 a 4,3 -0.003 0.02

 1) cosh( 1

which ranges from δ(0)=2/(1+β) to δ(∞)=1/β and is depicted in Fig. 8. The differences are important from an engineering point of view, showing that the non-planarity of Σ cannot be neglected, even for practical purposes. Note that the asymptotically attained difference is greater than the initial one.

To better visualize the effects of the non-planarity, we have drawn in a 3D section of the phase space (Fig. 9) the invariant manifold Σ with the homoclinic loop y 5,hom (τ), its projection on the plane (y, ), and its approximation y y 3,hom (τ), which, as said, is obtained by the classical Galerkin projection of the original infinite dimensional system onto the first mode. Note that the vertical direction is blown-up to highlight the non-flatness of the manifold. Once y is known, we can compute the integral (29) in closed form (although with a complex formula), but unfortunately we are not able to invert the expression τ=τ(y); thus we have no analytical expression of y 7,hom (τ), although it can be easily obtained numerically. However, the effects of the manifold curvature on the homoclinic solution evidenced in Fig. 11 are sufficient to remark its importance, so that we do not pursue the computation of τ=τ(y) further on.

Comparisons between ROMs

In this section we compare the results obtained with the three ROMs previously introduced. We consider only the fixed-fixed case. Furthermore, as we want to perform a comparison of the physical potentials V(x), we use the true master amplitude x instead of its abstract version y, however accounting for the relevant symmetry.

The potentials V 4 (x), V 6 (x) and V 8 (x) of the three ROMs are reported in Fig. 12a, which shows the corrections due to the higher order nonlinearities. The corresponding phase portraits are reported in Fig. 12b. Apart from the quantitative differences, which are of the order of 10÷15%, we see that there is a non-monotonic convergence toward the limit potential V ∞ representing the exact homoclinic solution of the buckled beam.

The softening behaviour of the potential V 6 is clearly visible in Fig. 12a, because the considered value of Γ (corresponding to β=0.458) is not so far from the value (Γ-Γ cr,1 )/(Γ cr,2 -Γ cr,1 )≅0.76 corresponding to the critical threshold β=0 (see Fig. 6a). For smaller values of Γ, the right hilltop is higher and out of the chosen window; when Γ tends to (Γ-Γ cr,1 )/(Γ cr,2 -Γ cr,1 )≅0.76, on the other hand, point B collapses onto point A on the zero potential line.

A c c e p t e d m a n u s c r i p t