Damien Jallas

Olivier Marquet

David Fabre

Linear and non-linear perturbation analysis of the symmetry-breaking in time-periodic propulsive wakes

Keywords: flapping foil, spatio-temporal symmetry, Floquet stability

The two-dimensional and time-periodic wake flows produced by a pitching foil are investigated numerically for a fixed flapping amplitude. As the flapping frequency is increased, three regimes are identified in the time-marching non-linear simulations. The first regime is characterised by nondeviated wake flows with zero time-averaged lift. In the second regime, the wake flow is slightly deviated from the streamwise direction and the time-averaged lift is slightly positive or negative. The third regime is characterised by larger deviations of the wake, associated with larger values of both the time-averaged lift and the thrust. The transition from the first to the second regime is examined by performing a Floquet stability analysis of the non-deviated wake. A specific method is introduced to compute the time-periodic, non-deviated wake when it is unstable. It is found that one synchronous anti-symmetric mode becomes unstable at the critical frequency where deviation occurs. Investigation of its instantaneous and time-averaged characteristics show that it acts as a displacement mode translating the non-deviated wake away from the streamwise direction. Finally, it is demonstrated that the transition from the second to the third regime is linked to non-linear effects that amplify both anti-symmetric and symmetric perturbations around the foil.

I. INTRODUCTION

The locomotion of living animals such as fish and birds is a source of inspiration for researchers and engineers designing underwater or aerial vehicles. By mimicking the kinematics of insect wings, they expect that flappingbased Micro Air Vehicles will benefit from the amazing capabilities of insects' flight: vertical take-off, hovering and slow-forward flight, and low-acoustic signature [START_REF] Jm Anderson | Oscillating foils of high propulsive efficiency[END_REF][START_REF] Sanjay | The aerodynamic effects of wing rotation and a revised quasi-steady model of flapping flight[END_REF][START_REF] Lentink | The scalable design of flapping micro-air vehicles inspired by insect flight[END_REF].

The aerodynamics of an insect's flapping wings is very complex, not only because the kinematics of the wing is highly unsteady and three-dimensional [START_REF] Lighthill | Aerodynamic aspects of animal flight[END_REF][START_REF] Taylor | Animal flight dynamics ii. longitudinal stability in flapping flight[END_REF], but also because the surrounding flow is described by the non-linear Navier-Stokes equations. Consequently, even when the kinematics of a two-dimensional wing is a purely harmonic two-dimensional pitching or heaving motion, the periodic wake-flow can spontaneously become three-dimensional [START_REF] Leontini | Three-dimensional transition in the wake of a transversely oscillating cylinder[END_REF][START_REF] Rafael | Floquet stability analysis of the flow around an oscillating cylinder[END_REF][START_REF] Deng | Three-dimensional transition after wake deflection behind a flapping foil[END_REF]. Various two-dimensional periodic flow-patterns have also been observed in the wake of the wing [START_REF] Manoochehr | Vortical patterns in the wake of an oscillating airfoil[END_REF], [START_REF] Schnipper | Vortex wakes of a flapping foil[END_REF], depending on the flapping frequency and amplitude. For small values of amplitude and frequency, the classic Bénard-Von Karman vortex street is observed in the wake and the aerodynamic force exerted on the wing is resistive (drag force). For high values, a reversed Bénard-Von Karman vortex street is observed in the wake and the aerodynamic force exerted on the wing is propulsive (thrust force). For intermediate values, many other wake patterns have been observed experimentally [START_REF] Manoochehr | Vortical patterns in the wake of an oscillating airfoil[END_REF][START_REF] Chk Williamson | Vortex formation in the wake of an oscillating cylinder[END_REF] and numerically [START_REF] Dong | Wake topology and hydrodynamic performance of low-aspect-ratio flapping foils[END_REF] (see [START_REF] Schnipper | Vortex wakes of a flapping foil[END_REF] for a detailed classification). The present paper investigates one of the periodic wake flow-patterns, characterised by the mean deviation of the reverse Bénard-Von Karman vortex street from the streamwise direction.

To our knowledge, the first experimental observation of an asymmetric pattern in the wake of pitching wings was made by Bratt [START_REF] Bratt | Flow patterns in the wake of an oscillating aerofoil[END_REF]. In this case, the wake-flow breaks the spatio-temporal symmetry of the wing kinematics. Later on, Jones et al. investigated experimentally [START_REF] Jones | Experimental and computational investigation of the knoller-betz effect[END_REF] and numerically [START_REF] Jones | Wake structures behind plunging airfoils: a comparison of numerical and experimental results[END_REF] the deviated reverse Bénard-Von Karman vortex street in the wake of a heaving foil. Results of the two-dimensional numerical simulations compare well with experimental results. Interestingly, both upward and downward deviation of the wake have been reported for the same values of flapping frequency and amplitude. The deviation direction was shown to depend on the initial flapping conditions [START_REF] Lai | Jet characteristics of a plunging airfoil[END_REF][START_REF] Ellenrieder | Piv measurements of the asymmetric wake of a two dimensional heaving hydrofoil[END_REF]. Zheng et al. [START_REF] Zc Zheng | Study of mechanisms and factors that influence the formation of vortical wake of a heaving airfoil[END_REF] observed that the deviation's direction results from a competition between the first three vortices that are emitted by the foil. The deviation of the propulsive wake behind a pitching wing was investigated experimentally by Godoy-Diana et al. [START_REF] Godoy-Diana | Transitions in the wake of a flapping foil[END_REF]. Using only two-dimensional measurements of the velocity field behind the wing, they identified the region of asymmetric wake in the amplitude and frequency parameter space. Later [START_REF] Godoy-Diana | A model for the symmetry breaking of the reverse bénard-von kármán vortex street produced by a flapping foil[END_REF], they proposed a criterion based on the relation between the experimentally determined phase velocity of the vortex street and an idealized self-advection velocity of two consecutive counter-rotating vortices in the near-wake. The disappearance of asymmetric wakes for low-flexibility pitching foils was first investigated by Marais et al. [START_REF] Marais | Stabilizing effect of flexibility in the wake of a flapping foil[END_REF] and then by Zhu et al. [START_REF] Zhu | How flexibility affects the wake symmetry properties of a selfpropelled plunging foil[END_REF], who also showed that highly-flexible pitching foils can instead enhance the deviation. Experimental investigations of the wake deviation behind heaving wings were performed by Cleaver et al. [START_REF] David J Cleaver | Vortex mode bifurcation and lift force of a plunging airfoil at low reynolds numbers[END_REF][START_REF] Dj Cleaver | Bifurcating flows of plunging aerofoils at high strouhal numbers[END_REF][START_REF] Dj Cleaver | Investigation of high-lift mechanisms for a flat-plate airfoil undergoing small-amplitude plunging oscillations[END_REF]. Using direct measurements of the aerodynamics forces, they reported that the deviation was associated with high lift production. They also showed the deviation occurs through a supercritical bifurcation by obtaining the same results for increasing or decreasing flapping frequencies without discontinuity.

It should be noted that symmetry-breaking bifurcations also exist in the related problem of oscillating objects in a quiescent fluid. Experimental [START_REF] Vandenberghe | On unidirectional flight of a free flapping wing[END_REF] and numerical [START_REF] Alben | Coherent locomotion as an attracting state for a free flapping body[END_REF][START_REF] Zhang | Effects of geometric shape on the hydrodynamics of a self-propelled flapping foil[END_REF] observations have shown that bodies of different shape oscillating vertically in a symmetrical way can create a non-symmetric flow, leading to a non-zero mean horizontal force. If the object is free to move in the horizontal direction, this force can set the object into motion with a non-zero mean velocity in either the left or right direction. For the simplest case of a two-dimensional cylinder, Elston et al. [START_REF] John R Elston | Two-dimensional floquet stability analysis of the flow produced by an oscillating circular cylinder in quiescent fluid[END_REF][START_REF] John R Elston | The primary and secondary instabilities of flow generated by an oscillating circular cylinder[END_REF] performed a Floquet stability analysis of the periodic, symmetric flow. They showed that the symmetry-breaking occurring in the flow is effectively due to the onset of a linear mode which breaks the spatial, left-right mirror symmetry of the base flow.

