1 Introduction

Fluid–structure interaction (FSI) is a topic spanning through many engineering disciplines from aeronautics and renewable energy to biomedical engineering and particle transport. Traditionally FSI has been treated as an aeroelasticity or hydroelasticity problem, where reduced-order modelling has been used along with an assumption of small deformation for the solid body [1]. Unsteady flow dynamics that is usually based on potential flow theory and sometimes with a boundary layer correction is then coupled with structural dynamics that is often modelled using a system of equivalent springs and dampers as for wind turbine blades [2]. A similar approach can also be used for acoustic–structure interaction, where a fluid-phase wave equation is coupled with a solid-phase wave equation by requiring continuity in the displacement and pressure at the solid–fluid interface [3].

Although the reduced-order modelling has been found successful in various problems ranging from wings flutter to feathering of wind turbine blades, it has limitations in terms of the simplified assumptions used in the fluid and solid dynamics. Extensions to include large deformations and computational fluid dynamic (CFD) results through reduced-order modelling have been suggested [1], but they also point to the need for a full computational FSI solver. Such solvers exist as monolithic or partitioned [4, 5]. The monolithic solver uses the same solver for the fluid and solid phases. An example is the use of the structural dynamics Y-Code based on the combined finite element method–discrete element method (FEM-DEM) to simulate the propagation of an acoustic pulse through the solid part of a hammer into liquid metal [6]. This was done by modifying the elastic properties in the solver for the liquid phase in order to agree with the wave equation for the liquid. Another example is the particle finite element approach, where an assumption of a quasi-incompressible flow is used for the fluid [7].

The monolithic solver can also yield a stable information transfer between the fluid and solid phases. However, it has the disadvantage of not using the well-established optimised algorithms that were already separately developed for the fluid and solid phases [5]. This is the advantage of the partitioned solver, using specialised solvers for the fluid and solid phases [5], and it is the solver to be pursued in this study.

Partitioned solvers can use deforming meshes that adjust to the displacement of the fluid–solid interface or stationary meshes. The arbitrary Lagrangian–Eulerian approach (ALE) is a common approach for the deforming mesh approach and is widely used in commercial packages. Its advantages are that it can provide a clear cut between the fluid and the solid and an optimal level of grid clustering near the fluid–solid interface. However, ALE can become complicated for bodies of large deformations, problems of contact between bodies and problems containing many complex geometry bodies.

The immersed boundary method (IBM) is commonly used for the stationary mesh approach [4, 5]. It provides a simple geometry grid that can deal with complex geometry bodies, large deformations and a large number of bodies. However, it can smear the interface between the fluid and the solid, and when implemented along with a partitioned solver, the problem of information transfer between the fluid and solid may lead to instability if implemented in a weak formation [4]. The latter means that no convergence procedure has been performed as part of the IBM in order to make sure that the fluid solver and the structural solver show the same displacement and stress at the interface at the same time [4]. This study will use the IBM approach with an iterative procedure to ensure the stability of the information transfer at the interface between the fluid and solid phases.

Two FSI topics will be dealt: aquatic vegetation in channel flow and hydro-/aerokinetic turbines of renewable energy. Both are similar by requiring simulation of flow around flexible cylinder-like slender bodies, which are the vegetation stems or the turbine’s blades. Both topics are of strong current engineering and environmental interest. They also present computational challenges of capturing the flow’s relevant spatial and temporal scales while accounting for the structural response of the flexible bodies. However, the two topics also exhibit differences in the physical behaviour of the flow and structure, most notably due to the difference in the level of the flow Reynolds number and the rigidity of the body. The typical Reynolds number of the vegetation stem is less than a few thousands when based on the stem’s diameter [8, 9], while the blade’s Reynolds number is at least about 50 k for small turbines when based on the blade’s chord length and going over 1 M for larger turbines [2, 10]. Vegetation stem can be highly flexible, while a tidal marine turbine blade is much more rigid.

Traditionally, aquatic vegetation has been seen as nuisance in waterways, slowing down the flow and increasing the risk of flooding. Thus, the interest was mainly in calculating the drag caused by vegetation in order to support a policy of dredging the waterways. However, studies in the last two decades have also pointed out the positive environmental effects of aquatic vegetation that include reducing the erosion of river bed by reducing sediment, acting as filters against pollutants and supporting an aquatic ecosystem [8, 9]. Aquatic vegetation can be classified as floating, emerging, i.e. piercing the water surface and submerged. In this paper, we will focus on the submerged vegetation.

