Numerical modelling of an undulating membrane tidal energy converter

The undulating membrane tidal energy converter is a device that uses the flutter instabilities occurring from the interaction between a slender body and a fluid flow. A new numerical model has been developed using a 2D corotational finite element method to represent the structure and the unsteady point-vortex method to compute the flow. These methods as well as the interaction process are presented. Trajectory and frequency of the undulating motion, hydrodynamic forces on the structure and velocity field in the wake are presented. Comparison shows a good agreement with experimental results obtained from a 1/20th scale prototype without power take off tested in flume tank.


I. INTRODUCTION
I NTERACTIONS between a semi-rigid membrane and a fluid flow above a certain current speed result in a periodic motion on a given wavemode with an associated frequency. This minimal current speed at which the membrane starts to undulate is called "critical flow speed". The membrane is prestressed by lateral cables that buckle it in order to lower the critical speed. Undulating membrane tidal energy converters use the undulating motion to convert tidal kinetic energy into electricity trough linear generators (Fig. 1). The dynamic deformation of the membrane activate electromechanical linear generators placed in the central line. They slow down the motion while converting deformation energy into electrical potential. Transverse rigid bars ensure a 2D motion and transmit forces from the membrane to the generators.
A fluid-structure interaction numerical model is developed in order to study this device. It follows the development of a 1/20 th scale experimental model [1] and a linear analytical model [2]. The purpose of the numerical model is to represent accurately the membrane motion, the flow around it, and the hydrodynamic forces on the structure. It should also simulate the flow in the wake and must be adaptable to three dimensions for a future side effects study. The fluid model is presented in section II-A. It uses the unsteady point-vortex method [3] to calculate the flow loading on the structure. The position of the membrane is evaluated with the co-rotational finite element method [4] detailed in section II-B. The interaction process is then explained in section II-C. It uses an implicit partitioned approach [5]. Results in term of structure behavior and flow properties are then analyzed in section III. The first comparisons with experimental results from [1] enable to validate these numerical developments.

A. Unsteady point-vortex method
The flow model is based on an unsteady point-vortex method. It uses a rotational formulation of the Navier-Stokes equation, which is discretized and solved in a lagrangian framework. The unsteady point-vortex method is a specific application of the panel method, which have been introduced for the first time in fluid dynamics by Rosenhead [6]. A review of the vortex methods can be found in Leonard [7]. This method gives good results for estimation of hydrodynamic forces and flow field in the wake. Moreover, vortex methods are reputed to be accurate and to require little computation time in 2D. Therefore, it is used in various fluid mechanics domains as animal locomotion modelization or aerodynamic flight [3]. This model is also used in marine renewable energy simulations by [8] and [9].
The velocity-vorticity formulation of the Navier-Stokes equation governing an incompressible fluid is given by Eq. (1), with: ω = ∇ ∧ u the vorticity, ∇ the gradient operator, u the velocity field, t the time and ν the fluid kinematic viscosity. The Helmholtz decomposition of the velocity field is used to decompose the velocity vector into a sum of rotational and irrotational flows (Eq. (2), where φ is the scalar potential and ψ is the vector potential).
The flow field u is decomposed into three components: • A constant flow u ∞ , corresponding to the upstream flow. • A potential flow u φ , representing the influence of the structure on the flow. • A rotational flow u ω in the wake. The vorticity field ω is discretized at time t in singular vorticity elements. Circulation is calculated along the structure, as point-vortices bound to the membrane. Particles carrying vorticity are also emitted at the trailing edge and conveyed in the wake (Fig.  2). They are generated so that Kelvin's theorem is respected. Upper case notation is used in the following for discrete variables and lower case for continuous variables. The velocity induced by a vortex is obtained by the two-dimensional Biot-Savart law. Flow velocity U in a position X = (x, y) is then the sum of the upstream velocity, of the velocity U φ induced by the bound vortices along the membrane and of the velocity U ω induced by the vortices emitted in the wake (Eq. (3)). Γ i is the circulation of the i th bound vortex and X b i is its position, whereas circulation and position of vortices in the wake are noted Γ w i and X w i . N is the number of bound vortices and N w is the number of vortices modelled in the wake. The ⊥ symbol indicates a +π/2 rotation.
The system is solved along the structure, where the normal component of flow speed is equal to the structure speed. This leads to Eq. (4), where X c i is a position on the membrane, the dot represent the time derivative and n i is the local membrane normal unity vector.
The problem is expressed by the matrix linear equation (5). α is the influence coefficient matrix, of which each element represents the influence of vortex j on control point i (Eq. (6), where τ i is the local membrane tangent unity vector). Γ is the bound vortices circulation vector and H the right-hand-side vector establishing the interface condition and calculated by Eq. (7), where U w is the wake-induced flow speed.
[α]Γ = H (5) On every panel the vortex position X b is placed at a quarter of the panel length, whereas the control point X c , where the system is solved, is placed at three quarter of the panel length. This ensure implicitly the Kutta condition.
Wake vortices are initially placed at a distance 0.25.U ∞ .∆t downstream of the trailing edge's position. Then the model computes the local flow speed at each particle's location and moves the wake vortices accordingly.
Once the vortices circulation is calculated by multiplying the inverse of the influence coefficient matrix α and the vector H (Eq. (5)), the pressure difference across the structure can be calculated by Eq. (8) and the hydrodynamic forces by Eq. (9). dl is the element length, i.e. the distance between two vortices.
Mechanical power transmitted from the flow to the membrane is then the product of hydrodynamic forces and the structure's speed (Eq. (10)).