In the present paper, it is proposed to apply a similar stability approach to explain the deviation of the wake of a flapping foil. Unlike in the case of objects oscillating in a quiescent fluid discussed in the previous paragraph, the symmetry which is broken by the bifurcation is not only spatial but spatio-temporal. Namely, in the non-deviated Bénard-Von Karman vortex street before the bifurcation, the half-period corresponding to a downward stroke is the mirror image of the half-period corresponding to an upward stroke. Thus, a difficulty lies in the computation of this time-periodic base flow in unstable cases. A review of methods to compute unstable base flows highlights several approaches. In the case of spatial symmetries as in the work of Elston et al., one can impose the appropriate boundary conditions on the symmetry axis. However, in the case of spatio-temporal symmetry, the corresponding boundary conditions are unknown. When the base flow is stationary, one can use the selective frequency damping (SFD) method [START_REF] Åkervik | Steady solutions of the navier-stokes equations by selective frequency damping[END_REF] to filter unstable temporal frequencies. None of these approaches is directly applicable to the present case. An innovative method was therefore designed for this purpose. The method is related to the SFD in the sense it is used to damp the component of the time-periodic base flow that breaks the spatio-temporal symmetry.

The objective of the present study is to investigate numerically the deviation of the propulsive wake behind a pitching wing, not only by using two-dimensional unsteady non-linear simulations, but also by performing linear stability analysis of the time-periodic nondeviated wake flow. Section II describes the flow con-figuration and results of unsteady non-linear simulations. They are performed for fixed flapping amplitude and a large range of flapping frequencies, so as to identify non-deviated, weakly deviated and strongly deviated wake-flow regimes. The linear stability analysis of the propulsive wake flow is presented in section III. Section III A introduces the new "flow symmetry preserving" method used to compute the time-periodic non-deviated flow above the critical frequency. Section III B presents the Floquet stability analysis of the nondeviated wake-flows and details linear results obtained around the critical flapping frequency. Section IV investigates the non-linear evolution of time-periodic perturbations. A decomposition of the non-linear perturbation into its spatio-temporal symmetric and anti-symmetric components is introduced in this section to understand better the occurrence of strongly deviated wakes for the largest flapping frequencies.

II. UNSTEADY NON-LINEAR SIMULATIONS

The flow around a foil immersed in a fluid of viscosity ν with an incoming uniform velocity U 0 is investigated. The geometry of the foil, shown in Figure 1, is similar to the one used in experimental studies [START_REF] Godoy-Diana | Transitions in the wake of a flapping foil[END_REF][START_REF] Godoy-Diana | A model for the symmetry breaking of the reverse bénard-von kármán vortex street produced by a flapping foil[END_REF]. The leading edge of the foil is a half-cylinder of diameter D, its central part has the form of a triangle which is closed at the trailing edge of the foil by a smaller half-cylinder of diameter d. The ratio of the cylinder diameters is fixed to D/d = 20 and the chord-to-diameter aspect ratio is c/D = 4. Hereafter, all the variables are made non-dimensional using D and U 0 as characteristic length and velocity. As the flow is incompressible, the Reynolds number Re = U 0 D/ν is the unique flow control parameter, which is fixed to Re = 255 for the entire study. A periodic pitching motion of the foil is imposed along the z-axis located at the center of the large half-cylinder. The imposed flapping rotation follows the sinusoidal law

θ(t) = θ m sin(2πf t) (1)
where f is the non-dimensional flapping frequency, T = 1/f the flapping period and θ m is the maximal rotational angle. The non-dimensional (peak-to-peak) flapping amplitude A = 2(c -D/2)sin(θ m ) = 1.07 is fixed accordingly. In the following, the flapping frequency f will be varied in the range 0 ≤ f ≤ 0.5.

A. Numerical model and method

Following Mougin et al. [START_REF] Mougin | The generalized kirchhoff equations and their application to the interaction between a rigid body and an arbitrary timedependent viscous flow[END_REF], the Navier-Stokes equations are written in a non-inertial frame of reference (e X , e Y ), depicted in Figure 1, which rotates at speed ω = dθ dt around the z-axis but does not translate in the laboratory frame of reference (e x , e y ). The vector field u = (u, v) T then represents the flow velocity in written (e X , e Y ) and the incompressible Navier-Stokes equations are written

∂ t u = R (u, p ; θ(t)) , ∇ • u = 0 (2)
with a right-hand-side operator defined by

R (u, p ; θ(t)) = -ω(t) e z × u -((u -w) • ∇) u -∇p + 1 Re ∆u (3)
Note that the spatial operators are defined with respect to the spatial coordinates X = (X, Y ). Compared with the Navier-Stokes equations written in an inertial frame of reference, two additional terms appear in the righthand side: the rotational acceleration ω(t) e z × u and a modification of the convective velocity by the velocity vector w which is defined by

w(X, θ) = -(cos(θ)e X -sin(θ)e Y ) + ω(t) e z × X (4)
The first term accounts for the translational velocity of the foil, written in the rotating frame of reference, while the second term accounts for the angular velocity of the foil. At the surface of the foil, indicated by X w , the flow velocity is equal to the velocity of the foil, i.e.

u(X w ) = w(X w , θ)

Note that the translational velocity is imposed at the surface of the foil, not at the inlet of the computational domain, since the rotating frame of reference does not translate with the foil. Consequently, the boundary condition u = 0 is imposed at the inlet of the computational domain. The no-stress boundary condition (-pn + 1/Re(∇u + ∇u T )n = 0) is used at the outlet. The inlet and outlet of the computational domain are indicated in Figure 2.

The spatial discretisation of the Navier-Stokes equations [START_REF] Sanjay | The aerodynamic effects of wing rotation and a revised quasi-steady model of flapping flight[END_REF][START_REF] Lentink | The scalable design of flapping micro-air vehicles inspired by insect flight[END_REF] is based on a Finite Element formulation [START_REF] Hecht | New development in freefem++[END_REF], using quadratic elements (P 2 ) for the velocity u and linear elements (P 1 ) for the pressure p. A semi-implicit second-order accurate temporal discretisation is used: the temporal derivative is approximated by a backwarddifferential formulae, the linear term in the right-handside (3) are implicit and the non-linear terms are explicit, using a second-order accurate extrapolation. This semi-implicit discretisation yields a non-homogeneous unsteady Stokes problem, which is efficiently solved at each iteration with a preconditioned conjugate gradient algorithm [START_REF] Cahouet | Some fast 3d finite element solvers for the generalized stokes problem[END_REF][START_REF] Glowinski | Finite element methods for incompressible viscous flow[END_REF]. This numerical method, implemented in the non-commercial software FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF], is validated in Appendix §A.

The computational domain and mesh used to discre-

inlet outlet FIG.
2. Computational domain and typical mesh used for the finite-element discretisation. The mesh is made of triangles, only one tenth of them are shown in the Figure .   tise the Navier-Stokes equations are displayed in Figure 2. The influence of the size of the computational domain on the numerical results is described in Appendix §A.