Much research has been directed into identifying the effects of vegetation density, categorising it as sparse or dense, its blockage area as relative to the waterway’s cross section and its proximity to the water surface. CFD and experiments where the vegetation stems were modelled as cylinders have shown sparse vegetation to have higher turbulence intensity and sediment activity than in dense vegetation [8, 9]. Inflection points in the mean axial velocity profile are also dependent on the vegetation density and height, affecting the flow instability and creation of vortices. Flexibility of the vegetation has been modelled using an immersed body approach, where in order to reduce computational cost the stems were modelled using one-dimensional (1D) structural dynamics [11, 12]. Good agreement was generally achieved for the time-averaged flow between the FSI computations and similar experiments, but agreement in the bulk drag prediction could be improved. One-dimensional structural dynamic modelling for the stems will also be used in this study due to its relatively low computational cost while using the FEM-DEM methodology.

Hydro- and aerodynamic modelling of kinetic turbines traditionally relies on blade element methods [2]. These low-order models can be coupled with beam models for the blades in order to deal with FSI. Recent CFD development along with the ALE approach has delivered several studies where FSIs of horizontal and vertical axis kinetic turbines have been studied. The blade was modelled using shell elements to reduce computational cost and it was found that the blade’s flexibility could reduce the power production [13, 14]. IBM has been less common for studying kinetic turbines because of its reduced ability to capture sharp interfaces as compared to ALE. Nevertheless, it has been successfully applied to marine tidal turbines, yielding excellent agreement with experimental results for the power production [10], while utilising the simplicity of the grid to achieve high computational efficiency in parallel computing infrastructure [15]. Recently, interest has risen in using the blade flexibly in order to actually achieve higher hydro-/aerodynamic performance through the concept of smart structure and thus getting optimal geometry of the blade for different flow scenarios [16]. FSI simulations can have a significant role in achieving this aim.

The computational methodology of the FSI simulations used in this study is presented in the next section. It is followed by results analysis of flexible submerged cylinders, modelling aquatic stems and FSI analysis of horizontal axis kinetic turbine blades of the common type based on an aero-hydrofoil section as well as of spinning (Magnus) blades.

2 Methodology

FSI methodology was developed and implemented for studying deformable bodies embedded in incompressible viscous flow and/or containing such flow. A partitioned FSI solver named CgLES-Y has combined three advanced numerical technologies and the corresponding computing codes. The fluid solver is based on the finite volume approach, while the structural dynamic solver is based on the FEM-DEM scheme [17]. An IBM algorithm ties the two solvers together. This FSI methodology has already been successfully implemented for a range of problems from bio-fluids of red blood cells flow [18, 19] and urine flow in the ureter [20] to sediment [21, 22] and marine tidal turbines [10]. In this section, the main points of the method are highlighted, where further details are given in Refs. [10, 15, 21, 23].

2.1 The flow field solver

The three-dimensional (3D) incompressible Navier–Stokes (NS) equations are simulated using a rectangular grid:

$$ \frac{{\partial u_{i} }}{{\partial x_{i} }} = 0 , $$
(1)
$$ \frac{{\partial u_{i} }}{\partial t} + \frac{{\partial \left( {u_{i} u_{j} } \right)}}{{\partial x_{j} }} = - \frac{1}{\rho }\frac{\partial p}{{\partial x_{i} }} + \upsilon \frac{{\partial^{2} u_{i} }}{{\partial x_{j}^{2} }} + \frac{1}{\rho }\frac{{\partial {\tau}_{ij} }}{{\partial x_{j} }} + f_{i} , $$
(2)

where \( t,u_{i} ,p,\rho ,\upsilon ,{\tau}_{ij} \) are time, a velocity component in the xi direction, pressure, density, kinematic viscosity coefficient and the large eddy simulation (LES) sub-grid scale (SGS) stress, respectively. The SGS stress was modelled using the mixed timescale (MTS) model that has a mechanism to turn itself off at laminar flow [24]. In case of direct numerical simulation of low Reynolds number, the SGS stress was taken as zero. fi is the body force in the xi direction, and this is the term where the IBM force is applied on the flow field.

