FSI analysis and simulation of flexible blades in a Wells turbine for wave energy conversion

In this paper a preliminary design and a 2D computational fluidstructure interaction (FSI) simulation of a flexible blade for a Wells turbine is presented, by means of stabilized finite elements and a strongly coupled approaches for the multi-physics analysis. The main objective is to observe the behaviour of the flexible blades, and to evaluate the eventual occurrence of aeroelastic effects and unstable feedbacks in the coupled dynamics. A series of configurations for the same blade geometry, each one characterized by a different material and mechanical properties distribution will be compared. Results will be given in terms of total pressure difference, supported by a flow survey. The analysis is performed using an in-house build software, featured of parallel scalability and structured to easy implement coupled multiphysical systems. The adopted models for the FSI simulation are the Residual Based Variational MultiScale method for the Navier-Stokes equations, the Total Lagrangian formulation for the nonlinear elasticity problem, and the Solid Extension Mesh Moving technique for the moving mesh algorithm.


Introduction
Renewable energy is consistently gaining more and more relevance in terms of shares of global electricity generation, growing worldwide from 14% in 2005 to 22% in 2013 [1], and its expected almost to double in 2040 [2]. In spite of the great potential of ocean waves energy as renewable energy source, the production of electricity from ocean waves is still limited if compared to more traditional clean and renewable energy sources such as wind, hydroelectric, solar and bioenergy, and the technology of the conversion system for waves energy recently reached a stagnation point. The reasons of the deceleration in the development of new technologies are multiple, but they all hinge around the predicted costs of wave power, evaluated as non-competitive with respect to other sources [3]. In analogy with others renewable -and non-energy sources and energy conversion systems, a reduction of the predicted costs of waves energy can be achieved by increasing the reliability of the equipment, increasing the efficiency of the conversion system, and reducing the capital investment costs [3]. In this framework, Wells turbines already represent a largely developed and consolidated technology for wave energy conversion in Oscillating Water Columns (OWC) devices. Traditional Wells turbines present limits in terms of average efficiency, because they mostly operate in off-design flow conditions. Furthermore, their peak efficiency is low, because of the symmetric, untwisted blades. It is trivial to observe that a nonsymmetric, twisted blade profile could improve the performance of the machine. A Wells turbine with flexible, polymeric blades could solve at once multiple issues: a flexible blade can be designed in order to improve efficiency in off-design conditions (i.e. low mass flow rates) and to improve peak efficiency due to a passive-adaptive blade morphing, in both flow directions, without any need of guiding vanes. Such blade could also improve the reliability of the machine in a highly corrosive environment such as saline aerosol, and it can drastically reduce the capital investment and operative costs due to the adoption of a cheaper material with respect to steel. This preliminary study presents a numerical analysis of a Wells turbine with flexible blades by performing a coupled 2D computational Fluid Structure Interaction (FSI) dynamic simulation using stabilized Finite Element Method (FEM), for the solution of the fluid field and the structural dynamics. The models adopted in the simulation are the Residual Based Variational MultiScale method (RBVMS) to model the flow aerodynamics [4] and the Total Lagrangian formulation (TL) to model the non-linear structure mechanics problem [5]. The consecutive movement of the mesh will be computed using the Solid Extension Mesh Moving Technique (SEMMT) [6]. To compensate the absence of a wall function and a turbulence closure at the wall, the Dirichlet boundary condition at the wall are imposed weakly [7], as suggested by the same authors of the RBVMS model. The three subsystems are strongly coupled in time. The strong coupling technique is not often used for dynamic problems in turbomachinery, due to the fact that in this kind of applications, the stiffness of the blade is very high and the aerodynamic feedback is secondary with respect to the inertia forces, and often a weak coupling approach [8] or a modal analysis [9] are preferred. However, in this application the high flexibility of the blade required a strongly coupling approach. Simulations are performed with an in-house build software, FEMpar, developed for previous CFD analysis in turbomachinery [10], and which represents a development of the software XENIOS++ [11], previously used by the authors for FSI [12,13] and particle tracking simulations in turbomachinery [14].
The analysis can provide useful insights about the deformations and stresses of the blade and about the performance improvements due to the flexibility of the blade in different operative conditions. Complex FSI dynamics and large deformations are expected, indeed, the low stiffness structure is selected such to be able to deform under the aerodynamic loads. In addition to that, the deformation and the variation of the shape must be evaluated with respect to the OWC period and to the dynamic response of the blade, and consequently, the choice of the material and the blade thickness can be accordingly designed.
In this paper, a series of alternative solutions for the material of the blade of the Wells turbine, in terms of elasticity and density distribution, are proposed. As a preliminary study for this application, no optimization for the choice of the material is performed, but rather there's a specific focus on the observation of the differences in the flow field and structural response of the blade, which is very thin, flexible, and constrained at the leading edge and trailing edge. An iterative, trial and error approach to define the material for each configuration is adopted, starting from a standard plastic material used for marine applications.
2 Turbine description