The selected domain is a circle of diameter 60, centered at the leading-edge half-cylinder of the foil. The mesh is composed of 42 • 10 3 triangles, yielding 170 • 10 3 degrees of freedom for each component of the velocity and 42 • 10 3 degrees of freedom for the pressure. A particular attention has been paid to create a symmetric mesh with respect to the x-axis, to avoid introducing artificial asymmetries in the flow. The size and spatial distribution of the triangles is adapted to the flow field, leading to a refinement of the mesh around the foil and in its wake. The smallest triangles are located at the trailing-edge of the foil and their size is 7.29 • 10 -3 . Recalling that the wake flow rotates periodically around the X-axis (while the foil is fixed), the far-wake region is refined accordingly, as shown in Figure 2. Finally, the semi-implicit time discretisation described above implies a Courant-Friedrich-Levy stability condition [START_REF] Kress | Time step restrictions using semi-explicit methods for the incompressible navier-stokes equations[END_REF]. The time-step is chosen accordingly , from ∆t = 8 • 10 -3 for f = 0 down to ∆t = 4 • 10 -4 corresponding to ∆t = T /5000 for large frequencies above the deviation.

B. Results

Non-linear simulations are performed for flapping frequencies in the range 0 ≤ f ≤ 0.5. For f = 0, the foil is in a fixed position with zero angle of incidence with respect to the incoming uniform flow. The wake flow pattern, computed for Re = 255 but not displayed here, is a classical Bénard-Von Karman vortex street, characterised by a (non-dimensional) natural frequency f 0 = 0.167. Let us now investigate the flow pattern obtained when the foil is forced to pitch. Snapshots of the vorticity field, obtained for four flapping frequencies, are displayed in Figure 3 at times t = T /4 (left) and t = 3T /4 (right), that correspond to the maximal and minimal angular positions of the foil, respectively. Corresponding movies are available as Supplemental Materials [37]. In these Figures, the vorticity is displayed in the laboratory frame of reference (e x , e y ). The four flapping frequencies used in Figure 3 correspond to four specific patterns observed in the wake. The flow frequency is always equal to the flapping frequency of the foil.

For the lowest flapping frequency f = 0.1 (see Figures 3 a-b), the wake flow pattern is similar to the classical Bénard-Von Karman vortex street obtained without pitching. Negative (black) and positive (white) vortices are alternately shed during the upstroke and downstroke phases of the flapping motion. In the far-wake, the core of negative and positive vortices lies above and under the x-axis, respectively. When the flapping frequency is equal to the natural flow frequency, i.e. f = f 0 , no vortex lock-in behaviour is observed. When further increasing the flapping frequency, the spatio-temporal structure of the wake flow changes dramatically.

For the frequency f = 0.35 (see Figures 3 c-d), negative and positive vortices of larger magnitude are shed during the upstroke and downstroke phase of the flapping motion. Unlike for the classic Bénard-Von Karman vortex steet, the cores of the negative and positive vortices now lie below and above the x-axis, respectively. This is the characteristic pattern of the reverse Bénard-Von Karman vortex street [START_REF] Jones | Wake structures behind plunging airfoils: a comparison of numerical and experimental results[END_REF]. Interestingly, the vorticity fields satisfy the spatio-temporal symmetry ω z (x, y, t) = -ω z (x, -y, t + T /2), at any position (x, y) downstream of the wing. Around the wing, this relation is not valid as some positions (x, y) belong to the flow or to the solid depending on the flapping motion. However, when considering the spatial position written in the rotating frame of reference, the spatio-temporal symmetry

(u, v, p)(X, Y, t) = (u, -v, p)(X, -Y, t + T /2) ω z (X, Y, t) = -ω z (X, -Y, t + T /2) (6)
holds for all points (X, Y ). This spatio-temporal symmetry breaks when the flapping frequency is increased. For the frequency f = 0.43 shown in Figures 3(e-f) the wake exhibits an upward deviation. Positive and negative vortices are still alternately shed from the trailing edge but they are now convected downstream along a direction not aligned with the streamwise axis y = 0.

For the frequency f = 0.45, the deviation is even more pronounced, as seen in Figures 3(g-h). Positive and negative vortices, instead of being arranged in an array of monopolar structures as for lower frequencies, now form an array of dipolar structures [START_REF] Godoy-Diana | A model for the symmetry breaking of the reverse bénard-von kármán vortex street produced by a flapping foil[END_REF].

Most of the results presented in this study exhibit an upward deviation. However, the downward deviation is also possible and illustrated in Figure 4(b) next to the upward deviation in Figure 4(a). In both cases, the spatio-temporal symmetry ( 6) is broken. During our simulations, upward or downward deviations were obtained by starting the pitching of the wing to the top or the bottom respectively. The selection of the deviation's direction seems to depend on the initial condition but was not investigate in our study. More details can be found on this matter in [START_REF] Lai | Jet characteristics of a plunging airfoil[END_REF][START_REF] Ellenrieder | Piv measurements of the asymmetric wake of a two dimensional heaving hydrofoil[END_REF][START_REF] Zc Zheng | Study of mechanisms and factors that influence the formation of vortical wake of a heaving airfoil[END_REF].

The deviation and antisymmetry of the wake flow can be analysed further by performing a temporal Fourier decomposition of the periodic vorticity. This decomposition writes

ω z (x, t) = ω z (x) + k≥1 (ω z k (x)e ik2πt/T + c.c.) (7)
where ω z is the (real) mean vorticity and ω z k is the k th (complex) Fourier component. In particular, the real and imaginary parts of the first Fourier component ω z 1 extract the component of the flow which oscillates at the flapping frequency f = 1/T . As long as the flow solution remains T -periodic, the Fourier components are easily computed as Positive and negative values are observed in the lower and upper parts of the wake, respectively. This is typical of dragging wake flow profiles, with a negative upper shear layer and a positive lower shear layer. For f = 0.35 in Figure 5(c), the mean vorticity field remains anti-symmetric with respect to the x-axis. However, negative and positive values are now observed in the lower and upper parts of the wake, respectively. This is typical of mean jet flow profiles, with a positive upper shear layer and a negative lower shear layer. Around the foil, the signs of the shear layers are opposite : a negative upper shear layer and a positive lower shear layer. The inversion of the sign of the two shear layers occurs around the trailing edge and is characteristic of propulsive flapping bodies [START_REF] Lai | Jet characteristics of a plunging airfoil[END_REF][START_REF] Yu | Experimental and numerical investigations on the asymmetric wake vortex structures around an oscillating airfoil[END_REF]. For both these frequencies, the vorticity of the first Fourier component in 5(b,d) is symmetric with respect to the x-axis, i.e. ω z 1 (x, -y) = ω z 1 (x, -y) unlike the mean vorticity. The vorticity of the first Fourier component exhibits its strongest values downstream of the foil.

ω z k (x) = 1 T t0+T t0 ω z (x, t)e -ik2πt/T dt (8)
In the f = 0.43 case (Figures 5(e-f)), the mean vorticity field is stronger in magnitude and clearly deviated upward. This deviation is also visible on the first Fourier component and on higher-order Fourier components, not shown here. The spatial anti-symmetry of the mean vorticity field as well as the symmetry of the vorticity of the first Fourier components are broken because of this lateral deviation of the wake. It is even more so the case for f = 0.45 (Figures 5(g-h)). The deviation of the wake- flow described above impacts the aerodynamic forces exerted on the foil. The instantaneous drag and lift coefficients are defined as C x (t) = 2F x (t) and C y (t) = 2F y (t), with F x and F y the components of the non-dimensional aerodynamic force. When the flow satisfies the spatiotemporal symmetry (6), the drag and lift coefficients satisfy