B. Co-rotational finite elements method
The structure model aims to describe the motion of the membrane under external dynamic loadings. The structure's displacement, velocity and acceleration generate stiffness, damping and inertial forces. The general form of the system's dynamic equation is Eq.
As noticed during previous experiments [1], the membrane undergoes small deformations and large displacements, which generate non-linear behaviours. This explains the use of the co-rotational finite elements method for the computation of the stiffness matrix. This formulation uses a fixed frame for the rigid body motion and expresses local deformations in a co-rotational frame attached to the element [4]. The membrane's thickness being small compared to its length, the Euler-Bernouilli beam hypothesis can be applied: the beam's cross-section is considered to remain normal to the element neutral axis after deformation [10].
The structure is divided into N parts. In its initial configuration, the beam element is defined by its nodes coordinates (x 1 , y 1 ) and (x 2 , y 2 ) in the global frame. There is 6 degrees of freedom by elements, namely two translations u and v and a rotation θ for each node. The global displacement vector X is defined by Eq. (12) for an element: The element is then defined in a new co-rotational local frame, with 3 degrees of freedom in which the local displacement vector is expressed by Eq. (13). Fig. 3 illustrates the parameters describing the element local and global displacements. It shows the relation between global and local variables. VectorX is obtained through Eq. (14). The element initial length is l 0 and the current length in the co-rotational frame l is calculated by Eq. (15). β 0 is the element initial angle and β its current angle in the co-rotational frame (Eq. (16)). θ 1 and θ 2 are the nodes global rotations,θ 1 andθ 2 are their counterparts in the local frame (Eq. (17)).
The relation between the local and the global displacement vectors for an element is obtained through Eq. (18). B is the transformation matrix calculated by Eq. (19). c and s respectively refer to cos(β) and sin(β), δ is the differential operator.
The transformation matrix B also links the local interior forcesq associated with local displacementsX and the global interior forces q associated with global displacements X : The relation between local forces and local displacement is obtained by Eq. (21), whereN is the local normal force,M 1 andM 2 are the nodes local moments, A is the area of the beam cross-section, E the Young modulus and I z the second moment of area.
The stiffness matrix The damping matrix is estimated as being proportional to the stiffness matrix (Eq. (24)). This kind of damping is a specific application of Rayleigh damping that is equivalent to Kelvin viscous damping model [11].
The lumped mass matrix is used, it is assembled from elementary mass matrices calculated by Eq. (25). C. Interaction process and system solver Fluid-structure coupling is here based on a partitioned approach. Fluid and structure models are solved iteratively using the results of each other. This approach is modular, which facilitates the calculations implementation and/or modification. However, the codes are solved successively and not simultaneously. This introduces an additional numerical error that can increase during the simulation. Several strategies are used to limit this error.
An iterative process on both models is implemented at each time step, with the respect of a convergence criteria (Fig. 4). This kind of coupling improves the interaction quality and enables to converge to the same solution as a monolithic approach [12]. It is used in high mass ratio applications such as offshore engineering. Further methods are used to ensure the convergence, such as the added mass implementation in the structure code [13] and the Aitken adaptive relaxation [14].
The added mass is the virtual mass of fluid matching the inertial part of the fluid forces, that is to say, the force component proportional to the body acceleration [15]. In the present study, the fluid density is of the same order as the solid one. The added mass is then important and affects significantly the structure's behaviour resolution [16].
A quick method for added mass estimation of slender body is used in this paper (Eq. (26)). The added mass term is proportional to fluid density ρ f and to the element surface S to the power of 3/2. It is defined by a coefficient k ma , that only depends on the structure's geometry [13]. Here, k ma = 0.38. The added mass term is then added to the diagonal terms corresponding to translations degrees of freedom in the solid mass matrix.
The system equations is then a second order differential equation. The selected method of numerical integration is the average acceleration method (Eq. (27) & (28)). It is widely used in structure dynamics because it is implicit and unconditionally stable [17].
The system equation can then be written as Eq. (29), is the system state matrix and R = F ext − F inertiel − F amort − F int is the numerical residual vector.
[A]δẌ = R The state matrix is then modified in order to take into account boundary conditions and imposed displacements. In this case, pre-stress cables are modelled as a bar element and the membrane is clamped at its leading edge.
The solution is obtained by the Newton-Raphson algorithm, by iterating onẌ (t+∆t) in order to minimize the residual. The iteration starts from an estimated solutionẌ g (t + ∆t) =Ẍ (t). This estimation is used to calculate the state matrix [A], the force vector F and the residual R.
Aitken method is also used to accelerate the convergence [14]. A relaxation factor h i calculated according to Eq. (30) is multiplied to the solution increment [A] −1 R. i being here the iteration number.
• If the norm of the residual is inferior to the convergence criteria , then the convergence is established:Ẍ (t + ∆t) ←Ẍ g (t + ∆t) (31) • Else, the solution is actualized and a new iteration is done: Initialization process includes a ramp on cables withdraw length and a temporary vertical loading to initiate the structure motion. This ramp facilitates the beginning of the modeling procedure. It is calculated by Eq. (33) for the cables withdraw length and by Eq. 34 for the temporary vertical loading, with f rampe (t < t 1 ) = f rampe (t 1 ) and f rampe (t > t 2 ) = f rampe (t 2 ). Here,