Blade geometry
In this work, a simulation of a peculiar design of a small Wells turbine will be presented. The design is directly descending from the DIMAFP-TW1.5, a prototype 1.5 Kw 3000 rpm Wells turbine for Mediterranean operations (i.e. designed to operate at an average wave height of 0.8m and an average period of the waves of 4s [15]), with a 9 bladed 0.5m diameter rotor, and two stators with 17 blades each. [16,17]. Here, the object of the simulation is a 7 bladed, rotor only Wells turbine, with a NACA 0002 profile: such a low thickness profile is required to provide the blade with the desired flexional behaviour with realistic materials. The chord distribution is constant along the blade span. Further details on the blade can be found in Tab. 1, while a 3D sketch of the Wells Turbine can be seen in Fig. 1. In this preliminary design of this flexible blade turbine, the blade is constrained to the hub as if two rigid clamps are holding the leading edge and the trailing edge along the entire span of the blade, as illustrated in Fig. 1.  Rotational speed (rpm) 3000

Material
For this paper, the objective is to simulate a flexible blade for a Wells turbine, without recurring to unrealistic materials, or unrealistic elasticity distributions. The reference material is Nylon 66/6 [18], a widely used polymer, reinforced with 30% glass fiber. As will be presented in section 4, several variations of this material will be explored, with the objective to observe the dynamic response, deformations, and stable configurations for the blade. Additional details about the material will be provided for each case. In particular, there will be a comparison between the flow around an homogenous isotropic Nylon 66/6 blade against other configurations with a distribution of elasticity and density differentiated in three zones: a zone next to the leading edge, a zone next to the trailing edge, and a central larger zone. This strategy allowed more control over the displacements of the blade, with the opportunity to compensate, for example, the extreme thinness at the trailing edge. A scheme of the subdivision in zones of the blade can be found in Fig. 2. The names of the zones, i.e. leading edge (LE), cord centre (MID), and trailing edge (TE), will be used in the upcoming sections to identify each zone.

Mathematical modelling
In this section is briefly presented an overview of the models adopted for the aerodynamic, structure mechanics and moving mesh problems.
To model the aerodynamic subsystem with the finite element method, the residual based variational multiscale formulation (RBVMS) proposed by Bazilevs et al. [4] is selected among other methods. The RBVMS is a particular declination of the family of the variational multiscale methods [19]. Such methods find their motivation in the fact that definition of large scales by spatial averaging leads to problems (in particular at the boundary). The idea behind them is to define scales by projection into model function spaces: the spaces and subsequently their functions are decomposed in coarse and fine scales by a multiscale direct sum decomposition. The large scales are defined by the finite element discretization and they are effectively resolved, while the small scales are modelled as function of the residuals of the large scales' equation.
RBVMS has been proven as a solid and attractive alternative to LES approach, as presented for example in [20], where the model is validated against DNS data and compared with the classical dynamic Smagorinsky LES [20]. Further details on the RBVMS models can be found in [4].
When performing structural analysis on highly flexible material in slender shapes, the linear elastic model is not sufficient to effectively describe the entire deformation range. In such cases, the dynamic is strongly non-linear, even when remaining in the elastic field. Some non-linearities, i.e. the geometric non linearities, requires formulations which accounts for the nonlinear elasticity, such as the Total Lagrangian formulation [5]. With this formulation it's possible to add the nonlinear elasticity contribution to the usual linear stiffness matrix, by taking into considerations the Green-Lagrange deformation tensor, and the Piola-Kirchoff stress tensor, evaluated with respect to the original configuration of the structure, at the beginning of the dynamics.
In a FSI simulation, the motion of the mesh is required to capture the movement of the interface between solid and fluid. This is a crucial aspect of the solution algorithm both in terms of computational effort, and in terms of mesh quality. In fact, the solution of remeshing is often not feasible, and the mesh at the interface must not be excessively distorted when dealing with large deformations (it is required, for example, to keep the wall distance of the first layer of elements up to a specified value). For those reasons, the selected algorithm for the moving-mesh is the SEMMT, proposed by Tezduyar et al. in [6]. In such a model, the element of the mesh is considered as a deformable solid, with a stiffness related to the Jacobian determinant of the element, and every point of the mesh is moved by solving a simple linear problem.

