Julien Favier
email: julien.favier@univ-amu.fr

Alistair Revell

Alfredo Pinelli

Numerical study of flapping filaments in a uniform fluid flow

Keywords: beating filaments, flapping flags, inverted hydrodynamic drafting, Immersed Boundary, Lattice Boltzmann

The coupled dynamics of multiple flexible filaments (also called monodimensional flags) flapping in a uniform fluid flow is studied numerically for the cases of a side-by-side arrangement, and an in-line configuration. The modal behaviour and hydrodynamical properties of the sets of filaments are studied using a Lattice Boltzmann -Immersed Boundary method. The fluid momentum equations are solved on a Cartesian uniform lattice while the beating filaments are tracked through a series of markers, whose dynamics are func-

Introduction

The scope of this work is the physical analysis of the dynamics of flapping filaments in a streaming ambient fluid, which has a large spectrum of applications in aeronautics, civil engineering or biological flows. From the theoretical side, this fluid structure interaction problem is particularly challenging as it involves non-linear effects as well as large structural deformations [START_REF] Païdoussis | Fluid-structure interactions: slender structures and axial flow[END_REF][START_REF] Shelley | Flapping and bending bodies interacting with fluid flows[END_REF]. The present study is particularly inspired by various experiments on flapping filaments realised in soap films [START_REF] Zhang | Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind[END_REF][START_REF] Zhu | Interaction of two flapping filaments in a flowing soap film[END_REF][START_REF] Ristroph | Anomalous Hydrodynamic Drafting of Interacting Flapping Flags[END_REF]. Indeed, soap film experiments associated to thin-film interferometry for flow visualisation can be considered as a reasonable approximations of 2D fluid structure interaction scenarios, thus suitable for the validation of the results obtained with our 2D numerical approach.

In our simulations, we consider a 2D incoming incompressible flow modeled using a Lattice Boltzmann method, coupled to a model of infinitely thin and inextensible filament experiencing tension, gravity, fluid forces and flexural rigidity (i.e. a bending term in the form of a 4th derivative with respect to the curvilinear coordinate describing the filament). Also, at all time instants tension forces are determined to maintain the inextensibility of the structure. In this simple model the energy balance of the system is driven by the bending forces and fluid forces, as the structure is controlled by an inextensibility constraint which prohibits stretching or elongation motions that would dissipate energy. This system encompasses all the essential ingredients of a complex fluid-structure interaction problem: large deformations, slender flexible body, competition between bending versus fluid forces, inextensibility and effect of the filament tips on the surrounding flow as vorticity generators.

To enforce the presence of the solid on the fluid lattice, we use a variant of the immersed boundary method previously developed by the authors [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume navier-stokes solvers[END_REF] on finite difference and finite volume Navier Stokes solvers. In this work we use the same algorithm as in [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume navier-stokes solvers[END_REF] to impose the immersed boundary forces, but we adapt it to a Lattice Boltzmann solver. This approach for imposing the forces has shown to be order 2 in space, computationally cheap and directly provides for the forces exerted on the fluid by the filaments without the introduction of any empirical parameter. Using the Lattice Boltzmann method in conjunction with an Immersed Boundary technique to solve the motion of an incompressible fluid also allows for a clean imposition of the boundary conditions on the solid since it does not suffer from errors originating from the projection step, as it is the case when associated with unsteady incompressible Navier Stokes solvers [START_REF] Domenichini | On the consistency of the direct forcing method in the fractional step solution of the navier-stokes equations[END_REF].

Making use of the outlined Lattice Boltzmann -Immersed Boundary approach, we consider the coupled dynamics of systems made of highly deformable flexible filaments, as introduced by [START_REF] Favier | A lattice boltzmann -immersed boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF]. No artificial contact force is introduced between the filaments, in order to preserve a purely hydrodynamical interaction. We focus in this work on the modal behaviour of a set of two and three side-by-side filaments, by varying the spacing between them. The obtained results confirm the theoretical predictions and experimental observations mentioned in literature. the so-called anomalous hydrodynamic drafting pointed out experimentally in [START_REF] Ristroph | Anomalous Hydrodynamic Drafting of Interacting Flapping Flags[END_REF] is recovered here numerically and a physical mechanism is proposed to explain this phenomenon.

Coupled Lattice Boltzmann -Immersed Boundary Method

This fluid-structure problem is tackled using an Immersed Boundary method coupled with a Lattice Boltzmann solver. In the following we provide a summary of the numerical technique while details of the methodology can be found in [START_REF] Favier | A lattice boltzmann -immersed boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF].

The fluid flow is modeled by advancing in time the Lattice Boltzmann equation which governs the transport of particles density distribution f (probability of finding a particle in a certain location with a certain velocity). It is often classified as a mesoscopic method, where the macroscopic variables, namely mass and momentum, are derived from the distribution functions f . An excellent review of the method can be found in [START_REF] Succi | The lattice Boltzmann equation[END_REF].

Using the classical BGK approach [START_REF] Bhatnagar | A model for collision processes in gases. i: small amplitude processes in charged and neutral one-component system[END_REF], and after an appropriate discretization process [START_REF] Malaspinas | Lattice boltzmann method for the simulation of viscoelastic fluid flows[END_REF], the Boltzmann transport equation for the distribution function f = f (x, e, t) at a node x and at time t with particle velocity vector e is given as follows:

f i (x + e i ∆t, t + ∆t) -f i (x, t) = - ∆t τ f (x, t) -f (eq) (x, t) + ∆tF i (1)
In this formulation, x are the space coordinates, e i is the particle velocity in the i th direction of the lattice and F i accounts for the body force applied to the fluid, which conveys the information between the fluid and the flexible structure. The local particle distributions relax towards an equilibrium state f (eq) in a single time scale τ . Equation 1 governs the collision of particles relaxing toward equilibrium (first term of the r.h.s.) together with their streaming which drives the data shifting between lattice cells (l.h.s of the equation). The rate of approach to equilibrium is controlled by the relaxation time τ , which is related to the kinematic viscosity of the fluid by ν = (τ -1/2)/3. Equation 1 is approximated on a Cartesian uniform grid by assigning to each cell of the lattice a finite number of discrete velocity vectors. In particular, we use the D2Q9 model, which refers to two-dimensional and nine discrete velocities per lattice node (corresponding to the directions east, west, north, south, center, and the 4 diagonal directions as given by equation 2), where the subscript i refers to these discrete particle directions. As is usually done, a convenient normalization is employed so that the spatial and temporal discretization in the lattice are set to unity, and thus the discrete velocities are defined as follows:

e i = c 0 1 -1 0 0 1 -1 1 -1 0 0 0 1 -1 1 -1 -1 1 (i = 0, 1, ..., 8) (2)
where c is the lattice speed which defined by c = ∆x/∆t = 1 with the current normalization. The equilibrium function f (eq) (x, t) can be obtained by Hermite series expansion of the Maxwell-Boltzmann equilibrium distribution [START_REF] Qian | Lattice bgk models for navier-stokes equation[END_REF]:

f (eq) i = ρω i 1 + e i • u c 2 s + (e i • u) 2 2c 4 s - u 2 2c 2 s (3)
In equation 3, c s is the speed of sound c s = 1/ √ 3 and the weight coefficient ω i are ω 0 = 4/9, ω i = 1/9, i = 1...4 and ω 5 = 1/36, i = 5...8 according to the current normalization. The macroscopic velocity u in equation 3 must satisfy the requirement for low Mach number, M, i.e. that | u | /c s ≈ M << 1. This stands as the equivalent of the CFL number for classical Navier Stokes solvers. The force F i in equation 1 is computed using a power series in the particle velocity with coefficients that depend on the actual volume force f ib applied on the fluid. The latter is determined using the Immersed Boundary method, originally introduced by [START_REF] Peskin | The immersed boundary method[END_REF], following the formulation described in [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] for classical Navier Stokes solvers [START_REF] Favier | A lattice boltzmann -immersed boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF] for Lattice Boltzmann solvers. In this approach, the flexible filaments are discretised by a set of markers X k , that in general do not correspond with the lattice nodes x i,j . The role of f ib is to restore the desired velocity boundary values on the immersed surfaces at each time step.

The global algorithm is decomposed as follows. The Lattice-Boltzmann equations for the fluid are first advanced to the next time step without immersed object (F i = 0), which provides the distribution functions f i needed to build a predictive velocity u p by ρu p = i e i f i and ρ = i f i . The predictive velocity is then interpolated onto the structure markers, which allows one to derive the forcing required to impose the desired boundary condition at each marker using:

F ib (X k ) = U d n+1 (X k ) -I[u p ](X k ) ∆t (4)
In equation 4, capital letters are used to identify variables defined on each marker and I[u p ](X k ) refers to the interpolated predictive velocity. The term U d n+1 (X k ) denotes the velocity value at the location X k we wish to obtain at time step completion. Adding this force term to the right hand side of the momentum equations allows to restore the desired velocity U d n+1 at the boundary (see [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] for instance). The value of U d n+1 is determined for each filament and at each Lagragian point by integrating in time the respective dynamic equation:

dU d n+1 dt = ∂ ∂s (T ∂X k ∂s ) -K B ∂ 4 X k ∂s 4 + Ri g g -F ib (5)
Here, the Richardson number is Ri = gL/U 2 ∞ , T is the tension of the filament and K B is the flexural rigidity. All variables are non dimensional and the reference quantities used for the normalisation are: the reference force tension

T ref = ∆ρU 2 ∞ , the reference bending rigidity K B ref = ∆ρU 2 ∞ L 2 and the reference Lagrangian forcing F ref = ∆ρ ερ f U 2 ∞ ,
where ε refers to the numerical thickness of the filament on the lattice. More details can be found in [START_REF] Favier | A lattice boltzmann -immersed boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF]. The closure of equation 5 is provided by the inextensibility condition that reads:

∂X k ∂s • ∂X k ∂s = 1 (6)
This condition, that ensures that the filament does not stretch (and thus its length remains constant), is satisfied using the tension values that effectively act as Lagrange multipliers. The boundary conditions for the system (5-6) are X = X 0 , ∂ 2 X k ∂s 2 = 0 for the fixed end and T = 0, ∂ 2 X k ∂s 2 = 0 for the free end.

Returning to equation 4, the term I[u p ](X k ) refers to the value of the predictive velocity field interpolated at X k . This provides the kinematic compatibility between solid and fluid motion, i.e. zero relative velocity on the solid boundary. At this stage, the required forcing is known at each marker by equation 4, and needs to be spread onto the lattice neighbours by:

f ib (x) = S(F ib (X k )).
More details on the interpolation operator I, spreading operator S and the filaments equations of motion can be found in [START_REF] Favier | A lattice boltzmann -immersed boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF]. With respect to literature, this approach of immersed boundary preserves an order 2 in space and ensures a conservation of the force and torque between Eulerian and Lagrangian space, for a relatively low computational cost. The forcing f ib is finally discretised on the lattice directions and reads [START_REF] Malaspinas | Lattice boltzmann method for the simulation of viscoelastic fluid flows[END_REF][START_REF] Guo | Discrete lattice effects on the forcing term in the lattice boltzmann method[END_REF]:

F i = 1 - 1 2τ ω i e i -u c 2 s + e i • u c 4 s e i • f ib (7)
Once the system of hydrodynamic forces has been determined following the outlined procedure, equation 1 is then solved once again with the forcing F i which impose the correct boundary condition at each marker X k . Finally, the macroscopic quantities are then derived from the obtained distribution functions f by ρu = i e i f i + ρ ∆t 2 F and ρ = i f i , which closes one time step of the solver.

One single flapping filament in an incoming fluid flow

Following the experiments of [START_REF] Zhu | Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method[END_REF], and the numerical study of [START_REF] Huang | Simulation of flexible filaments in a uniform flow by the immersed boundary method[END_REF], we start by considering the beating of a single filament fixed at one end, and subject to gravity and hydrodynamics forces. Let L be the length of the filament, we fix the density difference between solid and fluid ∆ρ = ρ s -ρ f L to ∆ρ = 1.5, the non-dimensional bending rigidity to K B = 0.001, and the value of the Richardson number to Ri = 0.5.

The inlet velocity imposed in the Lattice-Boltzmann normalization is set to U ∞ = 0.04 (aligned with gravity direction), with a relaxation time of τ = 0.524 and a filament length of L = 40. With these values, the simulation is run at a Reynolds number Re = U ∞ L/ν equal to 200. The size of the computational domain is set to 15L × 10L, in the streamwise and transverse direction respectively. The lattice discretization (600 × 400 nodes) has been determined as the result of a preliminary grid convergence study. The initial angle of the filament is set to θ = 18 o with respect to the gravity direction, and its fixed end is placed at the centerline of the domain, at a distance of 4L from the inlet. The L2 norm of the inextensibility error is kept below 10 -12 at all times.

(a) Periodic flapping pattern.

(b) y-coordinate of the free end. Figure 1a shows the periodic pattern of the beating in the established regime, characterised by sinuous traveling waves moving and amplifying downstream from the fixed end. The same behavior has been observed both in the simulations of [START_REF] Bagheri | Spontaneous symmetry breaking of a hinged flapping filament generates lift[END_REF] and in the experiments of [START_REF] Shelley | Flapping and bending bodies interacting with fluid flows[END_REF]. Figure 1b displays the time evolution of the y-coordinate (transverse direction) of the free end of the filament. After six beating cycles, a periodic orbit is established, with a period of 3 time units (the same value as the one found by [START_REF] Huang | Simulation of flexible filaments in a uniform flow by the immersed boundary method[END_REF]). The predicted amplitude of the beating compares well: the difference with reference data on the maximal excursion of the free end is less than 5%. Also, the peculiar trajectory of the free end exhibiting a characteristic figure-eight orbit (dashed line in 1a) is recovered, in agreement with the findings of the soap film experiments carried out by Zhang et al. [START_REF] Zhang | Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind[END_REF].

Figures 2a and2b show the effect of the bending rigidity coefficient on the beating pattern. Without bending rigidity (figure 2a), the filament is totally flexible and a rolling up of the free extremity is observed; this effect has been termed as kick following [START_REF] Bailey | Motion of a hanging chain after the free end is given an initial velocity[END_REF]. On the other hand, when the filament has a finite flexural rigidity (K B = 0.001 in this simulation), the rolling up of the free end is inhibited, the kick disappears and the flapping amplitude is reduced. Thus, the proposed slender structure model, incorporating both bending terms and tension, computed to enforce inextensibility, reproduces successfully the same phenomena as the ones observed in experiments.

Side-by-side flapping flaps

We now consider the case of two filaments in a side-by-side configuration at Re = 300. The non-dimensional values, the domain size and the initial angles (θ = 18 o ) are kept the same as in the case of the single beating filament. According to the experiments of [START_REF] Zhang | Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind[END_REF], varying the spacing between filaments d/L leads to the appearance of different filaments beating regimes. In particular, a symmetrical flapping is observed for distances d/L < 0.21. For higher values, a bifurcation towards a regime characterised by an out-of-phase flapping is detected. Additionally, the linear stability analyses carried out in [START_REF] Schouweiler | Coupled flutter of parallel plates[END_REF]; [START_REF] Michelin | Linear stability analysis of coupled parallel flexible plates in an axial flow[END_REF] have put forward the existence of three different modes for such configurations. Therefore, in this context we have considered various scenarios corresponding to different values of the spacing d/L.

Figure 3 displays the snapshots of iso-vorticity that we predict when considering three different spacings. The wakes are characterised by a periodic vortex shedding and by a flapping motion of the filaments (shown in fig- ure 4 for the three cases). The variations of the relative spacings between the filaments lead to different physical scenarios that are briefly reviewed hereafter.

• When the spacing is small (d/L = 0.1), we observe the mode M1, where the filaments are in close proximity and they behave almost as a single thick filament (see figure 3a), resulting in an in-phase beating of the filaments, as displayed in figure 4a.

• In contrast, when increasing the distance to d/L = 0.3, a different behaviour is observed. This mode (mode M2) is characterised by symmetrical out-of-phase oscillations, occurring after a transient period which occurs between t = 20 and t = 60 (see figure 4b). By increasing the filament spacing, the lock-in effect weakens but the interaction between the wakes generated by each filament still plays a dominant role, as shown in figure 3b. In this regime, the fluid enclosed between the filaments behaves like a flow generated by a pump due to the outof-phase flapping, cyclically being compressed when the two free ends approach (which is the case of the snapshot displayed in figure 3b), and released when they move apart.

• Further increasing the spacing to d/L = 1 results in a further wearkening of the wake interaction and a decoupling of the vortex streets behind the filaments (see figure 3c). However, beyond 5L downstream of the filaments tails, the vortices merge into a unique wake and the filaments reach the mode M2 characterised by an out-of-phase flapping (see figure 4c).

• If the spacing d/L is further increased, the two filaments eventually reach a totally decoupled dynamics with an in-phase flapping (mode M1).

The modal behaviour is consistent with the experimental observations of [START_REF] Zhang | Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind[END_REF] that report the onset of the anti-phase regime at d/L = 0.21, compared to our numerical predictions indicating a transitory regime occurring between d/L = 0.21 and d/L = 0.24. Note that for a sufficiently fine discretisation, the duration of this transcient is independent of the mesh refinement.

Retaining the same Reynolds number Re = 300, the configuration of three filaments placed side-by-side at an initial angle of 0 o is investigated. Figure 5 summarizes the different coupled dynamics obtained with the present simulations. The system follows the same behaviour as for the case of two filaments, except that an additional beating mode appears:

• for small spacings (d/L < 0.1), the mode M1 is observed, as in the case of two filaments, where the three filaments are in-phase (mode M1 in figure 5a);

• for d/L = 0.3, the two outer filaments flap out of phase while the inner filament is quasi stationnary (mode M2 in figure 5b); • for large spacing (d/L = 1.0) the outer filaments flap in-phase and the inner filament is out of phase (mode M3 in figure 5c);

• as for the case of two filaments, mode M1 is observed for very large spacing (d/L > 4.0) with an in-phase flapping of the three filaments.

Additionally, we observe a transition mode for d/L = 0.6 characterised by the same behaviour as mode M3 but with a low frequency modulation in the amplitude of the flapping of the filaments, as shown in figure 6. This transition mode has been reported in the numerical study of [START_REF] Tian | Coupling modes of three filaments in side-by-side arrangement[END_REF]. In their simulations at Re = 100, they also record a second transitional mode whereby the inner filament is observed to flap at a frequency reduced by half that of the outer filaments. Note that this second transitional mode is not observed in our simulations but the Reynolds number is different (Re = 300).

For cases where more than three filaments are considered (not reported here), the system is observed to exhibit further transitory modes resulting from the coupling between the baseline modes (M1, M2 and M3).

In-line flapping filaments

Next, we consider the collective behaviour of flapping filaments aligned with each other in the direction of an uniform fluid flow. This configuration is well know in the case of rigid bodies for the hydrodynamic drafting phenomenon which explains why one rigid body experiences a drag reduction downstream of another one, e.g. for cars or bicycles. It is of course related to the fact that the downstream body lies in the recirculating bubble of the upstream one, and is therefore subjected to lower velocities and thus lower fluid stresses. Peculiarly, these hydrodynamic properties due to the drafting effect are not straightforwardly inherited by flexible/flapping objects. While not formally attracting a great deal of attention to date, this issue has been recently put forward by the work of [START_REF] Ristroph | Anomalous Hydrodynamic Drafting of Interacting Flapping Flags[END_REF], who conducted soap film experiments in those configurations. In particular, they report an inverted drafing for in-line flapping filaments, as opposed to inline rigid bodies. A so called inverted drafting occurs where the upstream filament experiences a drag reduction instead of the downstream filament.

Here we aim to reproduce similar hydrodynamic conditions on aligned flexible filaments via a Direct Numerical Simulation considering a low Reynolds number as compared to the one used in experiments (i.e., Re = 300 versus almost 10 4 in experiments). The large difference in Reynolds number does not allow for direct comparisons between numerical and experimental results. Nonetheless, we are still able to give a qualitative analysis of the filaments interaction with special emphasis on the inverted drafting phenomenon. The non-dimensional values are kept the same as in the case of the single beating filament (and thus gravity is considered), except that the initial angles are set to θ = 0 o (aligned with the flow), and that the domain size is adapted in consequence, following the same requirements as the case of the single beating filament. In the following, the spacing s/L is defined as the distance between the tip of the leading filament and the pole of the trailing filament.

Figure 7 presents the beating patterns in the established regime at Re = 300 for two filaments spaced by a distance of s/L = 0.6. A good qualitative agreement with the soap film experiments is obtained although the flapping amplitudes are quite different. The amplitude of the upstream filament is found to be significantly smaller than the downstream one. Figure 8 shows for different spacings s/L, the drag ratio D/D 0 defined as the filament drag normalised by the reference drag obtained from an isolated filament at the same Reynolds number. For a spacing less than s/L = 1, we find for the upstream filament (circles) a reduced drag compared to the reference case, while the drag stays equal to the reference case for s/L > 1. Concerning the downstream filament (squares), apart from very close spacings (s/L < 0.3), a significant drag increase is observed compared to the reference case, up to 1.5D 0 . These numerical results confirm the anomalous drafting studied experimentally by [START_REF] Ristroph | Anomalous Hydrodynamic Drafting of Interacting Flapping Flags[END_REF]. For the set of parameters considered, the total drag force of the system is reduced (compared to two isolated filaments) for all spacings s/L < 0.5; and is otherwise increased. This trend is similar qualitatively to that obtained in figure 2 of [START_REF] Ristroph | Anomalous Hydrodynamic Drafting of Interacting Flapping Flags[END_REF] (although we recall that the Reynolds number is different as we focus here on the relative effect with respect to the isolated case).

For the small spacings (s/L < 0.5), the pair of filaments behaves as a single virtually longer filament, as the gap between the filaments enables a hydrodynamic link between the filament extremities (tail and pole). This effect tends to artificially increase the bending rigidity of the virtual longer filament. As a direct consequence the amplitude of the flapping of the upstream filament is reduced, as the amplitude of the flapping is lower near the filament pole. This assumption is confirmed by figure 9 which shows that the evolution of drag is correlated to the amplitude of flapping (measured as the maximum value of the excursion of the filament over time). Indeed a smaller upstream amplitude is associated on the figure to a smaller drag.

For large spacings (s/L > 0.5), as the upstream filament is hydrodynamically decoupled from the downstream one, its drag tends to the drag equivalent of the isolated case. However, the dowstream filament is found to have a higher drag as it experiences the perturbations induced by the wake of the upstream filament.

Following the experimental investigation of [START_REF] Ristroph | Anomalous Hydrodynamic Drafting of Interacting Flapping Flags[END_REF], we run the simulations for a chain of 6 in-line filaments with a spacing of s/L = 0.1, to assess if the proposed physical mechanisms hold true for more than two filaments. Figure 10 summarises the results showing that a drastic drag reduction is observed for the set of filaments: all six filaments experience drastic drag reductions compared to the reference drag of an isolated filament. The beating patterns observed experimentally are also qualitatively reproduced.

Conclusions

Using a model of flexible filament incorporating its flexural rigidity, the tension (enforcing inextensibility) and the added mass, we have successfully modeled numerically the dynamics of multiple flapping filaments immersed in a uniform flow.

When considering side-by-side filaments, the wake interactions and the modal behaviour of the system have been captured correctly, in agreement with experiments (and linear stability analysis). In the present work we have restricted our attention to the influence of the filament spacing, but the influence of the added mass ∆ρ plays also a significant role [START_REF] Michelin | Linear stability analysis of coupled parallel flexible plates in an axial flow[END_REF].

For the case of three filaments, a set of three baseline modes have been highlighted: in-phase flapping (M1), out-of-phase flapping with the inner filament at rest (M2), in-phase flapping with the inner filament flapping out of phase (M3). For the general case of a layer made of N filaments, we expect the system to be characterised by the appearance of N baseline modes originating from the combination of the M1, M2 and M3 baseline ones consistently with the theoretical prediction of [START_REF] Michelin | Linear stability analysis of coupled parallel flexible plates in an axial flow[END_REF].

For flapping filaments aligned with the flow, through our simulations, we have been able to confirm and partially characterise the inverted hydrody- namical drafting effect: the upstream filament is seen to experience a drag reduction compared to the downstream body. As far as drag is concerned, in agreement with the experimental observations of [START_REF] Ristroph | Anomalous Hydrodynamic Drafting of Interacting Flapping Flags[END_REF], we have found that for flexible objects, the paradigm "it is better to be chased than to chase" applies, as opposite to the rigid bodies case.

The physical mechanism behind this phenomenon is due to the passive shape adaptation to the flow of flexible filaments, tending to increase the bending rigidity of the upstream filament near the pole. The physical mechanism identified here is thus intrinsically passive, and as was also anticipated by [START_REF] Ristroph | Anomalous Hydrodynamic Drafting of Interacting Flapping Flags[END_REF], this effect does not necessarily extend to active flapping bodies (birds).

Figure 1 :

 1 Figure 1: Flapping motion of a single filament immersed in fluid at Re = 200, Ri = 0.5, ∆ρ = 1.5. Fluid flows from left to right. (a): beating pattern visualised by superimposed positions of the filament over one beating cycle. (b) Periodic time evolution of the ycoordinate of the free end.

  (a) Totally flexible filament (K B = 0). (b) Rigid filament (K B = 0.001).

Figure 2 :

 2 Figure 2: Comparison between instantaneous snapshots of the flapping filament without bending (a) and with bending (b) starting from a straight initial configuration at an angle of θ 0 = 18 o . The trajectory of the free end is shown in dashed line.

Figure 3 :

 3 Figure 3: Snapshots of iso-vorticity for the case of two beating filaments at Re = 300 and three different spacings. (a) mode M1 at d/L = 0.1, (b) mode M2 at d/L = 0.3, (c) mode M2 at d/L = 1.0.

Figure 4 :

 4 Figure 4: Time evolution of the y-coordinates of the free extremity of a system of two beating filaments at Re = 300. (a) Mode M1 at d/L = 0.1, (b) Mode M2 at d/L = 0.3, (c) Mode M2 at d/L = 1.0.

