Two-Way Fluid–Solid Interaction Analysis for a Horizontal Axis Marine Current Turbine with LES

: Operating in the harsh marine environment, ﬂuctuating loads due to the surrounding turbulence are important for fatigue analysis of marine current turbines (MCTs). The large eddy simulation (LES) method was implemented to analyze the two-way ﬂuid–solid interaction (FSI) for an MCT. The objective was to a ﬀ ord insights into the hydrodynamics near the rotor and in the wake, the deformation of rotor blades, and the interaction between the solid and ﬂuid ﬁeld. The numerical ﬂuid simulation results showed good agreement with the experimental data and the inﬂuence of the support on the power coe ﬃ cient and blade vibration. The impact of the blade displacement on the MCT performance was quantitatively analyzed. Besides the root, the highest stress was located near the middle of the blade. The ﬁndings can inform the design of MCTs for enhancing robustness and survivability. Formal analysis, J.G. and Y.Z.; Funding J.G. and Y.Z.; Investigation, J.G.; Methodology, J.G. and H.C.; F.C. and N.M.; Resources, N.M.; Software, J.G. and H.C.; Supervision, F.C.; Validation, J.G.; Visualization, J.G.; J.G.;


Introduction
As a kind of clean energy with huge reserves, marine current can be predicted and utilized accurately [1][2][3]. The marine current turbine (MCT) was adapted first from wind turbines, and considerable part of the technology can be directly applied. However, the density of water is about 800-times higher than that of air. Compared with wind turbines, water loads on MCTs are much larger. Besides, there are more challenges, such as biofouling, marine debris, and cavitation [4].
In the early stage, the experimental research on MCTs mainly focused on fluid dynamics and performance. Coiro et al. [5] obtained the characteristics of horizontal axis tidal current turbine power and thrust for a range of tip speed ratio (TSR) and hub pitch angle through towing tank tests. Bahaj et al. [6] investigated the effect of flow speed, rotation speed and hub pitch angle on MCT performance in a towing tank and a cavitation tunnel. They also studied the secondary effects of rotor depth, yaw and dual rotor interference. Song et al. [7] focused on the tip vortex generated by different raked blade tip shapes to optimize the performance of MCTs and verified their results by towing tank experiments with models. Gaurier et al. [8] studied the influence of wave-current interactions on the MCT performance and measured strains of blades in a flume tank.
An increasing number of researchers are interested in developing numerical methods that could simulate complex flows near MCTs. Those are also useful for field or laboratory tests because they Mass, momentum, and energy are transported mainly by large eddies, which are dictated by boundary conditions and geometries within the involved flow. Small eddies depend less on the geometries. With LES, larger eddies that can be captured by cells are resolved directly, whereas the effects of smaller eddies that cannot be captured by cells are accounted for with modeling, indirectly, using a sub grid model. This involves a filtering process, which makes it possible to simulate the transient flow with high-performance computers, although it still needs sufficient space and time resolution. A filtered variable (denoted by an overbar, such as Φ) at a point x is defined by the following convolution: where D is the fluid domain, and G is the filter function determining the scale of the resolved eddies. The finite-volume discretization itself implicitly affords the filtering procedure: where V is the volume of a computational cell. The filter function, G(x, x ), implied here is then For incompressible flows, filtering the continuity and momentum equations, one obtains and where u is the velocity vector, i and j range from 1 to 3 to denote three directions and velocity components respectively, ρ is the density, f is a volume forcing term, and σ ij is the stress tensor due to molecular viscosity expressed by and τ ij is the subgrid-scale stress expressed by where v is the kinematic viscosity, δ ij is Kronecker delta, which is 1 if i = j, and 0 otherwise. The subgrid-scale stress is unknown and requires modeling. Here, the Boussinesq hypothesis is adopted to compute the subgrid-scale stress where τ kk is the isotropic part of the subgrid-scale stresses, which is not modeled but added to the filtered static pressure term. S ij is the rate-of-strain tensor for the resolved scale expressed by and v t is the subgrid-scale turbulent viscosity. Here, the Wall-Adapting Local Eddy Viscosity (WALE) Model is adopted, it is modeled by where L s and S d ij are defined, respectively, as where κ is the von Kármán constant, d is the distance to the nearest wall, and C w , the WALE constant, is 0.325, which has produced consistently superior results under intensive validation [25].