The NS equations were discretised using the finite volume approach where the spatial derivatives were calculated using second-order central schemes. Grid clustering in areas requiring fine computational grid was achieved using grid mapping. Continuity was achieved using the projection method that was coupled with the time marching methods of third-order Runge–Kutta method or second-order Adams–Bashforth [15]. The Poisson equation for the pressure that results from the projection method was solved using rapid solvers based on the bi-stab approach [15].

2.2 The structural dynamic solver

The FEM-DEM approach was used to simulate the dynamic equation [17]:

$$ \frac{{d^{2} \vec{y}}}{{dt^{2} }} = \vec{F} + \frac{1}{\rho_{s} }\nabla \cdot \overline{{\bar{\sigma }}} , $$
(3)

where \( \vec{y},\vec{F},\rho_{s} ,\overline{{\bar{\sigma }}} \) are the element’s displacement vector, body force vector, the solid body density and the Cauchy stress tensor, respectively. Elastic behaviour was used to relate the stress to the strain for the problems discussed in Sect. 3. For the vegetation stems problem, Eq. (3) was modelled using 1D elements, reducing it to a flexural wave equation of a beam as in Refs. [11, 12].

2.3 FSI solver

Strong coupling through an iterative procedure was used to match the fluid and solid velocity at the interface while enforcing continuity. The stress was matched between the fluid and the solid every time step. The iterative procedure is based on the direct injection FSI approach, and only the main steps are illustrated here for the Adams–Bashforth time marching scheme. The reader is referred to Refs. [10, 15, 18,19,20,21,22,23] for further details.

Equation (2) can be marched from time n to time n + 1 using an intermediate velocity u* as follows:

$$ u_{i}^{ *} = u_{i}^{n} + \Delta t\left( {\frac{3}{2}h_{i}^{n} - \frac{1}{2}h_{i}^{n - 1} + \frac{1}{2\rho } \frac{{\partial p^{n - 1} }}{{\partial x_{i} }}} \right) , $$
(4)

where

$$ h_{i}^{n} \equiv - \frac{{\partial \left( {u_{i}^{n} u_{j}^{n} } \right)}}{{\partial x_{j} }} + \upsilon \frac{{\partial^{2} u_{i}^{n} }}{{\partial x_{j} \partial x_{j} }} + \frac{1}{\rho }\frac{{\partial {\tau}_{ij}^{n} }}{{\partial x_{j} }} $$
(5)

and h n−1i is similarly defined. The IBM f n+1/2i is found through an iterative procedure that yields from iteration k−1 to iteration k:

$$ f_{i}^{n + 1/2,k} \Delta t = D\left[ {v_{i}^{n + 1} - I\left( {u_{i}^{ *} - \frac{3}{2}\Delta t\frac{1}{\rho }\frac{{\partial p^{n,k - 1} }}{{\partial x_{i} }}} \right)} \right] , $$
(6)

where v is the desired velocity from the structural dynamic solver. D and I are the IBM distribution and interpolation around the IBM points illustrated in Fig. 1 and are specified in “Appendix”. The next intermediate velocity is found as \( \bar{u}_{i}^{k} = u_{i}^{ *} + f_{i}^{n + 1/2,k} \Delta t \), and one solves for pn,k in order to ensure continuity:

$$ \nabla^{2} p^{n,k} = \frac{2\rho}{3\Delta t}\frac{{\partial u_{i}^{k} }}{{\partial x_{i} }} . $$
(7)
Fig. 1
figure 1

Illustration of the immersed boundary method (IBM) layout

Full size image

The next iteration goes back to Eq. (6). The procedure stops when the change in \( I\left[ {u_{i}^{ *} - 1.5\left( {\Delta t/\rho } \right)\partial p^{n,k} /\partial x_{i} } \right] \) is small enough, leading to

$$ u_{i}^{n + 1} = \bar{u}_{i}^{k} - \frac{3\Delta t}{2\rho }\frac{{\partial p^{n,k} }}{{\partial x_{i} }} . $$
(8)

3 Results and analysis

As already noted, two types of FSI problems are investigated: submerged vegetation in a channel flow and kinetic turbine rotors. The FSI code has already been extensively verified in previous publications [10, 18,19,20,21,22,23], and thus, only verifications specific to the class of the investigated cases are shown followed by results of interest.

3.1 Vegetation in a channel flow