Figure 5 :

 5 Figure 5: Flapping patterns in the established regime for the beating of three filaments in a uniform flow for various spacing. (a) mode M1 at d/L = 0.05, (b) mode M2 at d/L = 0.3, (c) mode M3 at d/L = 1.0. The solid lines represent the time evolution of the y-coordinates of the free extremity of each filament at Re = 300.

Figure 6 :

 6 Figure 6: Transition mode observed between M2 and M3 for d/L = 0.6. The solid lines represent the time evolution of the y-coordinates of the free extremity of each filament at Re = 300.

Figure 7 :

 7 Figure 7: Beating patterns for a pair of flapping filaments in drafting configuration with a spacing of s/L = 0.6. The flow is going from left to right. Top: Experimental visualisations obtained by Ristroph and Zhang (2008) using thin-film interferometry. Bottom: present numerical results at Re = 300.

Figure 8 :

 8 Figure 8: Drag ratio D/D 0 as a function of the spacing s/L; circles: upstream filament; squares: downstream filament; triangles: sum of both drags.

Figure 9 :

 9 Figure 9: Drag ratio for upstream filament (plain line) and downstream filament (thick line) together with the amplitude of upstream filament (crosses) and downstream filament (circles), plotted as functions of the spacing s/L.

Figure 10 :

 10 Figure 10: Six flapping filaments in an in-line configuration. Top: experiments of Ristroph and Zhang (2008). Middle: present numerical results. Bottom: Relative drag ratio for each filament.

Acknowledgements

The financial help of the PELskin European project (FP7-TRANSPORT 334954) is greatly acknowledged.


Abdel-Hamid Soubra

Emilio Bastidas-Arteaga

A.-H Soubra

Advanced reliability analysis methods

come

Introduction

The most robust method used for the probabilistic analysis of geotechnical structures is the classical well-known Monte Carlo Simulation () methodology. It should be noted that the probabilistic analysis of an engineering system involves the computation of the PDF of the system response or the calculation of the failure probability for a prescribed threshold of this system response.

is not suitable for the computation of the small failure probabilities encountered in the practice of geotechnical engineering (especially when using a computationally-expensive finite element/finite difference deterministic model) due to the large number of simulations required to calculate a small failure probability. As an alternative to methodology, [START_REF] Au | Estimation of small failure probabilities in high dimentions by subset simulation[END_REF] proposed the subset simulation approach to calculate small failure probabilities. The basic idea of the approach is that the small failure probability can be expressed as a product of larger conditional failure probabilities.

Similarly, MCS is not suitable for the accurate determination of the PDF of a system response because of the great number of calls of the deterministic model, which are required for such a computation. The PCE methodology allows one to approximate a given system response by a polynomial chaos expansion (PCE) of a suitable order. Thus, the PCE methodology replaces the computationally-expensive deterministic model by a meta-model (i.e. a simple analytical equation). Once the PCE coefficients are determined, MCS can be applied on the obtained PCE to compute the PDF of the system response (and the corresponding statistical moments) with a quasinegligible computation time.

Subset simulation approach

Subset simulation approach was proposed by [START_REF] Au | Estimation of small failure probabilities in high dimentions by subset simulation[END_REF] as alternative to Monte Carlo Simulation () methodology to compute small failure probabilities. The basic idea of the approach is that the small failure probability can be expressed as a product of larger conditional failure probabilities. In this section, one presents a brief description of the steps of approach in case of two random variables, the extension to the case of several random variables being straightforward. A detailed description of the approach may be found in [START_REF] Au | Estimation of small failure probabilities in high dimentions by subset simulation[END_REF], in the chapters 1 and 4 of the book by [START_REF] Phoon | Reliability-based design in geotechnical engineering: Computations and applications[END_REF] and in [START_REF] Ahmed | Simplified and advanced approaches for the probabilistic analysis of shallow foundations[END_REF].

The steps of approach in case of two random variables ! , ! can be described as follows:

1. Generate a vector of two random variables ! , ! according to a target using direct Monte Carlo simulation. 2. Using the deterministic model, calculate the system response corresponding to ! , ! . 3. Repeat steps 1 and 2 until obtaining a prescribed number ! of vectors of random variables and the corresponding values of the system response. 4. Determine the value of the performance function corresponding to each value of the system response and then, arrange the values of the performance function in an increasing order within a vector ! where

! = ! ! , … , ! ! , … , ! ! ! . Notice
that the subscripts '0' refer to the first level (level 0) of the subset simulation approach. 5. Prescribe a constant intermediate conditional failure probability ! for the failure regions ! = 1,2, … , -1 and evaluate the first failure threshold ! which corresponds to the first level of approach (see Figure 1). The failure threshold ! is equal to the ! × ! + 1 !! value in the increasing list of elements of the vector ! . This means that the value of the conditional failure probability of the first level ! will be equal to the prescribed ! value. 6. Among the ! vectors of random variables, there are ! × ! ones whose values of the performance function are less than ! (i.e. they are located in the failure region ! ). These vectors are used as 'mother vectors' to generate ! new vectors of random variables (according to a proposal ! ) using Markov chain method based on the modified Metropolis-Hastings algorithm by [START_REF] Santoso | Modified Metropolis-Hastings algorithm with reduced chain-correlation efficient Probabilistic[END_REF]. This algorithm is presented in Appendix 1. 7. Using the deterministic model, calculate the values of the system response corresponding to the new vectors of random variables (which are located in level 1). Then, calculate the corresponding values of the performance function. Finally, gather the values of the performance function in an increasing order within a vector ! where

! = ! ! , … , ! ! , … , ! ! ! .
8. Evaluate the second failure threshold ! as the ! × ! + 1 !! value in the increasing list of elements of the vector ! . 9. Repeat steps 6-8 to evaluate the failure thresholds ! , ! , … , ! corresponding to the failure regions ! , ! , … , ! . Notice that contrary to all other thresholds, the last failure threshold ! is negative. Thus, ! is set to zero and the conditional failure probability of the last level ! !!! is calculated as:

! !!! = 1 ! ! ! ! ! ! ! !!! (1)
where ! ! = 1 if the performance function ! is negative and ! ! = 0 otherwise.

Figure 1: Nested Failure domain.

Finally, the failure probability is evaluated as follows:

= ! ! !!! ! !!! (2)
It should be mentioned that a normal was used herein as a target probability density function ! . However, a uniform was used as a proposal probability density function ! (for more details, refer to Appendix 1).

Example application

In order to illustrate the algorithm of methodology in a simple way, a numerical example is provided herein. In this example, approach was used to calculate the failure probability ! against bearing capacity failure of a strip footing of breadth .

The footing rests on a , soil and it is subjected to a service vertical load ! . The soil cohesion and the soil angle of internal friction were considered as random variables. The following formula was used for the computation of the ultimate bearing capacity:

! = 2 ! + ! + ! (3)
in which:

! = 2 ! -1 (4a) ! = !"#$% . ! 4 + 2 (4b) ! = ! -1 (4c)
where ! , ! and ! are the bearing capacity factors due to soil weight, surcharge loading and cohesion respectively. These coefficients are function of the soil friction angle. On the other hand, is the soil unit weight and is the surcharge loading.

The performance function used in the analysis is:

= ! ! -1 (5)
where ! is the ultimate footing load and ! is the footing applied load. As mentioned previously, only the soil cohesion and friction angle were considered as random variables. All the other parameters were considered as deterministic. These parameters are given in Table (1).

In this example, the intermediate failure probability ! of a given level = 1,2, … , -1 was arbitrary chosen equal to 0.2. A small number of samples per level ( ! = 10 samples) was used to facilitate the illustration.

Table 1: Data used for the probabilistic analysis of a strip footing against bearing capacity failure Table 2 presents (i) the values of and of each sample for the successive levels (ii) the corresponding values of the performance function and (iii) the values of the failure thresholds ! for the different levels. Notice that only the first two levels and the last level for which the failure threshold becomes negative were provided herein for illustration. Table (2) indicates that the failure threshold decreases with the successive levels until reaching a negative value at the last level. This means that the samples generated by the subset simulation successfully progress towards the limit state surface = 0. In order to select the failure threshold of a given level, the calculated values of the performance function of this level were arranged in an increasing order as shown in Table (2). Then, the failure threshold was selected as the ! × ! + 1 !! value of the arranged values of the performance function. Since ! = 10 and ! = 0.2, the failure threshold is equal to the third value of the arranged values of the performance function. The computation continues until reaching a negative value (or a value of zero) of the failure threshold. In this example, the negative value was reached in the sixth level (where ! = -0.0936) as shown in Table (2). Theoretically, the last failure threshold should be equal to zero. For this reason, ! was set to zero. This means that the last conditional failure probability ! ! is not equal to ! . In this case, the last conditional failure probability ! ! is calculated as the ratio between the number of samples for which the performance function is negative and the chosen number ! of samples (i.e. 10). According to Table (2), ! ! is equal to 3 10 = 0.3. Thus, the failure probability of the footing under consideration is equal to 0.2 ! ×0.3 = 9.6×10 !! . It should be emphasized that the failure probability calculated in Table (2) is not accurate due to the small value of ! . For an accurate computation of the failure probability, ! should be increased. This number should be greater than 100 to provide a small bias in the calculated ! value (see chapter 4 by Honjo in [START_REF] Phoon | Reliability-based design in geotechnical engineering: Computations and applications[END_REF]).