Solid Governing Equation
The conservation equation mainly used in solid solution is as follows: where, ρ s is the solid density, .. d s is the local acceleration vector, f s is the body force, including the centrifugal force and fluid field pressure, σ s is the Cauchy's stress tensor expressed by where, H s is the tensor of the elastic coefficients, ε s is the Cauchy's strain tensor expressed by where, d s is the solid displacement [26].

Two-Way FSI Iterative Computation
In two-way FSI analysis, the fluid pressure is applied to the solid, where, in turn, the deformation of the solid affects the flow field. The computational domain are generally divided into fluid field and solid field, both of which are defined by materials, boundaries, and initial conditions. The interaction occurs at the interface between two fields. Two-way FSI iterative computation is also called partition method because the fluid equation and the solid equation are solved separately and continuously, and always update the latest information provided by the other part of the coupling system. This iteration continues until the solution of the coupling equation converges.
The calculation process could be summarized as: getting the solutions at t + ∆t by iterating between the fluid model and the solid model. Assume the initial solution where, τ f is the fluid stress. With iteration k = 1, 2, . . ., the equilibrium solution for each time step X t+∆t is obtained. (1) The fluid solution vector X k f is obtained from where, F f represents the correlative fluid formulas of the terms in the brackets. In fluid analysis, the solution is obtained by specifying the solid displacement. Since the fluid and the solid are not solved in the same matrix, the solid displacement was solved using the relaxation factor 0 < λ d < 1, which is conducive to iterative convergence.
(2) The stress residuals are calculated and compared with the tolerance due to the need to meet the stress criteria. When the criteria are met, steps 3-5 are skipped.
(3) The solid solution vector X k s is obtained from where, F s represents the correlative solid formulas of the terms in the brackets. Likewise, stress relaxation factor λ τ (0 < λ τ < 1) is used in fluid stress.
(4) The fluid node displacement is determined by the specified boundary conditions (5) Similarly, the displacement standard needs to be met, and the displacement residual is calculated and compared with the tolerance. If there is no convergence, the process returns to step 1 and iterates again.
(6) The fluid and solid solutions is saved.

Test Case Description
Experimental data is available for turbine performance characteristics in terms of power and thrust coefficients for ranges of TSR and blade pitch settings. In this case, the data from the towing tank tests of Bahaj s group at the University of Southampton [6] were modeled because the geometry of the MCT was provided along with data for model verification. Although there was a lack of experimental data for the velocity distribution in the wake or for the pressure distribution along the blades, additional CFD simulations were conducted to study the sensitivity of these parameters to TSR and inflow turbulence. Calculations were carried out with both RANS [15,27,28] and LES models [19,29].
The MCT tested in [6] and modeled within this research was equipped with three blades and a rotor with diameter of 0.8 m. The diameter of the nacelle and of the support was 0.1 m. The blades were developed from NACA 63-8xx profiles, with chord, thickness, and pitch distributions given at 17 locations along the blade. The towing tank had a depth and width of 1.8 m and 3.7 m, respectively, resulting in a blockage ratio of 7.5%. The rotor was centered 0.96 m above the bed. Among the results of the power coefficient (C P ) and thrust coefficient (C T ) from the tank test, the MCT with 20 • hub pitch angle at zero yaw angle has the best performance. The 20 • hub pitch angle refers to the nose-tail angle of the blade at a radius of 0.2R (R is the rotor radius), while the pitch angle at the blade tip is 5 • . Figure 1 illustrates the computation domain as it extends a distance of 5D upstream and 15D downstream from the point of origin in the rotor, where the rotational axes for the blade pitch cross. The domain was divided into an external domain and an internal domain of the rotor with sliding mesh [16]. The internal domain had a diameter of 1.5D (where D is the diameter of the rotor) and extended to a distance of 1/4D upstream and 1/8D downstream from the point of origin in the rotor. For visualization purposes, only in Figure 1, the hexahedral mesh is given as discretized with 3 × 10 6 cells. Figure 2a shows the surface of the coarse 3 × 10 6 cells mesh on the rotor, hub and support. Figure 2b gives an enlargement of a part of the blade surface mesh around r = 0.4R. Figure 2c shows the blade surface mesh around r = 0.4R with a mesh three-times denser in all three directions, resulting in 81 × 10 6 cells. Because LES performance is improved with a finer mesh for resolving even smaller eddies, the actual simulation was performed with a mesh of 100 × 10 6 cells as shown in Figure 2d,e, where the wall-normal and wall-parallel (streamwise and spanwise) spacing was decreased in all three directions as a wall was approached (Figure 2e). The adaptive mesh method was used via HyperMesh to avoid high near-wall aspect ratio, which remained between 1 and 6. Finally, the blade surface mesh consisted of 145 nodes in the chordwise and 756 nodes in the spanwise directions. The y plus was about 0.4 to 3 at r = 0.4R.