Submerged vegetation stems are modelled as circular cylinders following Refs. [11, 12]. The channel is taken as straight where the case of a curved channel is left for a future study. The channel’s bed is taken as no-slip wall, while the top of the channel is taken as free-slip wall, and hence, free surface effects are neglected. This can be justified by assuming a low Froude number \( {\text{Fr}} = U_{\text{S}} /\sqrt {gH} \), i.e. less than 0.5, where US is the typical velocity at the top of the channel and H is the depth of the channel [25]. Periodic boundary conditions are taken in the streamwise and spanwise directions. In order to sustain the flow in the channel, a body force fi is assumed and it corresponds to a gravity force acting on a channel with a shallow slope angle α ≪ 1 rad. This yields the following velocity profile for the laminar channel flow without vegetation:

$$ \frac{u\left( y \right)}{{U_{\text{S}} }} = \frac{y}{H}\left( {2 - \frac{y}{H}} \right) , U_{\text{S}} = \frac{{\alpha gH^{2} }}{2\upsilon } . $$
(9)

Hence, the bed friction velocity is \( u_{\tau } = \sqrt {\alpha gH} \). uτ is used to normalise the velocity field. The slope α is taken as 1/1000 rad in the following results.

The computational domain size is (6, 1.5, 4) d in the streamwise, stream-normal and spanwise directions, respectively, where d is the characteristic length scale used to normalise all lengths in the following results. The stem’s height is taken as 0.5 d and its diameter is 0.02 d. The flexible stem is assumed to be of elastic material with a Poisson ratio of 0.3. The density ratio between the stem and the surrounding water is taken as 1.1. The flow computational domain was discretised using 963 grid points, and the stem was modelled using 32 1D FEM-DEM elements. This kind of discretisation was found to be sufficient for the following results.

3.1.1 A single stem

Laminar flow of low Reynolds number Red = 50 was initially simulated for verification purposes. A single flexible stem was placed at the middle of the channel, while it was clamped to its bed. The Young’s modulus normalised by \( \rho u_{\tau }^{2} \) was taken as 5 × 105 in order to achieve a noticeable deflection with such a low flow velocity. The deformation of the stem is shown in Fig. 2 as calculated by the FSI code after the flow reached a steady-state condition. An analytical prediction based on beam theory is also plotted where the load on the stem was taken from the FSI solver. Both predictions show excellent agreement, verifying at least the FSI structural solver that is based on the FEM-DEM approach. As the rest of the solver was extensively verified in other references including channel flows [21, 22], we proceed to the results of a stem embedded in a turbulent channel flow.

Fig. 2
figure 2

Single stem’s streamwise deflection in a laminar channel flow

Full size image

Turbulent channel flow for Red = 250 was achieved by initiating the flow with random disturbance and waiting until a fully developed turbulent channel flow was reached [21, 22]. A single stem was then placed in the middle of the computational domain again and the simulation continued until the flow settled to a steady periodic state. The normalised Young’s modulus for the flexible stem was raised to 2 × 108 to account for the higher flow velocity while achieving a reasonable deflection. Typical contours of the instantaneous streamwise velocity are shown in Fig. 3 for the cases of a rigid stem and a flexible stem. A similar pattern is revealed for both cases where the highest velocity is near the top of the channel as expected and the effect of the stem is similar in both cases while only mildly causing any noticeable velocity reduction behind the stem. This is because of the stem’s small size and the associated low Reynolds number flow. The time development of the drag force acting on the two stems is shown in Fig. 4, where the drag force was normalised by \( 0.5\rho u_{\tau }^{2} A \) and A is the cross-sectional area of a virtual square enclosing the stem, i.e. (0.02 d)2. The time was normalised by d/uτ. The rigid stem is seen to experience a higher drag force, showing that vegetation can reduce the drag force acting on it through bending towards the direction of the flow.

Fig. 3
figure 3

Instantaneous streamwise velocity contours for a single a rigid stem and a b flexible stem

Full size image
Fig. 4
figure 4

Time history of the normalised drag force acting on the single stem

Full size image

3.1.2 A vegetation patch

A circular vegetation patch was placed at the middle of the computational domain instead of the single stem. The patch had a diameter of 0.5 d, and the stems were uniformly distributed inside it with a volumetric ratio of 0.1188. The normalised Young’s modulus for the flexible stem was reduced to 2 × 107. Typical instantaneous streamwise velocity contours are shown in Fig. 5. As in the single-stem case, the peak velocity is near the top of the channel. However, a more profound wake is seen behind the patch than behind the single stem in Fig. 3, and in that instant of time, the wake behind the rigid stems’ patch is longer than behind the flexible stems’ patch. One can also notice that while the upstream side of the flexible stems’ patch shows some deformation, the stems on the downstream side of the patch do not, hence further contributing towards the wake generation.