C x (t + T /2) = C x (t) , C y (t + T /2) = -C y (t) (9)
i.e. the instantaneous drag coefficients are equal in the upstroke and downstroke phases of the foil, while the instantaneous lift coefficients are of opposite sign. This property is proved in Appendix C. Figure 6 displays the time-evolution of the system in a C x (t)/C y (t) phase diagram. The trajectory is a closed orbit, which is consistent with the fact that the flow is always periodic. Note that the lift coefficient is one order of magnitude larger than the drag coefficient as previously reported [START_REF] Martín-Alcántara | Vortex flow structures and interactions for the optimum thrust efficiency of a heaving airfoil at different mean angles of attack[END_REF]. For f = 0.35 (solid line) the trajectory has a butterfly shape, symmetric with respect to the C y = 0 axis, implying that the symmetry ( 9) is verified. Regarding the f = 0.43 case (dashed line), first we notice that during the cycle, both the lift and drag reach larger values. Most importantly, the butterfly shape of the trajectory in the (C x , C y )-plane is no longer symmetric with respect to the C y = 0 axis, implying that the symmetry ( 9) is broken. This symmetry breaking is even more pronounced for the f = 0.45 case (dash-dotted line).

Along with the instantaneous values of the lift and drag coefficients, Figure 6 also displays their values averaged over one cycle. For all three cases displayed, the mean drag coefficient is negative, implying that the flapping foil is in the propulsive regime. The mean lift is strictly zero for f = 0.35 and the symmetry relations (9) are satisfied. On the other hand, for f = 0.43 and f = 0.45, the spatio-temporal symmetry is broken, resulting in non-zero mean lift. Interestingly, the mean lift is found to be positive (resp. negative) when the wake is deviated upward (resp. downward), in agreement with results of previous studies [START_REF] James Cleaver | Lift enhancement by means of smallamplitude airfoil oscillations at low reynolds numbers[END_REF][START_REF] De Calderon | On the absence of asymmetric wakes for periodically plunging finite wings[END_REF][START_REF] Emblemsvag | Numerical simulation of flapping micro air vehicles[END_REF]. As explained by Cleaver et al. [START_REF] James Cleaver | Lift enhancement by means of smallamplitude airfoil oscillations at low reynolds numbers[END_REF], the sign of the mean lift is the same as the deviation direction because of a low (mean) pressure region existing on the same side of the foil as the deviation direction. This low-pressure region is induced by vortices detaching from the leading edge and travelling till the trailing edge, as visible on the top of the foil in Figure 3(g,h). The pressure lift and drag being two orders of magnitude larger than their respective viscous components, this low-pression region indeed explains why the lift force is oriented on the same side as the wake is deviated.

The evolution of the mean lift and drag coefficients is examined in Figure 7 as a function of the flapping frequency. Three successive transitions are identified in the range 0.3 ≤ f ≤ 0.5 . First there is a transition from a drag regime to a thrust regime, which occurs at f ∼ 0.31 according to Figure 7(b). This first transition was examined in [START_REF] Jones | Experimental and computational investigation of the knoller-betz effect[END_REF][START_REF] Godoy-Diana | Transitions in the wake of a flapping foil[END_REF] and is not the object of the present paper. The second transition is the spatio-temporal symmetry-breaking leading to non-zero C y occurring for f ∼ 0.39. This transition separates the regimes I and II in Figure 7, which correspond to non-deviated wakes and deviated wakes respectively. As already stated, both upward and downward deviations are equally possible and correspond to positive and negative C y respectively. Both cases are displayed in Figure 7(a) which strongly suggests that the onset of the wake deviation is due to a supercritical bifurcation of the periodic flow. This will be examined in section III by performing a linear Floquet analysis of the symmetric periodic flow.

A third transition is eventually observed in Figure 7 for f 0.435 where the mean lift reaches much higher values and the mean drag experiences a change of slope. This last transition separates the regimes II and III in Figure 7, the latter corresponding to strongly deviated wakes. This transition will be examined in section IV where it is argued that it results from non-linear effects.

III. FLOW BIFURCATION AND DEVIATION

To investigate the deviation of the periodic wake-flow as a bifurcation problem, it is first required to compute the symmetry-preserving base flow solution above the bifurcation threshold. A novel method is introduced in §III A which aims at preserving the spatio-temporal symmetry of the flow field to compute this unstable base flow. The linear stability of this symmetric base flow is then studied in §III B by performing a Floquet stability analysis.

A. Symmetry-preserving method

To present the method, let us first introduce the operator S which acts on a velocity field u and extracts the component satisfying the spatio-temporal symmetry [START_REF] Leontini | Three-dimensional transition in the wake of a transversely oscillating cylinder[END_REF]. This operator is defined by

u s = Su = 1 2 u(x, y, t) + u(x, -y, t -T /2) v(x, y, t) -v(x, -y, t -T /2) ( 10
)
where u s denotes the symmetric component. Similarly, the anti-symmetric component u a (with respect to the spatio-temporal symmetry ( 6)) can be extracted using the operator A defined by

u a = Au = 1 2 u(x, y, t) -u(x, -y, t -T /2) v(x, y, t) + v(x, -y, t -T /2) (11)
Then the flow can then be decomposed as

u = u s + u a ( 12
)
Note that the symmetric and anti-symmetric components satisfy Su s = u s ; Su a = 0 ; Au a = u a ; Au s = 0.

Inserting the above decomposition into the Navier-Stokes equations ( 2),( 3) and (5) yields the system of coupled equations governing the dynamics of the symmetric and anti-symmetric components

∂u s ∂t = R(u s , p s ; θ(t)) -(u a • ∇) u a ∂u a ∂t = L(u s ; θ(t))u a (13)
supplemented with the wall boundary conditions

u s (X w , t) = w(X w , θ(t)) , u a (X w , t) = 0 (14)
The non-linear equation governing the symmetric component is the original equation forced by an additional non-linear advection term. The anti-symmetric component is governed by a linear equation, where the linear operator L is defined as

L(u s ; θ)u a = -( dθ dt e z ) × u a -((u s -w(X, θ)) • ∇) u a -(u a • ∇) u s -∇p a + 1 Re ∆u a ( 15
)
and is the linearisation operator of the right-hand-side operator R around the symmetric velocity component u s . Interestingly, if the anti-symmetric component vanishes (i.e. u a = 0), the symmetric component satisfies the original equation ( 2). The symmetry-preserving method consists in driving the anti-symmetric component to zero, by adding a damping term to the second equation ( 13). The new system of equations writes

∂u s ∂t = R(u s , p s ; θ(t)) -(u a • ∇) u a ∂u a ∂t = L(u s ; θ(t))u a -χu a ( 16
)
where χ is a damping parameter. The theoretical choice of this parameter depends on the eigenvalues of largest growth rate of the linear operator L. χ should be chosen so as to stabilise all of these eigenvalues. The eigenvalues of (L -χI) are then all of negative growth rate and the anti-symmetric component is driven to zero. In the following, the solution satisfying the spatio-temporal symmetry (also called the base flow) is denoted U s (U a = 0) and satisfies

∂U s ∂t = R(U s , P s ; θ(t)) (17)
U s (X w , t) = w(X w , θ(t))

In practice, the eigenvalues of largest growth rate are not known and a trial-and-error method is used to determine values of the parameter χ. Its effect is discussed in Appendix B.