Verification of the Model
Before the FSI simulation, the fluid model was verified. Fluent 19.1 in ANSYS workbench was

Verification of the Model
Before the FSI simulation, the fluid model was verified. Fluent 19.1 in ANSYS workbench was

Verification of the Model
Before the FSI simulation, the fluid model was verified. Fluent 19.1 in ANSYS workbench was used for the CFD simulation. The above-described LES turbulence model was applied. The top of the computational domain was set as a slip wall, and the inlet velocity to the computational domain was uniformly U 0 = 1.4 m/s. The outflow boundary condition was used TSR = ωR/U 0 was varied in the range between 4 and 12 by changing rotational speed ω. The Reynolds number, based on the inlet velocity and diameter, was 1.12 × 10 6 . Convergence was specified as RMS residual of 1 × 10 −5 . The main nondimensional load parameters were power coefficient : thrust coefficient : where, Q is the torque, T is the thrust. Bahaj et al. [6] reported experimental values of C P , C T , and TSR with a blockage correction applied to identify them as 'equivalent open-water' measurements. The same correction was applied to the numerical results here, with the thrust coefficient calculated from numerical results. With the blockage ratio 7.5%, a typical C T = 0.8 resulted in multiplicative blockage correction factors of 92.2%, 94.7%, and 97.3% for the power coefficient, thrust coefficient, and TSR, respectively.
A complete mesh sensitivity study was performed with different mesh resolutions (see Table 1). The LES was grid-dependent, i.e., the finer the mesh, the more length-scales of turbulence were resolved explicitly. Hence, the fine mesh was chosen and time step size was set as 1 × 10 −5 s.  Figure 3 compares time averaged blockage-corrected power and thrust coefficients from pure CFD and FSI against the experimental data [6]. At each TSR, the simulation results were averaged over four rotations after the power and thrust coefficients achieved a repeatable periodic behavior. The numerical results both agree well with the experimental data. The only exception is the highest TSR, where the pure CFD and FSI underpredicted the C P by about 40% and 30%, respectively. The predicted C P results show the similar trend with the fitting curve of experimental data. With the increase of TSR, the C P increased to a maximum value and then decreased. Hence, the fluid model was deemed to be acceptable for further use in the numerical FSI simulation analysis. The C P and C T results at three different TSRs from FSI simulation are presented here in advance. It can be seen that the C P and C T at TSR = 4.9, 8.3 and 11.3 predicted by FSI are closer to the experimental data, while the C P and C T at TSR = 5.67 predicted by pure CFD are closer to the experimental data.
predicted results show the similar trend with the fitting curve of experimental data. With the increase of TSR, the increased to a maximum value and then decreased. Hence, the fluid model was deemed to be acceptable for further use in the numerical FSI simulation analysis. The and results at three different TSRs from FSI simulation are presented here in advance. It can be seen that the and at TSR = 4.9, 8.3 and 11.3 predicted by FSI are closer to the experimental data, while the and at TSR = 5.67 predicted by pure CFD are closer to the experimental data.