Fig. 5
figure 5

Instantaneous streamwise velocity contours for a circular patch of a rigid and b flexible stems

Full size image

In order to look at the time history of the patch, the time variation of the drag force felt by the patch is plotted in Fig. 6, where the drag force was normalised by the cross-sectional areas of the virtual squares enclosing the stems, i.e. 297*(0.02 d)2. Obviously, the drag forces for the rigid and flexible patches are not the same, but unlike in the case of the single stem, both forces show similar level when in some instant of time the drag force of the flexible patch is higher and sometimes the opposite. Increasing the normalised Young’s modulus back to 2 × 108 further reduced the difference between the two drag forces, making them almost identical. This is because the stems acted in a synchronised way strengthening each other as fibres of a rope strengthen each other.

Fig. 6
figure 6

Time history of the normalised drag force acting on the stems’ patch

Full size image

3.2 The horizontal axis kinetic turbine rotor

While the vegetation stem is dominated by the drag force, the blade of the horizontal axis kinetic turbine is dominated by the lift force. The current FSI code was already successfully used to investigate the flow around a three-blade horizontal axis marine current turbine (HAMCT), where the lift was generated due to the hydrodynamic profile of the blade [10]. However, lift can also be generated by a spinning cylinder as in the Magnus blade with an ability to achieve high lift coefficient, but at a cost of a reduced lift-to-drag ratio [26]. The Magnus blade is also a good problem to demonstrate the ability of the current FSI formulation to accurately capture the shear stress at the interface between the solid and the fluid, and hence simulate correctly the generated circulation in the flow, leading to the lift. Therefore, a single Magnus blade is investigated next, followed by a flexible three-blade rotor with the same geometry of the rigid rotor that was studied by Bai et al. [10].

3.2.1 The spinning cylinder

A rigid spinning cylinder was placed in a computational domain of [− 17.8, 47.8] × [− 32.8, 32.8] × [0, 0.5] in the streamwise, stream-normal and spanwise directions, respectively, where the cylinder’s diameter was normalised to one. The cylinder spanned all along the domain’s spanwise length and was placed at the origin of the streamwise-normal coordinates. The Reynolds number was set to ReD = 200 as relative to the free-stream velocity U and the cylinder’s diameter. Hence, the flow is laminar by its nature. The ratio between the tangential velocity of the cylinder’s circumference and U was set to 1.5. Dirichlet inflow boundary condition was used together with free-slip top and bottom boundary conditions. A Neumann-type boundary condition was adopted at the outflow, and periodicity was imposed in the spanwise direction to mimic a cylinder with an infinite length. The computational domain was discretised using a gradual stretched Cartesian mesh with a grid resolution of 384 × 384 × 16 with a gradual stretching mesh spacing in order to cluster grid points near the cylinder. Along the circumference of the cylinder, 256 immersed boundary points (IBPs) were evenly distributed in each x–y layer and 16 layers of IBPs in the spanwise direction were adopted.

Figure 7 shows the contours of the instantaneous pressure as relative to the ambient pressure and the vorticity component pointing spanwise, together with the streamlines. One should note that the pressure was normalised by ρU2. Clearly, a low-pressure region can be observed near the cylinder’s top surface, which attributes to the increase in lift as compared with a stationary cylinder. This is also evident by the vortex wake shedded downwards, thus causing a lift force acting on the cylinder. One should note that increasing the ratio of the tangential velocity to the free-stream velocity U from 1.5 to 2 will suppress the vortex shedding at this Reynolds number as was found by Mittal and Kumar [27].

Fig. 7
figure 7

Contours of a instantaneous pressure and b z-vorticity, together with streamlines flow past a spinning cylinder with ReD = 200 and tangential-to-free-steam velocity ratio of 1.5

Full size image

Table 1 compares the statistical features of the present results with those of Mittal and Kumar [27] who studied laminar flow past a spinning cylinder at various rotation rates using a body conformal (BC) grid with a very fine mesh near the cylinder. The Strouhal number St relates the vortex shedding main frequency, and hence, the oscillations of the force act on the cylinder. The lift and drag coefficients are normalised as relative to the ρU2A/2, where A is the cross section of the cylinder. While the mean lift coefficient shows excellent agreement between two results, a mild discrepancy is found in the averaged drag coefficient. Our mean drag coefficient is about 8% larger than of Mittal and Kumar [27].