Software description
FEMpar is a in-house software recently developed, built as a second iteration of the software XENIOS++, for 2D-3D CFD-FEM simulations of aerodynamic and fluid structure interaction problems. FEMpar relies on its own parallel finite element and methods library and has been developed to improve scalability (speedup and memory management), and customizability (new physical models or FEM data structures addition). FEMpar adopts the same solution algorithm for the solution of the FSI problem of XENIOS++, a strong coupling block iterative solution for the multiphysics problem (aerodynamic and structural). Further implementation details can be found in [10]. For the solution of the linearized systems, FEMpar uses the GMRES solver for sparse matrices from the linear solver algebra PETSc [21], and it relies on the Generalized-α time integration scheme proposed by Chung et al. [22].

Numerical setup
A 2D triangular mesh is set up with boundary layer inflation at the wall. The computational domain is a rectangle obtained extracting the cascade section of the turbine at radius r=246mm, corresponding to the 95% of the blade span. The very small ratio between the axial and tangential velocity components at the inflow and the high solidity of the turbine, allowed to use a relatively small domain. The rectangle extends, in the axial direction, for one chord upstream and one chord downstream. The resulting mesh is shown in Fig. 3, it contains about 5000 fluid nodes and 1500 solid nodes.
An inlet relative velocity condition of v=78.3 m/s is imposed at the inlet boundary. The relative velocity is obtained by considering the nominal flow rate of Q = 0.85m 3 /s and the rotational speed of 3000 rpm. The simulated wells turbine is operating at max load conditions, to obtain the maximum deformation of the blade. At the outflow boundary, a zero gradient condition is imposed. The remaining boundaries are considered as periodic boundaries, so that the domain can be considered as a periodic cascade. At the wall, the no-slip condition is imposed in a weak form, according to Bazilevs [7]. This technique will result in an effective relaxation allowing a slip velocity when a not sufficient minimum wall distance is locally encountered to resolve the flow at the boundary layer. The first cell is at a distance of 3.3 × 10 m and Re ~ 3 × 10 . The simulation is carried on for 4000 time steps, with a dt = 1 × 10 s, corresponding to two revolutions of the turbine. The FSI simulation starts from the full resolved dynamic flow field solution, obtained by keeping the blade rigid. The boundary conditions are highlighted in different colours (yellow for inlet boundary, black for wall boundary/interface, green for periodic boundaries).

Test cases
Starting from the NACA 0002 symmetric blade, the analysis begins from a CFD study where the performance of a series of blade profiles deformed a priori are evaluated, where the curvature degree varies by modifying the camber line, considering a curvature coefficient m as a percentage of the chord length, with m respectively equal to 2, 5, 7, 9, 11 and 15%, by assigning such value to the m coefficient in the expression of the camber line (1): Where yc and x are the coordinate of the camber line, p is the position where the blades has maximum curvature and m is the maximum curvature express as percentage of the chord. From the simulations, a series of characteristic curves are obtained as shown in Fig.4, in terms of total pressure drop. When looking at those lines, it is almost a natural consequence to consider the possibility to obtain an envelope curve, as shown in Fig. 5, only by the passive morphing of the flexible blade, as if the material were able to deform in such a way to assume those positions as operative configurations.  With such background, the main objective is to evaluate the possibility of adopting an elastic material, characterized by low stiffness, to obtain a comparable effect to the one observed in Fig. 4, by exploiting the passive morphing capabilities of such a material. By doing so, a complex series of new vibrational dynamics -that must be considered when observing the total dynamic-are introduced, and for such reason both frequency analysis tools and FSI simulations are performed.
To reach this objective, it is useful to compare five different flexible blades characterized by a different material. All the blades where simulated under the same inflow conditions given in Sec. 4.2, so that the test cases are illustrated in Tab. 2 in terms of Young's Modulus E, density ρ and Poisson's ratio ν. Case A refers to the Nylon 66/6 blade, and case B, C, and D refer to the blades with modified elasticity and density distributions, according to the scheme of Fig. 2.