Two-Way FSI Model
The two-way FSI numerical simulation was also carried out in ANSYS workbench, where the fluid domain and the solid domain were established separately. The fluid domain was the computational domain described above, with the rotor surface set as the FSI interface to transmit the pressure information and receive the deformation information of the rotor. Because the blades will always be deformed, dynamic mesh method was adopted to avoid cells with negative volumes.
To facilitate the analysis of interaction details, only the rotor deformation was simulated, whereas the deformation of the support was left for further studies in the future. Due to its long-time operation in seawater, the MCT is generally manufactured from composite materials. However, the blades in the tank test are solid bodies manufactured from T6082-T6 aluminum alloy [6]. To be consistent with the experiment, aluminum alloy was chosen as the rotor material in this study. The properties of T6082-T6 aluminum alloy are provided in Table 2. The trailing edge of blade was very thin, which was difficult to discretize into hexahedral cells. Therefore, and to reduce the error caused by the information exchange and interpolation of nodes on the FSI interface, the quadrilateral cells on the rotor surface were extracted from the fluid domain and split into triangular cells with a diagonal line. Then, the tetrahedral cells in the rotor interior were generated freely, obtaining the solid domain with 2.5 × 10 5 cells. For visualization purpose, the cell size in Figure 4 was the same as in Figure 2a. Like the surface of the fluid domain on the rotor, the rotor surface of the solid domain was set to FSI interface. Fluent and Transient Structure are mutually iterative computations through a two-way coupled solver. TSR = 5.67 was chosen in the FSI simulation, i.e., ω = 20.43 rad/s before the blockage correction. These simulations were performed on high-performance computers using 128 processors, with a single FSI simulation taking three months and a single LES taking one month.    Figure 5 presents the power coefficient and thrust coefficient versus the angle of rotation during the fifth to seventh rotation obtained from CFD simulation and FSI simulation at the same TSR, respectively. Both plots reveal that the C P and C T obtained from pure CFD decayed to a minimum value every 120 • corresponding to a blade passing by the support. The support affected the flow field of the wake and reduced C P and C T by about 1.5% and 1.7% from the maximum value every time a blade passed by.  Figure 5 presents the power coefficient and thrust coefficient versus the angle of rotation during the fifth to seventh rotation obtained from CFD simulation and FSI simulation at the same TSR, respectively. Both plots reveal that the and obtained from pure CFD decayed to a minimum value every 120° corresponding to a blade passing by the support. The support affected the flow field of the wake and reduced and by about 1.5% and 1.7% from the maximum value every time a blade passed by. The and values obtained from FSI were about 2.4% and 1.6% lower than those from CFD due to the displacement of the blade. However, as mentioned above, in Figure 3, it is worth noting that the FSI did not always underpredict the compared to the pure CFD. Figure 6 illustrates the blade total displacement at 207°. The maximum total displacement is 1/27 , at the tip of the blade (Point 3). For each section of the blade, the total displacement of the nose and tail (e.g., Points 1 and 2) relative to each other was less than 1 × 10 −5 m, indicating that the twist of the blade can be viewed as negligibly small. The C P and C T values obtained from FSI were about 2.4% and 1.6% lower than those from CFD due to the displacement of the blade. However, as mentioned above, in Figure 3, it is worth noting that the FSI did not always underpredict the C P compared to the pure CFD. Figure 6 illustrates the blade total displacement at 207 • . The maximum total displacement is 1/27R, at the tip of the blade (Point 3). For each section of the blade, the total displacement of the nose and tail (e.g., Points 1 and 2) relative to each other was less than 1 × 10 −5 m, indicating that the twist of the blade can be viewed as negligibly small. due to the displacement of the blade. However, as mentioned above, in Figure 3, it is worth noting that the FSI did not always underpredict the compared to the pure CFD. Figure 6 illustrates the blade total displacement at 207°. The maximum total displacement is 1/27 , at the tip of the blade (Point 3). For each section of the blade, the total displacement of the nose and tail (e.g., Points 1 and 2) relative to each other was less than 1 × 10 −5 m, indicating that the twist of the blade can be viewed as negligibly small.       Figure 9 illustrates a local velocity triangle at the blade tip. It is well-known that the approaching flow velocity is less than the initial free-stream velocity In Figure 10, = was applied, and the tangential and radial components were not considered. For the rigid non-displaced blade, the inflow velocity angle = arctan = 9.72°, the pitch angle = 5°, and the angle of attack    Figure 9 illustrates a local velocity triangle at the blade tip. It is well-known that the approaching flow velocity u is less than the initial free-stream velocity U 0 In Figure 10, u = U 0 was applied, and the tangential and radial components were not considered. For the rigid non-displaced blade, the inflow velocity angle θ = arctan u ωr = 9.72 • , the pitch angle β = 5 • , and the angle of attack α = 4.72 • . To obtain the flow angle θ versus the angle of rotation in Figure 10, θ = arctan (u+∆u) ωr was used, where ∆u represents the displacement velocity in flow direction of Point 3 given in Figure 8. The majority of the points in Figure 10 fell in the range of 8 • to 11 • , indicating the high-frequency change of inflow velocity angle causing the higher-frequency fluctuations of C P predicted by FSI simulation in Figure 5. Cyclic behavior cannot be observed in both Figures 8 and 10. It is suspected that the flow separation was responsible for this phenomenon, which require further investigation.      Figure 11 illustrates the contour of equivalent (von Mises) stress on blade suction side, where the stresses were higher than on the pressure side. At 207°, the maximum value of 130 MPa occurred on the root in the Point 6. It was less than the tensile yield strength of 270 MPa and tensile ultimate strength of 330 MPa. Beyond the root, a maximum value of 81 MPa occurred near 0.5 . The root can be thickened to avoid damage. However, it is worth noting that the maximum stress occurred around 0.5 beyond the root at every moment, in the range where blades were broken [1].   Figure 12 shows the equivalent (von Mises) stress normalized by 130 MPa at the root of each blade versus the angle of rotation, where Points 7 and 8 were at the same location as Point 6 and on the same blade as Points 4 and 5, respectively. All blades experienced a maximum stress variation of 5.8% during one rotation due to the support. Bahaj et al. considered the influence of the support on MCT performance and placed it as far as at least one rotor radius downstream of the rotor. The variation would be larger a velocity profile of 1/7th power law was used instead of uniform velocity [22].