For the frequency f = 0.43, the solution of the original equation (2,5) is a deviated wake as seen in Figure 3. Starting from this solution, the new system of equations ( 16) is solved with the damping parameter χ = 0.3. As a consequence, the antisymmetric component vanishes, as shown in Figure 8 where its norm is given as a function of time. Snapshots and first Fourier components of the solution are displayed in Figure 9. A movie of the non-deviated wake for f = 0.43 is available as a Supplemental Material [37]. The spatio-temporal symmetry of the solution is clearly visible in the two snapshots of the vorticity separated by T /2. The suppression of the deviation is clearly observed on the mean vorticity field and on the first Fourier component. Using this symmetry-preserving method for several frequencies, the mean drag and lift coefficients of the solutions are computed and shown in Figure 10 with black circles. Results obtained without the symmetry preserving method, already shown in Figure 7, are recalled in the Figure with white circles. The left plot displaying C y demonstrates that the symmetrypreserving method is efficient to compute a symmetric base flow even in the range f 0.39 where it is unstable and cannot be accessed using temporal integration of the starting equations which only give access to the asymmetric states. The right plot displaying C x shows that for 0.39 f 0.43 the symmetry-breaking has almost no impact on the drag (the open circles being superposed to the black ones), while it has a strong impact above the transition for f 0.44 which will be examined in section IV.

For f ≥ 0.39, three branches of periodic solutions exist: one branch of symmetric (or non-deviated) periodic solutions characterised by zero mean lift and two branches of asymmetric (or deviated) periodic solutions characterised by mean lift of opposite sign. To complete this bifurcation diagram, the stability of the symmetric solutions is addressed in the next paragraph using a Floquet stability analysis.

B. Floquet stability analysis of the non-deviated wake flow

The stability of the time-periodic base flow U s is determined by investigating the long-term dynamics of infinitesimal perturbations u . The flow field is first decomposed as u = U s + εu with ε 1 and this decomposition is introduced into the Navier-Stokes equations [START_REF] Sanjay | The aerodynamic effects of wing rotation and a revised quasi-steady model of flapping flight[END_REF]. Recalling that the base flow is governed by [START_REF] Ellenrieder | Piv measurements of the asymmetric wake of a two dimensional heaving hydrofoil[END_REF] and neglecting the quadratic terms in ε, we obtain the equation governing the dynamics of linear perturbations ∂u ∂t = L(U s )u [START_REF] Zc Zheng | Study of mechanisms and factors that influence the formation of vortical wake of a heaving airfoil[END_REF] with the wall boundary condition u (X w , t) = 0. Note that symmetric perturbations u s = Su and antisymmetric perturbations u a = Au are both solutions of the above equation. The linear operator L(U s ) is Tperiodic because the base flow U s is T -periodic. Thus, according to Floquet theory [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF], any solution of ( 18) can be decomposed into the sum

u (x, t) = k ûk (x, t)e λ k t (19)
where The numerical method used to compute the Floquet modes and multipliers is similar to that used in [START_REF] Chowdhury | Asymmetry and hopf bifurcation in spherical couette flow[END_REF] and [START_REF] Barkley | Threedimensional floquet stability analysis of the wake of a circular cylinder[END_REF]. The evolution of the linear perturbation u over one period T is formally rewritten

u (t 0 + T ) = Pu (t 0 ) ( 20
)
where P is the propagator over one period, also known as the linearised Poincaré map. The action of this propagator on the perturbation u (t 0 ) at the arbitrary time t 0 is obtained by integrating the linearised equation ( 18) in time from t 0 to t 0 + T . The eigenvalues of the propagator P are precisely the Floquet multipliers µ k of L, and the eigenvectors of P correspond to the Floquet modes û(x, t 0 ) of L for the arbitrary time t 0 . The time-periodic evolution of the Floquet mode is then determined by temporal integration of equation ( 18) over one period starting with û(x, t 0 ) as the initial condition. The Floquet mode is the solution of this temporal integration, corrected by the multiplicative factor e -λ k t to account for the instantaneous growth or decay. Arnoldi method was used to compute the eigenvalues of largest amplitude [START_REF] Saad | Numerical methods for large eigenvalue problems[END_REF]. A serial implementation using the modified Gram-Schmidt algorithm for the orthogonalization process is used to generate an approximation of n eigenvectors. All computed modes were normalised by their total kinetic energy.

First, results of the stability analysis are shown for the flapping frequency f = 0.4. The largest Floquet multipliers are displayed in Figure 11(a). Most of the eigenvalues are stable, but one eigenvalue marked by the black circle is unstable as it lies outside the unit circle. As this eigenvalue is real (µ 1 = 1.046 + 0.0i), the Floquet mode is synchronous with the base flow and therefore does not modify its periodicity. The evolution of this eigenvalue with the flapping frequency is shown in Figure 11(b). The Floquet mode becomes unstable for the critical frequency f c ∼ 0.385. This critical frequency is consistent with the observation of wake deviation and the apparition of mean lift in Figure 7 which define the transition between regimes I and II. For every flapping frequency tested, the leading Floquet mode was found to be real and therefore synchronous with the base flow. No other unstable modes were detected for frequencies below f = 0.46. It can also be seen that w z (X, Y, T /4) = w z (X, -Y, 3T /4). This property has been verified for all t over one flapping period and demonstrates that the leading Floquet mode breaks the spatio-temporal symmetry [START_REF] Leontini | Three-dimensional transition in the wake of a transversely oscillating cylinder[END_REF] and is therefore of the form u a . To explain the link between the computed Flo- quet mode and the deviation of the wake, we provide, in Figure 13, a snapshot of (a) the base flow, (b) the leading Floquet mode and (c) their superimposition inside the dotted square (the base flow is represented with isolines and the mode with isocontours). First, a certain synchronisation between the two flows can be observed.

In the whole wake and at all times, a vortex of the base flow is synchronised with a dipole of the linear mode.

A dipolar perturbation of a monopolar vortex is a classic structure called a displacement mode [START_REF] Leibovich | Bending waves on inviscid columnar vortices[END_REF][START_REF] Fabre | Kelvin waves and the singular modes of the lamb-oseen vortex[END_REF][START_REF] Vincent Brion | Linear dynamics of the lamb-chaplygin dipole in the twodimensional limit[END_REF]. It is well known that this structure is associated with a displacement of the vorticity centroid. More specifically, consider the negative vorticity monopole located at the left of the dotted square. The corresponding dipole in the Floquet mode has negative vorticity (black) in the upper part and positive vorticity (white) in the lower part. The superposition strengthens the top part of the monopole and weakens its lower part, resulting in a net displacement of the monopole in the upward direction. The same argument can be applied to the positive vorticity monopole located at the right of the dotted square, which is also displaced upward. Thus, the structure of the Floquet mode is able to explain the deviation of the whole wake in the upward direction. Note that as the amplitude of the linear mode is arbitrary, changing its sign leads to the equally probable deviation of the wake in the downward direction. The structure of the Floquet mode can also be analysed using temporal Fourier decomposition, just as was done for the flow obtained through non-linear simulations in §II. Figure 14 displays the mean vorticity (a) and the first Fourier component (b) which are respectively symmetric and anti-symmetric with respect to the x-axis. This is the opposite of the base flow presented in Figure 9 confirming that the mode does break the symmetry [START_REF] Leontini | Three-dimensional transition in the wake of a transversely oscillating cylinder[END_REF]. Three alternate vorticity layers are present in the mean wake, with negative vorticity on the x-axis. By superimposing it with the mean wake of the base flow (in Figure 9c), and using the same arguments as for the instantaneous superimposition, the mean wake of the mode will strengthen both the positive upper and negative lower shear layers on their tops, and weaken their bottoms. This will result in a deviated jet-like mean wake, like the bottom-left frame of Figure 5.

IV. NON-LINEAR EFFECTS