Results and discussion
The main objective for this paper was to evaluate the possibility to realize flexible blades for a Wells turbine. The dynamic behaviour of the five blades defined above, plus a fixed blade considered as a reference, in a cascade configuration is obtained by the observation of the CFD results. Each blade deformed in a different way, some of them barely moved around their original position, others presented an unstable dynamic behaviour because of a positive feedback between solid mechanics and fluid dynamics of the field.
In order to evaluate the difference in performance for the different blades, it is possible to see in Tab. 3 the total-to-total average pressure drop, computed in two different sections, upstream and downstream with respect to the blade, both at the same distance of half of the chord length. From the table stands out that the larger pressure drop is obtained with the configuration D. In addition to that, it is possible to observe that all flexible blades present a larger pressure total difference with respect to the fixed blade, encouraging further optimization of the design and of the choice of the material. A primary remark of this work is that the choice of the right material is crucial, because the final shape of the blade, the stability of its dynamic response, and consequently the performance of the turbine, are all affected by the choice of the material more than any other design variable. In order to highlight the occurrence of aeroelastic coupling, it's useful to see in Fig. 6 the comparison of the displacement of the center of the blade (e.g. in configuration D) when it is subjected to a uniformly distributed 100N impulse in the ydirection, and during the FSI simulation. In Fig. 7 the main frequencies of the dynamics are compared by plotting the normalized Fast Fourier Transformations (FFT) of the computed displacement and the static pressure evolution. It's easy to notice from Fig. 6 both the variation in frequency of the main dynamics and the occurrence of new modes (due to the aeroelastic coupling). In addition to that, in the same Fig.6 stands out that the main frequency of the pressure is perfectly overlapping one of the secondary modes of the emerged from the FSI simulation, and therefore it's possible to conclude that such secondary mode, not present in the structure-only dynamic, it's entirely originated by the aeroelastic coupling.  A flow survey obtained with the 2D FSI simulation is reported. For the sake of brevity, only images from the configuration D, which resulted the most promising configuration according to Tab. 4., are reported for the flow survey. Fig. 8 presents the total pressure and the velocity magnitude, both averaged in time. Fig. 9 shows two instantaneous fields of static pressure and velocity magnitude, both describing the same time step, at which the displacement of the central probe point was maximum.  For the sake of completeness, in Fig. 10 it is possible to see the Von Mises stress distribution for the test Case A (i.e. characterized with homogeneous Nylon 66/6). In the picture, only the stress distribution for Case A is shown, because it is only possible to compare the simulation results with real yield stress values for the Nylon 66/6.
The yield stress value considered for the material is Y=162 MPa , and as shown in the picture the entire blade is below the threshold value set as 0.2Y for safety reasons, with the exception of the zones near the constraints at the leading edge and trailing edge. Those peaks in the stress distribution suggest to modify the fixed constraints with a different constraint, i.e. a hinge. In any case, the distribution still proves the correctness of the initial choice of the material.   It is easy to observe that the average displacement varies by considerable amount from one material to the other, and that the obtained curvature of the blade by the adoption of those materials configurations, are not large enough to be comparable to the a priori imposed curvatures of the blades tested to obtain the curves of Fig. 4, but the effect on the flow of the passive morphing blade is still considerable.

Conclusions
In this work is presented a design and 2D CFD FSI analysis for a small Wells turbine with flexible blades using an in-house built software. The study begins from the consideration that an a priori imposed curvature on the blade produces different total pressure drop across the blade with respect to the straight blade. In order to observe the dynamic behaviour of a lowstiffness blade when the flow is at nominal condition, five different blades configurations are set up, each one with its peculiar material. The materials chosen for the test cases are different in terms of density and Young's Modulus, and the reinforced polymer Nylon 66/6 is considered as the reference material. With the aid of FSI-CFD simulations, the insurgence of new dynamics effect is highlighted, in terms of vibrations modes rising from the coupling of the fluid dynamics and structural dynamics, thanks to the strong coupling technique of the physical sub-systems implemented in the software. These modes of vibration are entirely excited by the aeroelastic effects. A comparison of the CFD results for each blade in terms of total pressure drop is given, and the most promising configuration is selected, as the obtained value of the total pressure drop is almost three times the value obtained with a fixed, straight blade. The results are encouraging further studies, supported for example by 3D CFD and possibly a new design iteration with new constraints and new materials, supported by an optimization process, to select an optimal material and possibly to reach stable configurations with high curvature of the blade.