Results and Discussion
Relative to the time average, LES can not only predict the effect of the support on the performance of the MCT, but also capture the tip vortices disturbed by the support, as shown in Figure 13. Q-criterion was used: = 0.5( − ), where is the vorticity tensor and is the rate-of-strain tensor. The helicoidal tip vortices were generated by the fast-moving blade tip surface and then transported downstream. These helical structures closely attached the shear layer, rolling up between the fast flow right above the tip of each blade and the low-velocity rotor wake. The radius, at which the tip vortices traveled downstream in the wake, expanded with the radius of the sheer layer and increased beyond the radius of the rotor. Figure 13 also shows an additional three vortices, which originated from the root of the blades. They flowed downstream around the hub and were maintained by the sheer layer roll-up between the low-velocity hub boundary layer and the higher-velocity wake of the rotor. The root vortices rotated in the direction opposite to the tip vortices. Figure 14 shows the eddy viscosity in the central cut plane. The high values correspond to vortices in Figure 13.    Figure 12 shows the equivalent (von Mises) stress normalized by 130 MPa at the root of each blade versus the angle of rotation, where Points 7 and 8 were at the same location as Point 6 and on the same blade as Points 4 and 5, respectively. All blades experienced a maximum stress variation of 5.8% during one rotation due to the support. Bahaj et al. considered the influence of the support on MCT performance and placed it as far as at least one rotor radius downstream of the rotor. The variation would be larger a velocity profile of 1/7th power law was used instead of uniform velocity [22].
Relative to the time average, LES can not only predict the effect of the support on the performance of the MCT, but also capture the tip vortices disturbed by the support, as shown in Figure 13. Q-criterion was used: = 0.5( − ), where is the vorticity tensor and is the rate-of-strain tensor. The helicoidal tip vortices were generated by the fast-moving blade tip surface and then transported downstream. These helical structures closely attached the shear layer, rolling up between the fast flow right above the tip of each blade and the low-velocity rotor wake. The radius, at which the tip vortices traveled downstream in the wake, expanded with the radius of the sheer layer and increased beyond the radius of the rotor. Figure 13 also shows an additional three vortices, which originated from the root of the blades. They flowed downstream around the hub and  Figure 12 shows the equivalent (von Mises) stress normalized by 130 MPa at the root of each blade versus the angle of rotation, where Points 7 and 8 were at the same location as Point 6 and on the same blade as Points 4 and 5, respectively. All blades experienced a maximum stress variation of 5.8% during one rotation due to the support. Bahaj et al. considered the influence of the support on MCT performance and placed it as far as at least one rotor radius downstream of the rotor. The variation would be larger a velocity profile of 1/7th power law was used instead of uniform velocity [22].
Relative to the time average, LES can not only predict the effect of the support on the performance of the MCT, but also capture the tip vortices disturbed by the support, as shown in Figure 13. Q-criterion was used: Q = 0.5 Ω ij Ω ij − S ij S ij , where Ω is the vorticity tensor and S is the rate-of-strain tensor. The helicoidal tip vortices were generated by the fast-moving blade tip surface and then transported downstream. These helical structures closely attached the shear layer, rolling up between the fast flow right above the tip of each blade and the low-velocity rotor wake. The radius, at which the tip vortices traveled downstream in the wake, expanded with the radius of the sheer layer and increased beyond the radius of the rotor. Figure 13 also shows an additional three vortices, which originated from the root of the blades. They flowed downstream around the hub and were maintained by the sheer layer roll-up between the low-velocity hub boundary layer and the higher-velocity wake of the rotor. The root vortices rotated in the direction opposite to the tip vortices. Figure 14 shows the eddy viscosity in the central cut plane. The high values correspond to vortices in Figure 13.