The linear stability analysis detailed in the previous section explains the onset of the small deviation observed in regime II (0.385 f 0.43). However it does not provide any explanation for the large deviation, the sudden increase of thrust and lift obtained in regime III (f 0.44). To better understand the large deviation and corresponding increase of thrust/lift obtained in regime III, the effects of non-linearity are now analysed. The non-linear perturbation is defined as the difference between the asymmetric flow and the symmetric base flow, i.e. u = u -U s . Linear and non-linear perturbations are first compared for frequencies f = 0.43 and f = 0.45 which correspond to regime II and III respectively.

For f = 0.43, snapshots of the linear and non-linear perturbations are shown in Figures 15(a) and 15(c). The linear perturbation is an array of dipolar structures aligned with the x-axis and centered in the monopolar structures of the base flow, as described previously. These dipolar structures grow in time, and in space when moving downstream. The spatial growth is quantified by examining Figure 15(e), which shows the streamwise evolution of the kinetic perturbation energy integrated in y and averaged over one flapping period. The spatial growth of the linear perturbation (dashed line) is quasiexponential. Note that its value depend on the normalisation of the Floquet mode (unitary total kinetic energy in our study). In order to be the same order of magnitude as the non-linear perturbation, the y-integrated kinetic energy of the linear perturbation has been mutliplied by 10 4 . Regarding the pattern of the non-linear perturbation shown in Figure 15(c), we observe that it is very similar to that of the linear perturbation in the near-wake. The spatial growth of the linear and non-linear perturbations are also very similar for x ≤ 10 (Figure 15(e)). In this region, the effect of non-linearities is therefore negligible. However, in the far-wake (for x ≥ 15), the kinetic energy of the non-linear perturbation tends towards a constant value, unlike the linear perturbation which still grows in space. Clearly the effect of non-linearities is to saturate the spatial growth of the perturbation. This non-linear saturation comes with a very specific pattern in the far-wake region (see Figure 15(c)). Two arrays of vortical structures are now visible. One array of monopolar structures remains aligned with the x-axis while the other is slightly deviated upward. In this far-wake region, the monopolar structures in the non-deviated array are of opposite sign to the monopolar structures of the base flow (compare Figure 15(c) with 9). Therefore, when adding the non-linear perturbation to the base flow, the array of monopolar structures disappears in the far-wake and only the deviated array of dipolar structures remains in the instantaneous solution. As such, a second effect of non-linearities is to increase the deviation of the wake. Let us now examine the non-linear effects for the larger frequency f = 0.45 in regime III where large deviation and thrust are obtained. Snapshots of the linear and non-linear perturbations are shown in Figures 15(b) and 15(d), and the streamwise evolution of the kinetic energy is displayed in Figure 15(f). The region where the linear and non-linear spatial growths match, is limited to the very near-wake region. For x ≥ 7, the saturation occurs and the two arrays of vortical structures are visible in the wake. The array of dipolar structures is more deviated than for f = 0.43 and the dipolar structure is more pronounced.

When increasing the flapping frequency, non-linear saturation occurs closer to the wing. This may impact the aerodynamic forces exerted on the wing and explain the increase of thrust observed in regime III. First, we note that a linear anti-symmetric perturbation generates a mean lift force but does not create a mean drag/thrust force (see details in Appendix §C). Therefore, the additional thrust observed for f ≥ 0.44 is necessarily associated to symmetric flow components. To better understand how non-linear terms produce a symmetric flow component, the equations governing the dynamics of the symmetric and anti-symmetric components are examined. To that aim, the non-linear perturbation is decomposed into its symmetric u s and anti-symmetric u a components. They are defined respectively by Su s = u s and Su a = 0. Introducing this decomposition into (13) yields the following system of coupled non-linear equations

             ∂u s ∂t = L(U s ; θ(t))u s -[u s • ∇] u s -[u a • ∇] u a ∂u a ∂t = L(U s ; θ(t))u a -[u a • ∇] u s -[u s • ∇] u a (21)
where the first and second equations govern the nonlinear dynamics of the symmetric and anti-symmetric components, respectively. In addition to the linear operator L around the periodic base solution U s , quadratic non-linear terms appear in the right-hand-sides and couple the two dynamics. In the linear regime, for which all quadratic terms are negligible, the symmetric component decays in time while the anti-symmetric component grows in time. For large enough amplitude of the anti-symmetric component, its quadratic interaction with itself in the first equation drives the development of a non-zero symmetric component. The existence of such a symmetric component influences in return the anti-symmetric dynamics, via the quadratic interactions between u s and u a in the second equation of [START_REF] Marais | Stabilizing effect of flexibility in the wake of a flapping foil[END_REF]. Let us examine now the spatial evolution of the symmetric and anti-symmetric components of the non-linear perturbation computed for the flapping frequencies f = 0.43 (regime II) and f = 0.45 (regime III). The mean kinetic energy integrated in the y direction is displayed in Figure 16 for the non-linear perturbation (black solid line) and its symmetric (red dashed line) and anti-symmetric (blue dash-dotted line) components. Results obtained for f = 0.43 and f = 0.45 are shown in Figures 16(a) and 16(b) respectively and present similar tendencies when progressing downstream. The anti-symmetric perturbation is amplified first while the symmetric one is negligible, as can be expected in the linear regime. Then the amplitude of the antisymmetric perturbation saturates while the symmetric one increases. As explained above, the transfer of energy from the anti-symmetric component to the symmetric component is ensured by the quadratic interaction of u a with itself. For a critical streamwise position x c , the amplitude of the symmetric perturbation becomes larger than that of the anti-symmetric one. Further downstream, the anti-symmetric component decreases in amplitude, while the symmetric component continues to increase, then saturates before decreasing. The evolution To our knowledge, this is the first time those two distinct deviated regimes are clearly highlighted. Godoy-Diana's work [START_REF] Godoy-Diana | A model for the symmetry breaking of the reverse bénard-von kármán vortex street produced by a flapping foil[END_REF] also presented a change of slope in the evolution of the angle of deviation with increasing flapping frequency. However, this change of slope was associated with significant error bars, which makes it questionable whether this is the same phenomenon.

Finally, the transition between the two regimes is also related to a change in the wake structure. As can be seen in Figure 3(g,h) for f = 0.45, successive counterrotating vortices approach each other in pairs to form dipolar structures. For lower unstable frequencies such as f = 0.43, this is much less notable and occurs only far from the foil as seen in Figure 3(e,f). This dipolar structure was previously noted by Godoy-Diana et al. [START_REF] Godoy-Diana | A model for the symmetry breaking of the reverse bénard-von kármán vortex street produced by a flapping foil[END_REF]. They derived a quantitative criterion based on the phase velocity produced by two successive counter-rotating vortices to identify the deviation threshold. They showed that the dipolar organisation of the wake promotes the deviation, without ruling on any causal relationship. Later, Zheng et al. [START_REF] Zc Zheng | Study of mechanisms and factors that influence the formation of vortical wake of a heaving airfoil[END_REF] derived a similar model based on the competition between successive phase velocities and attributed the wake deviation to the vortex pairing in the wake. Our study shows the dipolar aggregation develops in the wake, far from the foil in regime II. For increasing flapping frequency, the aggregation is visibly closer to the foil due to non-linear amplification. The clear aggregation of vortices directly after the trailing edge marks the transition with regime III. These observations imply that the dipolar aggregation, though shown in these studies to be favorable to deviation, does not cause it. It is a secondary effect of the deviation, enhancing it through the non-linear coupling of symmetric and anti-symmetric perturbations. Ultimately, this non-linear coupling results in a strong increase of the deviation angle, the mean lift and thrust.

V. CONCLUSION