D. Experimental set-up
The 1/20 th scale prototype of an undulating membrane tidal energy converter is described in [1] and summarised here. Its results are used for comparison with the numerical model. Geometric and dynamic scaling have been used with respect to structural similitude. The prototype is made of a polyacetal plate of dimension: L × L × h = 0.8 × 0.8 × 0.003 m 3 . It is stiffened by 6 transverse carbon/epoxy bars to ensure a 2D motion. Upstream and downstream rigid flaps lengthen the membrane by L f lap = 0.15L. The membrane is attached by three polyacetal bars of stiffness k f ix = 6250 N/m, themselves clamped to the frame that maintain the prototype in the middle of the 2meters-deep water column.
The structure is pre-stressed by compression cables that buckle it in the axial direction (Fig. 5). Compression cables connect the downstream end to the frame. They are shorter than the distance at rest between the attachment points and are characterised by the withdraw length ratio d (Eq. (35)). The cables put the membrane in buckling before it is submitted to axial flow.
The prototype has been tested in the IFREMER flume tank of Boulogne-sur-Mer [18]. Instrumentation of tank tests is composed of a 6-components load cell attached to the prototype's frame and a motion tracking system using 7 LEDs as targets that are located on the membrane's side. Sampling frequency is f s = 100 Hz and measurement time is ∆t = 80s. A PIV wake characterisation has also been carried out and presented in [19]. The membrane motion obtained with these parameters is plotted in Fig. 6 and 7 and compared with experimental results. The trajectory, as well as the amplitude and the frequency, are really close to those observed during experiments. Indeed, the simulation gives an undulation amplitude of a = 0.38L and a frequency of f = 0.25u ∞ /L. The differences with the experimental model are respectively +4.2 % and -4.3 %. The shift in the membrane average vertical position is due to the absence of gravity in the numerical model.  The evolution of lift and drag forces is represented in Fig. 8 and 9. The numerical model manages to reproduce the hydrodynamic forces on the membrane. However, there is a high variability of numerical results, with what looks like a parasite signal on a frequency at about ≈ 3 Hz and amplitude ≈ 100 N.
This problem could be due to the mass and damping matrices that are oversimplified. Indeed, Rayleigh damping is supposed to be a linear combination of stiffness and mass matrices (see, for example [20]). The linear coefficients of this combination are calculated from the damping of two deformation modes. It is then possible to specify one damping rate for the first unstable mode and another for the highest mode acting in the calculus.
Furthermore, the added mass matrix is here considered to be a diagonal matrix where the terms corresponding to rotational degrees of freedom are cancelled. This formulation is certainly too simplistic and should be improved (see for example [13]). The value of the added mass coefficient k ma could also be adjusted to reflect more closely the inertial component of fluid forces. If it is increased, it will limit velocity variations at high frequency and consequently restrict the hydrodynamic forces variations visible in Fig. 8.