Conclusions
As a form of ocean energy conversion, MCT has been widely accepted by the industry and researchers. The most concern is reducing design, installation, and maintenance costs. Based on this purpose, a two-way FSI method, suitable for MCTs, was described and implemented. LES was chosen instead of the time-averaged model to capture accurate pressure information from the flow field at any time, blade vibration, passing by the support and the tip vortices, which expand outward.
Bahaj's group expected the reductions of the and of the MCT every time a blade passed by the support to be very small. The results of this numerical simulation quantified these reductions at 1.5% and 1.7%, respectively. The oscillation and vibration of the blade reduced the and by an additional 2.4% and 1.6%, respectively. The maximum total displacement was 1/27 at the tip of the blade. With the blade displacement, the displacement velocity of the blade tip changed the velocity triangle continuously and periodically.
The maximum equivalent (von Mises) stress occurred on the root on the suction side, and a high value occurred near 0.5 at every moment. During every rotation, the stresses varied around 5.8% for each blade, which would be larger if the profiled inflow was considered. This can considerably impact the fatigue loading with an MCT, typically experiencing in the order of 1 × 10 8 rotational cycles over a 20-year life span [30].
Compared with pure CFD simulation without considering solid deformation and one-way fluid-solid interaction simulation, two-way fluid-solid interaction simulation can well-predict the stress variation experienced by blades under more complex loads.
In the future, we expect to explore the impact of the deformation of the support and the influence of a free surface and waves on the performance, and to simulate the MCT in a velocity shear environment.

Conclusions
As a form of ocean energy conversion, MCT has been widely accepted by the industry and researchers. The most concern is reducing design, installation, and maintenance costs. Based on this purpose, a two-way FSI method, suitable for MCTs, was described and implemented. LES was chosen instead of the time-averaged model to capture accurate pressure information from the flow field at any time, blade vibration, passing by the support and the tip vortices, which expand outward.
Bahaj's group expected the reductions of the C P and C T of the MCT every time a blade passed by the support to be very small. The results of this numerical simulation quantified these reductions at 1.5% and 1.7%, respectively. The oscillation and vibration of the blade reduced the C P and C T by an additional 2.4% and 1.6%, respectively. The maximum total displacement was 1/27R at the tip of the blade. With the blade displacement, the displacement velocity of the blade tip changed the velocity triangle continuously and periodically.
The maximum equivalent (von Mises) stress occurred on the root on the suction side, and a high value occurred near 0.5R at every moment. During every rotation, the stresses varied around 5.8% for each blade, which would be larger if the profiled inflow was considered. This can considerably impact the fatigue loading with an MCT, typically experiencing in the order of 1 × 10 8 rotational cycles over a 20-year life span [30].
Compared with pure CFD simulation without considering solid deformation and one-way fluid-solid interaction simulation, two-way fluid-solid interaction simulation can well-predict the stress variation experienced by blades under more complex loads.
In the future, we expect to explore the impact of the deformation of the support and the influence of a free surface and waves on the performance, and to simulate the MCT in a velocity shear environment.