In this paper, we investigated the deviation of the wake behind a flapping foil. A novel flow symmetry-preserving method was specifically developed to compute the unstable non-deviated wake. This method is based on the decomposition of the governing equations into a set of two equations governing the dynamics of the non-deviated wake and the perturbations evolving around it. Then, the second equation is stabilised to suppress the perturbations that break the spatio-temporal symmetry inherent to the flapping motion. A Floquet stability analysis of the time-periodic non-deviated wake showed the existence of a synchronous anti-symmetric mode becoming unstable at the critical flapping frequency where deviation occurs. This Floquet mode is an array of counterrotating dipoles that act as a succession of displacement modes. Both the instantaneous and the averaged-in-time effects of the mode are to displace the non-deviated wake away from the streamwise direction, resulting in a deviated wake. Finally, non-linearities are associated with a dipolar aggregation of vortices in the wake. As non-linear effects develop closer to the foil with increasing flapping frequency, they ultimately occur at the trailing edge, which results in the transition to a third regime with a much stronger deviation of the wake. This third regime is also associated with a strong increase of mean thrust and lift which are associated respectively with symmetric and anti-symmetric perturbations around the foil. used to describe this case : the Keulegan-Carpenter number KC = U/f D and the Reynolds number Re = U D/ν, where U is the maximum velocity of the cylinder, D is the cylinder diameter and ν is the fluid kinematic viscosity. We reproduce some results obtained experimentally by Dütsch et al. [START_REF] Dütsch | Lowreynolds-number flow around an oscillating circular cylinder at low keulegan-carpenter numbers[END_REF] and numerically by Guilmineau et al. [START_REF] Guilmineau | A numerical simulation of vortex shedding from an oscillating circular cylinder[END_REF] and Hosseinjani et al. [START_REF] Akbar | Numerical simulation of the wake structure and thrust/lift generation of a pitching airfoil at low reynolds number via an immersed boundary method[END_REF], at Re = 100 and KC = 5.

In Figure 19 In a second time, we show a large enough computational domain has been chosen to ensure good values of the aerodynamic forces. On this study case, we test four different domains, which size is specified by their radius R. The same refinement has been used for the four meshes. Three frequencies f = 0.35, f = 0.42 and f = 0.44 were tested, corresponding to the three regimes defined in the paper respectively.

In Table I is shown the mean drag coefficient C x for each domain and frequency. It is barely not affected by the domain size for f = 0.35 and f = 0.42. For f = 0.44, it changes a little more with 3% relative difference between the smallest and largest domains.

In Table II is shown the mean lift coefficient C y for each domain and frequency. For f = 0.35 (regime I), the mean lift is zero for every domain, which is consistent with the fact the wake is not deviated. For f = 0.44 (regime III), the mean lift is practically independent of the size of the domain for R > 30 (0.5% relative difference). For f = 0.42 (regime II), the mean lift is more sensitive to the size of the domain with 3% relative difference for the largest domains (R > 30). However, this magnitude of error does not affect the distinction between the different regimes. The lift coefficients are not distinguishable when they are superimposed in a ( C y ,f ) map as in Figure 7. For these reasons, we decided to use the mesh of size R = 30 since it is large enough to ensure good values of aerodynamic forces and allow us to keep a satisfactory computational cost.

Appendix B: Choice of damping coefficient χ in the symmetry preserving method

As explained in §III A, one needs an appropriate choice of damping coefficient in order to stabilise the second equation of system [START_REF] Lai | Jet characteristics of a plunging airfoil[END_REF].

Figure 21 displays the time evolution of the antisymmetric component's norm for four values of the damping coefficient. All four simulations have been initialised with the same deviated wake at f = 0.45. For χ = 0.02, the norm converges towards T-periodic os- cillations (of very small magnitude with respect to its mean magnitude) indicating the anti-symmetric component is not damped and the wake remains deviated. For the three other values of χ, and after a small transitory stage, the norm decreases at almost constant speed towards the zero-machine value. As can be observed, convergence speed increases with the value of χ. Tests showed that there is a minimum threshold of the damping coefficient above which the anti-symmetric component is suppressed. With this flapping frequency, the threshold is between χ = 0.02 and χ = 0.1. Without the damping term, the linear operator L governs the anti-symmetric component dynamics, which asymptotically, will have a non-zero solution if and only if the real part of its leading eigenvalue is positive. The damping coefficient counteracts this so that if χ is bigger than the leading eigenvalue's real part, the complete linear operator L -χI will only have stable eigenvalues. Therefore, the only asymptotic solution of the equation will be nil.

Let us now assume that the flow field is symmetric. We can introduce ( 6) and (C7) into (C6)

F p x (t) = Γ -p(x, -y, t-T /2) n x (x, -y, t-T /2) dΓ (C8)

By introducing a change of variables z = -y into equation (C8) and using the fact that the foil is symmetric with respect to the X-axis, we find the following spatiotemporal property

F p x (t) = F p x (t -T /2) (C9)
As explained previously, this procedure can be reproduced for the viscous component of the streamwise force and for the transverse force to obtain the spatio-temporal properties in (C1) for any symmetric flow field u s . Finally, the property of zero time-averaged lift for any u s is obtained by integrating the second property of (C1) over one flapping period. The properties for any antisymmetric flow field u a in equation (C3) are obtained in the same way by using the property of anti-symmetry (C2).

FIG. 1 .

 1 FIG.1. Sketch of the pitching foil pivoted at the center of the leading edge half-cylinder.

FIG. 3 .

 3 FIG. 3. Instantaneous vorticity fields for the flapping frequency f = 0.1 (a,b), f = 0.35 (c,d), f = 0.43 (e,f) and f = 0.45 (g,h) at t = T /4 (left) and t = 3T /4 (right). Black and white denote negative and positive values respectively.

FIG. 4 .

 4 FIG. 4. Instantaneous vorticity fields for f = 0.43 at t = T /4 showing an upward (a) and downward (b) deviation of the wake. Black and white denote negative and positive values respectively.

Figure 5

 5 Figure 5 displays the mean (left) and first Fourier component (right) of the vorticity field for the flapping frequencies f = 0.1 (a,b), f = 0.35 (c,d), f = 0.43 (e,f) and f = 0.45 (g,h). For f = 0.1 in Figure5(a), the mean vorticity field is anti-symmetric with respect to the x-axis i.e. ω z (x, -y) = -ω z (x, -y). Positive and negative values are observed in the lower and upper parts of the wake, respectively. This is typical of dragging wake flow profiles, with a negative upper shear layer and a positive lower shear layer. For f = 0.35 in Figure5(c), the mean vorticity field remains anti-symmetric with respect to the x-axis. However, negative and positive values are now observed in the lower and upper parts of the wake, respectively. This is typical of mean jet flow profiles, with a positive upper shear layer and a negative lower shear layer. Around the foil, the signs of the shear layers are opposite : a negative upper shear layer and a positive lower shear layer. The inversion of the sign of the two shear layers occurs around the trailing edge and is characteristic of propulsive flapping bodies[START_REF] Lai | Jet characteristics of a plunging airfoil[END_REF][START_REF] Yu | Experimental and numerical investigations on the asymmetric wake vortex structures around an oscillating airfoil[END_REF]. For both these frequencies, the vorticity of the first Fourier component in 5(b,d) is symmetric with respect to the x-axis, i.e. ω z 1 (x, -y) = ω z 1 (x, -y) unlike the mean vorticity. The vorticity of the first Fourier component exhibits its strongest values downstream of the foil.

FIG. 5 .

 5 FIG. 5. Vorticity field of the mean flow (left) and the real part of the first Fourier component (right) for the flapping frequency f = 0.1 (a,b), f = 0.35 (c,d), f = 0.43 (e,f) and f = 0.45 (g,h). Black and white denote negative and positive values respectively. The foil at zero incidence is shown in red.