F. Wake results
Axial and vertical speeds are presented in Fig. 10 for six instants of the undulation cycle from top to bottom. The membrane is plotted as a white line. Areas of acceleration on the concave sides of the membrane and deceleration on the convex sides are visible in the axial speed map, on the left of Fig 10. These areas alternate and propagate along the membrane with variable size according to the membrane motion. As observable on the right part of the figure, the vertical motion of the membrane produces areas of vertical flow speed of the same direction. These areas also propagate in the current direction. Fig. 11 presents the flow field at a given position of the membrane. On this picture the deceleration and acceleration areas are clearly visible. We can also see how the vertical speed is impacted by the motion of the membrane. These flow perturbations then separate from the membrane at its trailing edge and shape the wake structure as they move downstream.
The evolution of the circulation emitted in the wake is presented in Fig. 12. As expected, emitted circulation oscillates around zero, it is positive during one half of the undulation cycle and negative during the other half, which explains the coiling motion of the particles in the wake and the composition of the wake maps ( Fig. 13 and 14). Here also, there is a high variability of the signal. This error can affect the displacement of wake vortices and induce an error in the following wake maps. Therefore, the following results must be considered with care until this variability problem is solved   Wake vortices positions are plotted as blue dots in Fig. 13 on a instantaneous map of the flow speed magnitude. The motion of these vorticity particles according to computed local flow speed enables the model to calculate the wake flow field more precisely and especially its vertical expansion. The wake presents a ribbon of flow perturbations with a sinusoïdal shape. It is composed of an alternation of faster and slower flow speed areas. Point-vortices positions in the wake also show an alternation of positive and negative vertical speed areas. These areas are generated by the vertical motion of the membrane and correspond to half a cycle each. Indeed, undulation period being T = 4 s and flow speed U ∞ = 0.8 m/s, each area is about 1.6 m long, that is ≈ T U ∞ /2. Mean horizontal flow speed map is presented in Fig 14. It shows that the vertical expansion of the wake seems to have a parabolic shape according to the horizontal position. There is an acceleration of the flow around y = 0 and on the top and the bottom of the wake. Between these areas, the flow is decelerated. This is consistent with the experimental results presented in [19] for a similar, but different configuration.
These results are satisfying but there are still some differences with experimental results. This can be due to the absence of flaps and fixation bars in the numerical model. Diffusion and dissipation of the wake particles vorticity could also improve the results. The use of more elaborated damping and added mass matrix should correct the circulation of wake particles, hence their displacement.

IV. CONCLUSION AND PERSPECTIVES
The 2D fluid-structure interaction numerical model presented in this paper is based on an implicit partitioned coupling between an unsteady point-vortex method and a co-rotational finite element formulation. Its results are close to the experimental model's ones for an undulating membrane tidal energy converter without power take off. The high quality of the results despite various hypotheses enables to validate the numerical model for this specific application.
Furthermore, the differences are easily explainable and will be corrected in the near future. For example, the addition of gravity into the numerical simulations could adjust the membrane's trajectory mean vertical position. The use of more elaborate damping and added mass matrix should also improve the results' quality. A numerical model validation on several configurations and a wake comparison with experimental results are also planned.
Next step for the numerical development will be the addition of a power take off model. Then, flaps and fixation arms will be modelled. Different fixations for the membrane could possibly be tested through the addition of pivot and/or spring liaisons.
In a second phase, it will be possible to add an upstream turbulence model [21], to add structure thickness in order to test different geometries and finally to develop a 3D model. This 3D model would enable to study the side effects, to represent better the force repartition and the wake and to adapt the tidal converter to its environment. Once this interaction model validated, it can be used for system optimization.