Parameter

In order to determine the optimal number of samples ! to be used per level, different values of ! were considered to calculate ! and its coefficient of variation ! ! as shown in Table (3). The thresholds corresponding to each ! value were calculated and were shown in this table. Table (3) indicates (as was shown before when ! = 10) that for the different values of ! , the failure threshold decreases with the successive levels until reaching a negative value at the last level. ! =2200 and 2400 samples are quasi similar as it may be seen from Table (3) and (ii) the corresponding final ! values are too close (they are respectively equal to 2.60×10 -3 and 2.63×10 -3 ).

Figure (2b)

shows the effect of ! on the coefficient of variation of the failure probability ! ! . As expected, ! ! decreases with the increase of ! . Notice that the values of ! ! for ! =2200 and 2400 samples are equal to 12.8% and 12.4% which indicates (as expected) that the ! ! decreases with the increase in the number of realizations.

It should be mentioned here that for ! = 0.2, four levels of subset simulation were found necessary to reach the limit state surface = 0 as may be seen from Table (3). Therefore, when ! =2200 samples, a total number of ! =2200×4=8800 samples were required to calculate the final ! value. Remember that in this case, the of ! was equal to 12.8%. Notice that if the same value of (i.e. 12.8%) is This means that, for the same accuracy, the approach reduces the number of realizations by 56%. On the other hand, if one uses with the same number of samples (i.e. 8800 realizations), the value of of ! would be equal to 19.6%. This means that for the same computational effort, the approach provides a smaller value of ! than .

!
0.1573 ------------------------------------------------ ! -0.0936 ------------------------------------------------ ! ×10 !! 0.

Polynomial Chaos Expansion (PCE) methodology

The basic idea of this method is to approximate a given system response by a polynomial chaos expansion (PCE) of a suitable order. In other words, the PCE methodology replaces the computationally-expensive deterministic model by a meta-model.

N s

In order to achieve this purpose, all the uncertain parameters (which may have different PDFs) should be represented by a unique chosen PDF. Table (4) presents the usual PDFs and their corresponding families of orthogonal polynomials. Within the framework of the present methodology, the response of a system that involves n random variables can be expressed by a PCE as follows:

Γ !"# = ! Ψ ! !!! !!! ξ (6) where ( )
ξ ψ i
are multi-dimensional polynomials defined as the product of onedimensional polynomials φ ! ! , ξ ! , ξ ! , … , ξ ! are independent random variables, ! , ! , … , ! are unknown coefficients to be evaluated and P is the size of the PCE.

The size P of the PCE (which is equal to the number of the unknown PCE coefficients) depends on the number n of random variables and the order p of the PCE. It is given as follows:

( )

! p ! n ! p n P + = (7)
It should be mentioned here that in this chapter, the random variables are represented in the independent standard normal space. Thus, the suitable corresponding bases are the multidimensional Hermite polynomials as may be seen from Table (4). The expressions of the multi-dimensional Hermite polynomials are given as follows:

∏ = = n 1 i i ) ( i ξ φ ψ α α (8)
where α=[α 1 , ….., α n ] is a sequence of n non-negative integers and

) ( i i ξ φ α
are onedimensional Hermite polynomials. More details on the one-dimensional Hermite polynomials are given in Appendix 2.

For the determination of the PCE unknown coefficients, a non-intrusive technique (in which the deterministic model is treated as a black-box) is used (see [START_REF] Ahmed | Simplified and advanced approaches for the probabilistic analysis of shallow foundations[END_REF] among others). In this chapter, the regression approach is employed. In this approach, it is required to compute the system response at a set of collocation points in order to perform a fit of the PCE using the obtained system response values.

As suggested by several authors (e.g. [START_REF] Huang | Geotechnical probabilistic analysis by collocation-based stochastic response surface method: An EXCEL add-in implementation[END_REF]), the collocation points can be chosen as the result of all possible combinations of the roots of the one-dimensional Hermite polynomial of order (p+1) for each random variable. For example, if a PCE of order p=2 used to approximate the response surface of a system with n=2 random variables, the roots of the one-dimensional Hermite Polynomial of order 3 are chosen for each random variable. These roots are (-√3, 0, √3) for the first random variable and (-√3, 0, √3) for the second random variable. In this case, 9 collocation points are available. These collocation points are (-√3, -√3), (-√3, 0), (-√3, √3), (0, -√3), (0, 0), (0, √3), (√3, -√3), (√3, 0), (√3, √3). In the general case, for a PCE of order p and for n random variables, the number N of the available collocation points can be obtained using the following formula:

N=(p+1) n (9)
Referring to Equations (7 and 9), one can observe that the number of the available collocation points is higher than the number of the unknown coefficients. This leads to a linear system of equations whose number N of equations is greater than the number P of the unknown coefficients. The regression approach is used to solve this system. This approach is based on a least square minimization between the exact solution Γ and the approximate solution Γ PCE which is based on the PCE. Accordingly, the unknown coefficients of the PCE can be computed using the following equation:

a = (Ψ T Ψ) -1 . Ψ T . Γ (10)
in which a is a vector containing the PCE coefficients, Γ is a vector containing the system response values as calculated by the deterministic model at the different collocation points and Ψ is a matrix of size N×P whose elements are the multivariate Hermite polynomials. It is given as follows:

Ψ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - - ) ( ....... .......... ) ( ) ( ) ( ) ( ....... .......... ) ( ) ( ) ( ) ( ....... .......... ) ( ) ( ) ( N 1 P N 2 N 1 N 0 2 1 P 2 2 2 1 2 0 1 1 P 1 2 1 1 1 0 ξ ψ ξ ψ ξ ψ ξ ψ ξ ψ ξ ψ ξ ψ ξ ψ ξ ψ ξ ψ ξ ψ ξ ψ ! ! ! ! ! ! ! ! (11)
Notice that in order to calculate the system response corresponding to a given collocation point, the standard normal random variables ξ i should be expressed in the original physical space of random variables as follows:

1 () i ix i xF ξ - =Φ ⎡⎤ ⎣⎦ (12)
in which, x i is a physical random variable, F is the CDF of the physical random variable and Ф is the CDF of the standard normal random variable. Notice also that if the original physical random variables are correlated, the standard normal random variables should first be correlated using the following equation: Once the PCE coefficients are determined, MCS can be applied on the obtained PCE (called meta-model) to compute the PDF of the system response. This is achieved by (i) generating a large number of realizations of the vector (ξ 1 , ξ 2 , … ξ n ) of standard normal random variables and (ii) calculating the system response corresponding to each realization by substituting the vector (ξ 1 , ξ 2 , … ξ n ) in the meta-model.

1 1 22 . . . . . . . c c n nc H ξ ξ ξξ ξ ξ ⎡⎤ ⎡⎤ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ = ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎣⎦ ⎣⎦ (13) in which { } nc c c ξ ξ ξ ..., , , 2

Optimal number of collocation points

The number of the available collocation points significantly increases with the increase in the number of random variables (cf. Eq. 9) and it may be very large with respect to the number of the unknown PCE coefficients. This makes it necessary to determine the optimal number of collocation points which is needed by the regression approach to solve the linear system of equations (Eq. 10). Sudret [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF] has proposed to successively increase the information matrix A [where A=(Ψ T Ψ)] until it becomes invertible as follows: (a) the N available collocation points are ordered in a list according to increasing norm, (b) the information matrix A is constructed using the first P collocation points of the ordered list, i.e. the P collocation points that are the closest ones to the origin of the standard space of random variables and finally (c) this matrix is successively increased (by adding each time the next collocation point from the ordered list) until it becomes invertible.

Accuracy of the obtained PCE

For a given PCE order, the accuracy of the approximation of the system response by a PCE can be measured by the coefficient of determination. Two types of coefficients of determination exist in literature. These are the coefficient of determination R 2 and the leave-one-out coefficient of determination Q 2 .

Let us consider J realizations of the standard normal random vector ξ as follows:

( ) ( ) ( ) { ( ) ( ( ) ) } J n J 1 ) J ( 1 n 1 1 ) 1 ( ..., , ..., , ..., , ξ ξ ξ ξ ξ ξ = =
, and let us assume that the vector

( ) ( ) ( ) ( ) { } J 1 ..., , ξ Γ ξ Γ Γ =
includes the corresponding values of the system response determined by deterministic calculations. The coefficient of determination R 2 is calculated as follows:

PCE R Δ - = 1 2 (14)
where PCE Δ is given by:

( ) ( ) ( ) ( ) ( ) ( ) 2 1 1/ J ii PCE i PCE J Var ξξ = ⎡⎤ Γ-Γ ⎣⎦ Δ= Γ ∑ (15) and ( ) ( ) (
)

2 1 1 1 J i i Var J ξ = ⎡⎤ Γ= Γ -Γ ⎣⎦ -∑ (16)
Note here that J is the number of collocation points used to evaluate the unknown coefficients of the PCE. The value 1 R 2 = indicates a perfect approximation of the true system response Γ, whereas 0 R 2 = indicates a nonlinear relationship between the true model Γ and the PCE model Γ PCE .

The coefficient of determination R 2 may be a biased estimate since it does not take into account the robustness of the meta-model (i.e. its capability of correctly predicting the model response at any point which does not belong to the collocation points. As a consequence, a more reliable and rigorous coefficient of determination, called the leave-one-out coefficient of determination, was proposed in literature. This coefficient of determination consists in sequentially removing a point from the J collocation points. Let Γ ξ/i be the meta-model that has been built from (J-1) collocation points after removing the i th observation from these collocation points and let

) ( ) ( ) i ( /i ) i ( i ξ Γ ξ Γ Δ ξ - =
be the predicted residual between the model evaluation at point ξ (i) and its prediction at the same point based on Γ ξ/i . The empirical error is thus given as follows:

* 2 1 1 ()

J i PCE i J = Δ= Δ ∑ (17)
The corresponding coefficient of determination is often denoted by Q 2 and is called leave-one-out coefficient of determination. It is given as follows:

* 2 1 ()

PCE Q Var Δ =- Γ (18)

PCE-based Sobol indices

A Sobol index of a given input random variable is a measure by which the contribution of this input random variable to the variability of the system response can be determined. Sobol indices are generally calculated by MCS methodology. This method is very time-expensive especially when dealing with a large number of random variables. [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF] proposed an efficient approach to calculate the Sobol indices based on the coefficients of the PCE. This method is based on ranking the different terms of the PCE and gathering them into groups where each group contains only one random variable or a combination of random variables. For more details on the computation of Sobol indices, the reader may refer to [START_REF] Ahmed | Simplified and advanced approaches for the probabilistic analysis of shallow foundations[END_REF] among others.

Conclusion

This chapter first presented the subset simulation approach which is an efficient alternative to for the computation of a small failure probability. An example application was provided. It aims at showing the practical implementation of the approach. It was found that for a prescribed accuracy, the approach significantly reduces the number of realizations as compared to Monte Carlo simulations methodology (the reduction was found equal to 93.3% in the present chapter). In other words, for the same computational effort, the approach provides a smaller value of the coefficient of variation of ! than . It should be mentioned that the Matlab code used for the example application is provided in http://www.univnantes.fr/soubra-ah for practical use.

In a second stage, the Polynomial Chaos Expansion methodology was presented. It was shown that the PCE method replaces the computationally-expensive deterministic model by a meta-model (i.e. a simple analytical equation). Once the PCE coefficients are determined, MCS can be applied on the obtained PCE to easily compute the PDF of the system response with a quasi-negligible computation time.

Figure 2 :

 2 Figure 2: ! and ! ! versus the number of realizations ! .

1

  is the vector of correlated standard normal random variaof uncorrelated standard normal random variables and H is the Cholesky transformation of the correlation matrix of the physical random variables.

Table 2 :

 2 Results of algorithm when ! =10 and ! =0.2

			Mean and coefficient
		Type of parameter	of variation of the pa-
			rameter
	Breadth	Deterministic	2m
	Surcharge loading	Deterministic	10kPa
	Soil unit weight	Deterministic	20kN/m 3
	Service vertical load	Deterministic	1000kN/m
	!
	Cohesion	Random normal varia-ble	! = 20kPa ! = 0.3
	Friction angle	Random normal varia-ble

! = 30° ! = 0.1

Table 3 :

 3 Evolution of the failure threshold ! with the different levels and with the number of realizations ! when ! = 0.2 Figure (2a) shows the effect of ! on the failure probability. It indicates that for small values of ! , the failure probability largely changes with ! . However, for high values of ! , the failure probability converges to an almost constant value. Figure (2a) also indicates that 2200 samples per level are required to accurately calculate the failure probability. This is because (i) the ! values corresponding to

Table 4 :

 4 Usual probability density functions and their corresponding families of orthogonal polynomials

	Probability density functions	Polynomials
	Gaussian	Hermite
	Gamma	Laguerre
	Beta	Jacobi
	Uniform	Legendre

APPENDIX 1 Modified METROPOLIS-HASTINGS algorithm

The Metropolis-Hastings algorithm is a Markov chain Monte Carlo () method. It is used to generate a sequence of new realizations from existing realizations (that follow a target called ' ! '). Refer to Figure (1) and let ! ∈ ! be a current realization which follows a target ' ! '. Using a proposal ' ! ', a next realization !!! ∈ ! that follows the target ' ! ' can be simulated from the current realization ! as follows:

a. A candidate realization is generated using the proposal ( ! ). The candidate realization is centered at the current realization ! . b. Using the deterministic model, evaluate the value of the performance function corresponding to the candidate realization . If < ! (i.e. is located in the failure region ! ), set !!! = ; otherwise, reject and set Notice that in step , if the candidate realization does not satisfy the condition < ! , it is rejected and the current realization ! is repeated. Also in step , if the candidate realization does not satisfy the condition ≥ 1 (i.e. is not distributed according to the ! ), it is rejected and the current realization ! is repeated. The presence of several repeated realizations is not desired as it leads to high probability that the chain of realizations remains in the current state. This means that there is high probability that the next failure threshold !!! is equal to the current failure threshold ! . This decreases the efficiency of the subset simulation approach. To overcome this inconvenience, [START_REF] Santoso | Modified Metropolis-Hastings algorithm with reduced chain-correlation efficient Probabilistic[END_REF] proposed to modify the classical M-H algorithm as follows:

a. A candidate realization is generated using the proposal ! . The candidate realization is centered at the current realization ! . b. Calculate the ratio ! = ! ! ! and the ratio ! = ! ! ! ! , then compute the value = ! ! . c. If ≥ 1, set !!! = ; otherwise, another candidate realization is generated.

Candidate realizations are generated randomly until the condition ≥ 1 is satisfied. d. Using the deterministic model, evaluate the value of the performance function !!! of the candidate realization that satisfies the condition ≥ 1.

If !!! < ! (i.e. !!! is located in the failure region ! ), one continues to retain the realization !!! obtained in step ; otherwise, reject and set !!! = ! (i.e. the current realization ! is repeated).

These modifications reduce the repeated realizations and allow one to avoid the computation of the system response of the rejected realizations. This becomes of great importance when the time cost for the computation of the system response is expensive (i.e. for the finite element or finite difference models).

APPENDIX 2

The one-dimensional Hermite polynomials of orders 0, 1, 2, 3, …, p+1 are given by:


Mouchira Labidi
email: mouchira.labidi@promes.cnrs.fr

Julien Eynard
email: julien.eynard@promes.cnrs.fr

Olivier Faugeroux
email: olivier.faugeroux@promes.cnrs.fr

Stéphane Grieu
email: stephane.grieu@promes.cnrs.fr

Predictive Control and Optimal Design of Thermal Storage Systems for Multi-energy District Boilers

Keywords: Multi-energy district boiler, thermal storage tank, optimal design, model predictive control, management strategy

As part of the second phase of the OptiEnR research project, the present work deals with improving the operation of a multi-energy district boiler by adding to the plant an optimally designed and controlled thermal storage tank. Previous study focused on both a design approach, based on a parametric analysis, and a non-predictive control strategy. The aim of the present work was to develop a Model Predictive Controller (MPC) to improve the management of the tank in real time. The proposed controller generates optimal command sequences dealing with the amount of thermal energy to be stored or released. As a result, both the fossil energy consumption and CO 2 emissions are significantly reduced while the economic gain is increased.

INTRODUCTION

According to the International Energy Agency (2013), 31.6 billion tons of carbon dioxide were emitted worldwide last year, what represents an increase of 1.4%. As a result, the organisation has made three main recommendations: (i) efforts should be focused on improving energy efficiency in buildings, transport and industry, (ii) fossil fuels must be replaced by low-emission sources of energy, (iii) the release of methane has to be reduced (unburned natural gas with high greenhouse effect) in oil and gas industries.

Because of the global energy crisis, the French government supports renewable energy production. As buildings account for about 40% of total final energy consumption (more than half of this consumption is used for heating), France makes a specific effort in this sector. In addition, using biomass materials such as wood in industrial and residential heating can significantly reduce the reliance on fossil fuels and limit CO 2 emissions [START_REF] Kitzing | Renewable energy policies in Europe: Converging or diverging?[END_REF]. In financial terms, biomass is cheaper than many fossil fuels commonly used for heating. As a key point, advanced control techniques and management strategies are needed to improve the operation of multi-energy district boilers. In this sense, the OptiEnR research project focuses on optimizing the performance of the plants operated by Cofely GDF-Suez, our industrial partner, by adding thermal storage systems.

Thermal energy storage is an attractive technology used in several industrial plants such as Combined Heat and Power (CHP) plants [START_REF] Taljan | Optimal sizing of biomass-fired Organic Rankine Cycle CHP system with heat storage[END_REF], Central Solar Heating (CSH) plants [START_REF] Rodríguez-Hidalgo | Domestic hot water consumption vs. solar thermal energy storage: The optimum size of the storage tank[END_REF] or multienergy district boilers (Eynard et al., 2011a[START_REF] Eynard | Predictive control and thermal energy storage for optimizing a multi-energy district boiler[END_REF]. It has been highlighted in a previous work [START_REF] Labidi | Optimal design of thermal storage tanks for multi-energy district boilers[END_REF] that once optimally designed and managed (we proposed a non-predictive strategy), a storage tank can improve in a significant way the overall efficiency of a plant. So, the present paper deals with the optimal management of a thermal storage tank using a Model Predictive Controller (MPC). We focused on the analysis of the energy savings one can achieve thanks to such an advanced control approach. First, the multi-energy district boiler we selected is described (section 2). Next, the reliable (nonpredictive) management strategy which was presented in the above-mentioned work is presented. The advantages and limitations are outlined (section 3). They account for the use of a predictive strategy. Then, the design of the MPC controller allowing the thermal storage tank to be efficiently managed is carried out (section 4). Lastly, a comparison between both strategies is performed (section 5). As a result, it can be highlighted that the way the tank is controlled is a key factor. Indeed, the efficiency of such a system is mainly related to its design and the way it is managed. The paper ends with a conclusion and an outlook to future work (section 6).

MULTI-ENERGY DISTRICT BOILER

We selected as a case study a multi-energy district boiler managed by Cofely GDF-Suez and located in the northeast of France, in the Alsace region (Haut-Rhin). Alsace has a semi-continental climate with cold and dry winters and hot summers. In addtion, there is little precipitation. The plant is connected to a heat network for thermal energy distribution (Fig. 1).

WB

Network GB1 The biomass (wood) boiler (WB) is characterized by a minimum (P min W B ) and a maximum heat power (P max W B ) as well as a variable efficiency (η). Efficiency, for each operation point, can be defined using interpolation. Usually, P min W B is higher or equal to 25% of P max W B . The biomass boiler is designed in order to meet the major part of the power demand but is not able to cover the peak loads. As a result, a 7000 kW gas boiler (GB1) is used as a backup unit. In addition, a 9000 kW gas boiler (GB2) is switched on in case of malfunction or during maintenance phases. Unlike biomass boilers, efficiency is the same whatever the power. P min GB is about 2% of P max GB . Table 1 summarizes the characteristics of the boiler units: and P max W B to meet the demand. When it is lower than P min W B , the biomass boiler is switched off and the backup boiler (GB1) is used. During the coldest periods of winter, both boilers (WB and GB1) operates jointly. In this case, WB operates at full load.

GB2 Storage

The main purpose of the previous study we carried out [START_REF] Labidi | Optimal design of thermal storage tanks for multi-energy district boilers[END_REF] was to optimize the operation of multienergy district boilers by adding optimally-sized thermal storage tanks and proposing an adequate (sequential) management strategy. As a key point, the main goal of thermal energy storage is to avoid switching on the gas boiler (GB1) as far as possible. Such a process leads to significant changes in fossil and renewable energy consumption. So, a multi-energy district boiler equipped with a thermal storage tank can be operated using the strategy described as follows, on the basis of an estimate of the power demand (P net ) and the characteristics of the boiler units:

• During the coldest months of winter (P net < P max W B ), instead of modulating its power, the biomass boiler operates at full load to meet the power demand and charge the thermal storage tank. Once such a demand is upper than P max W B , the stored energy is released. In this way, the gas boiler is only switched on when the tank is empty and the power demand still exceeds P max W B .

• During the hottest months of summer (P net < P min W B ), most of the buildings do not need to be heated and, as a consequence, only domestic hot water is required. Generally, biomass boilers are oversized in order to be able to operate during this period of the year and a gas boiler is used to meet low power requirements. A biomass boiler combined with a thermal storage tank can be used during such a period as follows: first the biomass boiler operates at minimum power and allows both the power demand to be met and the tank to be charged. Once the storage tank is completely filled, the boiler is shut down and the stored energy is released to afford domestic hot water. The boiler is switched on again when the tank is empty. This operating mode prevents the use of gas and favours the use of renewable energy.

Advantages and limitations

The main advantages of such a strategy are listed below. First, it is based on logical conditions and, as a result, no complex algorithm involving extensive calculation is needed. Secondly, it allows the specificities of each plant (through profiles of power demand and technical characteristics) to be taken into account. As a result, this strategy can be applied to any (existing or under-construction) multi-energy district boiler.

The main drawback of the strategy lies in not taking into account the future power demand. Consequently, the thermal storage tank cannot be used in an optimal way and, in some cases, it is unable to cope with peak loads. In addition, it can be sometimes full and not used for a long time. In this case, thermal losses occur. That is why a predictive strategy based on a MPC controller is likely to improve operation and performance.

MODEL PREDICTIVE CONTROL (MPC)

Principles of MPC

It is somewhat curious to note that the concept of Modelbased Predictive Controller (MPC) has a long history that began during the 1970's when Engineers at Shell Oil developed their own dependent MPC technology with an initial application in 1973 [START_REF] García | Model predictive control: Theory and practice -a survey[END_REF]. Nowadays, this concept is widely used in the control of industrial processes. Its popularity in industry is mainly due to the possibility it offers to treat operating specifications and constraints jointly during the development phase of the controller. MPC is commonly used to manage thermal comfort [START_REF] Castilla | Thermal comfort control using a non-linear mpc strategy: A real case of study in a bioclimatic building[END_REF][START_REF] Pravara | Role of MPC in building climate control[END_REF] and energy resources [START_REF] Ma | Demand reduction in building energy systems based on economic model predictive control[END_REF][START_REF] Kim | Building demand-side control using thermal energy storage under uncertainty: An adaptive multiple model-based predictive control (mmpc) approach[END_REF] in buildings. [START_REF] Eynard | Predictive control and thermal energy storage for optimizing a multi-energy district boiler[END_REF] developed a predictive controller in order to optimize the operation of a multi-energy district boiler located in northwest France.

As it is well known, the philosophy of MPC is down to use a model to forecast the behavior of the system to be controlled and choose the best decision in the sense of some objective function (J) while satisfying the constraints. Usually, the aim is to ensure the desired set-point regardless of disturbances with minimal effort. Constraints deal with physical limitations and are introduced for economic or security reasons. The forecast horizon is the time interval during which the objective function is minimized thanks to an optimization algorithm [START_REF] Manenti | Considerations on nonlinear model predictive control techniques[END_REF].

Design of the model predictive controller

In this section of the paper, the MPC problem is formulated. The aim is to optimize the control input that minimizes the consumption of energy (or operational costs) while meeting power demand requirements. The proposed model predictive controller defines the amount of energy to be stored or released through the storage tank at each time step. The architecture of the predictive control approach we propose to manage the amount of thermal energy stored in the tank is depicted by Fig. 2: Model of the multi-energy district boiler.

P TS (k/k, …, k + 23/k)
The MPC controller we designed uses the global model of the multienergy district boiler described in section 2 [START_REF] Labidi | Optimal design of thermal storage tanks for multi-energy district boilers[END_REF]. The gas power (P gas ) as well as the wood power (P wood ) (i.e. the model outputs) are defined as functions of both the tank charging/discharging power (P T S ) and the power demand (P net ) (i.e. the model inputs). P net is forecasted using a wavelet-based multi-resolution analysis and feedforward artificial neural networks [START_REF] Eynard | Waveletbased multi-resolution analysis and artificial neural networks for forecasting temperature and thermal power consumption[END_REF]. The main idea behind such an approach is to replace the forecasting of an original time series whose variability can be high by the forecasting of its wavelet coefficients (of lower variability). We decided for a forecast horizon (H f ) of 24 hours, as an interesting compromise between forecasting accuracy and periodicity in power demand. In addition, such an horizon is well adapted to the charging and discharging cycles of the tank.

Optimization variables. The proposed controller defines at each time step k the optimal values of the tank charging/discharging power (P T S ), along the forecast horizon.

Let us note that a positive value (i.e. P T S (k) > 0) stands for the "charging (storage) mode" while a negative value (i.e. P T S (k) < 0) is for the "discharging (release) mode".

Objective function.

The main goal of the control approach is to minimize the use of fossil energy by optimizing the storage of renewable energy during low-demand periods and releasing the energy stored when the power demand is high. Thus, the objective function J is defined as the quadratic sum of the gas power consumed at each time step k along the forecast horizon. J is depicted by equation 1:

J = H f k=1 P gas(k) 2
(1)

Constraints. Equations 2 and 3 allow the optimization constraints to be formulated. Such constraints ensure that thermal energy is stored or released adequately. The first one ( 2) is introduced in order to limit the interval of the possible values for P T S . This constraint is related to the characteristics of the storage tank feed pumps. The second constraint (3) makes reference to the design of the tank.

In other words, it is related to the capacity of the tank (E max ). E max = ρ ⋅ C p ⋅ V ⋅ ∆T with ρ (kg m3) the water density, C p (kJ kg⋅K) the specific heat of water, ∆T (K) the difference in temperature between cold and hot water and V the volume of the tank (m 3 ). At each time step k, the amount of energy stored has to be positive and lower than E max :

-P max T S ≤ P T S (k + j k) ≤ P max T S ∀j ∈ 0; H f -1 (2) 0 ≤ E init + j i=0 P T S (k + i k) ≤ E max ∀j ∈ 0; H f -1 (3)
Optimization problem. The optimization problem comes down to find, at each time step k (t = k × T s with T s the sampling time), the value of the manipulated variable P T S such that the function J is minimized and the constraints are satisfied (4 and 5): min (J)

[P T S (k k),...,P T S (k+H f -1 k)] (4) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ District boiler model with thermal storage tank -P max T S ≤ P T S (k + j k) ≤ P max T S 0 ≤ E init + ∑ j i=0 P T S (k + i k) ≤ E max ∀j ∈ 0; H f -1 and H f = 24 (5) MPC algorithm.
Fig. 3 depicts the MPC algorithm used. At each time step k, a simulation over the forecast horizon based on the non-predictive strategy, the current amount of energy stored in the tank and the forecasted values of P net is performed in order for the values of [P T S (k k), . . . , P T S (k + H f -1 k)] to be initialized. These values are then optimized using both the non-linear optimization algorithm "fmincon" from Matlab ® and the developed model of the multi-energy district-boiler. The first optimized value is applied to the model which stands for the real system and so on until the end of the simulation process. Then, various performance criteria are computed for an off-line analysis of the results.

RESULTS

This section of the paper deals with the control results we obtained. Three different scenarios are considered in order to check performance:

• Scenario 1: District boiler without thermal storage tank. This is the "Reference Scenario" (RS). • Scenario 2: District boiler with thermal storage tank managed with the "Non-Predictive Strategy" (NPS) discussed in section 3 of the paper. • Scenario 3: District boiler with thermal storage tank managed with the MPC strategy discussed in section 4 of the paper.

Energy, economic and environmental criteria are proposed as performance indicators. Thanks to these indicators and considering a thermal storage tank whose size can vary from 0 m 3 to 1000 m 3 , one can evaluate the performance of the proposed strategies (predictive and non-predictive) as well as the impact of the volume of the tank on the multi-energy district boiler operation.

Performance indicators and simulation parameters

Because the main purpose of thermal energy storage is to decrease gas consumption in multi-energy district boilers, the gas coverage rate (C gas ) is proposed as an energy indicator. It is calculated from the thermal energy produced from the combustion of gas (E gas ) and wood (E wood ) during the considered period (6).

C gas = E gas E wood + E gas (6)
The wood coverage rate (C wood ) is computed in the same way than C gas and is subject to contract (7):

C wood = E wood E wood + E gas (7)
With the aim of highlighting the economic benefits of energy savings, a criterion (Ec) is defined from E gas , E wood , U P gas (the unitary price of gas), and U P wood (the unitary price of wood) (8):

Ec = E gas × U P gas + E wood × U P wood (8)
In order to put in perspective the economic benefits of energy savings, the economic gain G is evaluated (9). It is defined as the difference between Ec (V ), the economic cost related to energy consumption, considering a storage volume V , and Ec (V = 0), the economic cost related to energy consumption without storage of thermal energy (V = 0 m 3 ):

G(V ) = Ec (V ) -Ec (V = 0) (9)
The environmental impact of such a technology is evaluated thanks to criterion L CO2 , which is about CO 2 emissions (10). L CO2 is expressed from Ec gas , Ec wood and the Life-Cycle Assessment of CO 2 emissions from gas (U gas CO2 ) and wood (U wood CO2 ):

L CO2 = Ec gas × U gas CO2 + Ec wood × U wood CO2 (10)
It should also be noted that new buildings connected to the heat network as well as future expansions of existing buildings are factors to be taken into account in order to evaluate the proposed strategy accurately. As a result, we considered an increase in the power demand (P net ) up to 30% and studied the impact of such an increase on the plant operation. Finally, table 2 summarizes the main simulation parameters.

Simulation results

First, simulation results show that the biomass boiler is sized to ensure around 85% of the power demand during the simulation period without thermal storage tank (V = 0 m 3 ). Fig. 5 depicts the way the gas coverage rate evolves according to the volume of the storage tank and the strategy used (NPS or MPC). Let us remember that NPS is for "Non-Predictive Strategy" while RS stands for "Reference Scenario". One can highlight that, regardless of the strategy used, adding to a plant a thermal storage tank reduces the use of the gas boiler significantly. Moreover, whatever the size of the tank, the predictive strategy is always more efficient. For example, when adding to the plant a 200 m 3 thermal storage tank managed using the non-predictive strategy, the gas coverage rate decreases down to 7.3%. The same tank managed using the predictive strategy allows this rate to be decreased until 6%.

Fig. 4 presents the dynamics of the district boiler during one week. Taking a look at the figure, it can be highlighted that the predictive controller anticipates very well the storage of thermal energy even if the current power demand is very low, knowing 24 hours in advance that such a demand will increase sharply. Figures 6 and 7 depict, for both strategies, the way the size of the tank impacts on Ec and G, respectively. One can highlight that Ec decreases with the size of the thermal storage tank. Furthermore, in comparison to the non-predictive strategy, the MPC strategy allows a considerable economic gain to be realized.

With the non-predictive strategy, the optimal volume of the tank is about 200 m 3 . A more important volume is not reflected by a bigger economic gain. With the MPC strategy, the optimum is reached for a volume of about 300 m 3 . This highlights that such a strategy takes advantage of an increase in the size of the tank and demonstrates its ability in managing energy resources.

From an environmental point of view, Fig. 8 clearly shows for the same thermal storage tank that the MPC strategy allows a higher decrease in CO 2 emissions than the nonpredictive strategy. A possible increase in power demand is also considered in order to complete the analysis. Regarding the ability of the thermal storage tank to cope with such an increase, one can observe (Fig. 9) that the predictive strategy allows the wood coverage rate to be about 2 points higher than with the non-predictive strategy, for a volume of 200 m 3 , and from 4 to 10 points higher than with the reference scenario (i.e. no thermal storage tank).

CONCLUSION

The present paper deals with optimizing the operation of a multi-energy district boiler by adding to the plant an optimally-sized thermal storage tank. First, we proposed a non-predictive strategy in order to manage the plant, on the basis of the power demand and the characteristics of the different boiler units. Next, a model predictive controller has been developed in order to optimize, over a forecast horizon, the use of the storage tank. As a key point, one can highlight that the control scheme allows both the fossil energy consumption and CO 2 emissions to be significantly reduced. In addition, the economic gain is increased. Ongoing research activities will now focus on improving the tool used to forecast the power demand.

Other multi-energy district boilers and configurations will also be considered in order to validate the predictive approach. Future work will also deal with implementing in situ the developed algorithm.

Fig. 1 .

 1 Fig. 1. Synoptic of the multi-energy district boiler.

  Fig. 2. Structure of the model-based predictive controller.

  Fig. 3. MPC algorithm.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Weekly dynamics of the multi-energy district boiler.

Fig. 6 .

 6 Fig. 6. Impact of the size of the tank on Ec.

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Impact of the size of the tank on G.

Table 1 .

 1 Characteristics of the boiler units.

	Boiler unit	P [kW]	η
		1050	0.7
	WB	2100 3150	0.75 0.88
		4200	0.95
	GB1	140 7000	0.97 0.97
	GB2	180 9000	0.97 0.97
	3. NON-PREDICTIVE STRATEGY (NPS)
	3.1 Overview of the strategy
	Usually, the biomass power is modulated between P min W B

Table 2 .

 2 Simulation parameters.

	Parameter	Value
	Simulation period	From September 21 to April 16
	Sampling time (Ts)	1 hour
	Forecast horizon (H f )	24 hours
	P max T S	Emax 2	[kW]
	U Pgas	40 e⋅MWh -1
	U P wood	17 e⋅MWh -1
	U gas CO 2	234 kgCO 2 ⋅MWh -1
	U wood CO 2	13 kgCO 2 ⋅MWh -1