FIG. 6 .

 6 FIG. 6. Instantaneous drag versus lift coefficients over one flapping period for f = 0.35 (solid line), f = 0.43 (dashed line) and f = 0.45 (dash-dotted line). The mean positions of the drag and lift coefficients are indicated by the black circle f = 0.35), diamond (f = 0.43) and square (f = 0.45).

FIG. 7 .

 7 FIG. 7. Mean lift and drag coefficients as a function of the flapping frequency.

FIG. 8 .

 8 FIG. 8. Evolution of the symmetric (solid) and antisymmetric (dashed) component as a function of time with χ = 0.3 at f = 0.43. Simulation is started from converged deviated wake.

FIG. 9 .

 9 FIG. 9. Vorticity snapshots (top) of the non deviated wake at f = 0.43 obtained using χ = 0.3. Vorticity field averaged in time and first temporal Fourier component (bottom). Black and white denote negative and positive values respectively. The foil at zero incidence is shown in red.

FIG. 10 .

 10 FIG. 10. Mean lift and drag coefficients as functions of the flapping frequency f . •: asymmetric solutions, •: symmetric solutions.

  ûk are T -periodic functions, called the Floquet modes of L, and the complex numbers λ k are the Floquet exponents. The Floquet multipliers µ k = e λ k T correspond to the temporal growth or decay of the Floquet modes over one period T . The stability of the base flow is determined by examining the spectrum of Floquet multipliers. A Floquet mode is stable (resp. unstable) when the corresponding Floquet multiplier lies inside (resp. outside) the unit circle |µ k | < 1 (resp. |µ k | > 1) in the complex plane. When one Floquet mode becomes unstable, the time-periodic base flow becomes unstable.

FIG. 11 .

 11 FIG. 11. Spectrum of the Floquet multipliers at f = 0.4. The dominant mode's Floquet multiplier, in black, is µ1 1.046 (a). Evolution of Floquet multiplier magnitude |µ| with the flapping frequency f (b). The line |µ| = 1 corresponds to marginal stability of Floquet modes. The dominant mode is always real.

Figure 12

 12 Figure 12 displays two snapshots of the dominant Floquet mode's vorticity separated by T /2. A movie of this mode dynamics is available as a Supplemental Material[37]. Vorticity dipoles are emitted at the trailing edge and grow spatially and temporally whilst convected downstream. Negative vortices (in black) are aligned with the central line whereas positive vortices (in white) are alternately above and below the negatives ones. It can also be seen that w z (X, Y, T /4) = w z (X, -Y, 3T /4). This property has been verified for all t over one flapping period and demonstrates that the leading Floquet mode breaks the spatio-temporal symmetry[START_REF] Leontini | Three-dimensional transition in the wake of a transversely oscillating cylinder[END_REF] and is therefore of the form u a . To explain the link between the computed Flo-

FIG. 12 .

 12 FIG. 12. Vorticity snapshots of the dominant Floquet mode at flapping frequency f = 0.45 and time t = T /4 (a) and t = 3T /4 (b). Black and white denote negative and positive values respectively.

FIG. 13 .

 13 FIG. 13. Vorticity snapshots of (a) the base flow, (b) the Floquet mode and (c) a zoom of their superimposition at f = 0.43 and t = 0. Black and white colours (in (a),(b) and (c)), and dotted and continuous lines (in (c)) denote negative and positive values respectively.

FIG. 14 .

 14 FIG. 14. Vorticity field averaged in time over one period and first temporal Fourier component at f = 0.43. Black and white denote negative and positive values respectively. The foil at zero incidence is shown in red.

FIG. 15 .

 15 FIG. 15. (a-d) Snapshots of the vorticity field for (a-b) linear and (c-d) non-linear perturbations at time t = T /4 and frequency (a,c) f = 0.43 and (b,d) f = 0.45. (e-f) y-integrated kinetic energy of the linear (dashed) and non-linear (solid) perturbations as a function of the streamwise coordinate for (e) f = 0.43 and (f) f = 0.45.

FIG. 16 .

 16 FIG. 16. Evolution of the non-linear perturbations' mean kinetic energy (integrated in the transverse direction) in the streamwise direction for f = 0.43 (a) and f = 0.45 (b) : u (black solid), u s (red dashed), u a (blue dash-dotted).

FIG. 18 .

 18 FIG.[START_REF] Zc Zheng | Study of mechanisms and factors that influence the formation of vortical wake of a heaving airfoil[END_REF]. Mean lift and drag coefficients for the symmetric u s (red circles) and anti-symmetric u a (blue diamonds) components as a function of the flapping frequency f .

  we show vorticity isolines, which are compared to the numerical results of Guilmineau et al. at two different moments t = 0 and t = 19T /72 of the oscillation. The results compare very well.

FIG. 19 .

 19 FIG. 19. Vorticity isolines for two different instants t. Guilmineau et al. (2002) at t = 0 (a) and t = 19T /72 (c). Present method at t = 0 (b) and t = 19T /72 (d).

FIG. 20 .

 20 FIG. 20. The x-velocity component in the transverse direction y at x = -0.6D for two different instants t of the flapping period T . (a) t = T /2; (b) t = 7T /12. Dütsch et al. (1998) : black squares. Hosseinjani et al. (2015) : blue dashed line. Present method : black solid line.

  FIG. 21. Norm of the anti-symmetric velocity component as a function of time for increasing damping coefficients χ = [0.02, 0.1, 0.2, 0.3] at flapping frequency f = 0.45.

ACKNOWLEDGMENTS

This work was partially funded by the French Ministry of Defence -Direction Générale de l'Armement.

Appendix A: Validation of the non-linear code and convergence tests

Different tests are used in order to validate the present study. First of all, we validate the numerical method by reproducing the case of a laminar flow generated by an oscillating circular cylinder in an initially quiescent fluid. This case was described and studied experimentally by Dütsch et al. [START_REF] Dütsch | Lowreynolds-number flow around an oscillating circular cylinder at low keulegan-carpenter numbers[END_REF]. The periodic oscillation of the cylinder follows the law x(t) = -Asin(2πf t) where A and f stands for the oscillating amplitude and frequency respectively. Two non-dimensional numbers are Appendix C: Symmetry properties for the forces In this section, more details are given concerning the determination of the following symmetry properties respected by the aerodynamic forces. For a symmetric flow field u s , which respects the spatio-temporal symmetry property (6), the aerodynamic forces verify the following properties

The first two equations imply that the instantaneous drag forces are equal in the upstroke and downstroke phase of the foil, while the instantaneous lift forces are of opposite sign. The third equation states that the time-averaged lift coefficient is zero. The same can be done for an anti-symmetric flow field u a which respects the following symmetry properties instead

In that case, the aerodynamic forces verify the following properties

This time, the instantaneous lift forces are equal in the upstroke and downstroke phase of the foil, while the instantaneous drag forces are of opposite sign. The timeaveraged drag is nil in this case.

In order to show how properties (C1) and (C3) are obtained, we will focus on the first property of (C1) : F x (t -T /2) = F x (t). The others can be derived using the same method. The aerodynamic forces F exerted by the fluid on the foil are expressed as

where Γ is the foil boundary with the fluid, n = (n x , n y ) denotes the unit outward normal. We take its horizontal component

Let us focus on the streamwise pressure forces F p x (t); the following reasoning can be applied the same way to the viscous component F v x (t).

The outward unit normal respects the following symmetry property, as a consequence of the symmetry of the foil w.r.t the X-axis and its periodic motion (n x , n y )(x, y, t) = (n x , -n y )(x, -y, t + T /2) (C7)