Table 1 Comparison of drag coefficient, lift coefficient and Strouhal number for flow past a spinning circular cylinder at ReD = 200 and tangential-to-free-steam velocity ratio of 1.5
Full size table

The mildly over-predicted drag coefficient is mostly related to the side-wall effects. Mittal and Kumar [27] did a convergence study on the drag and lift coefficients by changing the wall–cylinder distance from 15 to 150 times of the cylinder’s diameter. They stated that the drag coefficient showed a close dependency on the cylinder’s length and the side-wall effects were negligible when the length is larger than 100 times the cylinder’s diameter, while the lift force was insensitive to the distance.

3.2.2 The three-blade HAMCT flexible rotor

The rigid HAMCT studied by Bai et al. [10] was further simulated, but by making the blades flexible and composed of a material with a Young’s modulus of 80 MPa, which is about 100 times smaller than the modulus of the standard steel in order to hasten the process of the blade deformation. The Poisson ratio was set as 0.3, and the density of the solid body was set as seven times of water. The computational domain size was set as [− 3, 7]D, [− 1.95, 1.05]D and [− 2.5, 2.5]D in the streamwise stream-normal and spanwise directions, respectively, where D is the diameter of the rotor. Inflow velocity condition and outflow condition were used in the streamwise direction, and free-slip walls were used on the bottom of the computational domain and the spanwise sides. Free surface was assumed at the top of the computational domain [10]. A stretched Cartesian mesh was used with a grid resolution of (640, 320, 512) points, while clustering points were near the rotor. The Reynolds number based on the free-stream velocity and the rotor’s diameter is 1.4 M.

The rigid HAMCT simulation result for tip speed ratio TSR = 6 was used as the initial condition, where TSR = ΩR/U, Ω is the rotational speed of the rotor and R is its radius. This means the Reynolds number based on the overall speed and chord length was about 200 k, and hence, the blade section is dominated by a laminar boundary layer going to transition through a laminar separation bubble on the suction side of the blade. Such process was well captured by the current FSI solver for rigid blades [10, 15]. The y + resolution on the blade was about 30 to 50.

The deformation of the blade at different time stages is shown in Fig. 8, where the blade in its initial rigid condition is shown in Fig. 8a. The blade starts to deform after about 20 rotations as demonstrated in Fig. 8b, where the highest stress is found at the root of the blade and near the trailing edge of the blade’s profile. The deformation process accelerates, resulting in a clear twist in the blade after 50 rotations as shown in Fig. 8c. After 100 rotations, the blade lost its shape and also bended opposite to the direction of the rotation. In that situation, the HAMCT power has been reduced by about 40% as compared to the rigid rotor’s produced power.

Fig. 8
figure 8

Blade deformation of HAMCT for a initial condition as of a rigid blade, b after 20 revolutions, c after 50 revolutions and d after 100 revolutions

Full size image

4 Conclusions

The presented FSI methodology coupled a FEM-DEM structural solver with a LES incompressible flow solver using an IBM algorithm. The FSI solver has been applied to submerged vegetation stems and small kinetic turbine blades embedded in laminar and turbulent flows. Both bodes have slender cylinder-like shapes, but the stem experiences a dominant drag force, while the blade experiences a dominant lift force due to its design. Verification cases of a single stem and a spinning Magnus blade in laminar flows were successfully pursued following previous verification studies which used this solver for bio-fluids and sediment problems. It was shown that the flexibly of the stem could reduce the drag force experienced by the stem, but when the stem was part of a dense patch, the patch could act as a group and thus considerably reduced the effect of the stem’s flexibility. A small horizontal axis marine turbine rotor was studied for the effect of its blade’s flexibility by significantly reducing the Young’s modulus as compared to steel. This resulted in a significant loss of power due to unwanted twisting and bending of the blades in less than 100 revolutions of the rotor.

This study has demonstrated the capability of a partitioned FSI solver based on IBM and FEM-DEM to deal with slender bodies. It has pointed to the growing importance of FSI in water flows as of vegetation and renewable energy flow systems as of the kinetic turbine to help in enhancing our understanding and ability to control and use such flows for the benefit of society.