Next Article in Journal
Isotopic Signatures as an Indicator of Long-Term Water-Use Efficiency of Haloxylon Plantations on the Dried Aral Sea Bed
Previous Article in Journal
Assessment of Effective Dose from Radioactive Isotopes Contained in Mineral Waters Received by Patients During Hydrotherapy Treatments
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China
2
Department of Mechanical Engineering, Michigan State University, East Lansing, MI 48824, USA
3
College of Energy and Electrical Engineering, Hohai University, Nanjing 210098, China
4
College of Agricultural Engineering, Hohai University, Nanjing 210098, China
*
Author to whom correspondence should be addressed.
Water 2020, 12(1), 98; https://doi.org/10.3390/w12010098
Submission received: 2 September 2019 / Revised: 23 December 2019 / Accepted: 23 December 2019 / Published: 27 December 2019
(This article belongs to the Section Hydraulics and Hydrodynamics)

Abstract

:
Operating in the harsh marine environment, fluctuating 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 fluid–solid interaction (FSI) for an MCT. The objective was to afford insights into the hydrodynamics near the rotor and in the wake, the deformation of rotor blades, and the interaction between the solid and fluid field. The numerical fluid simulation results showed good agreement with the experimental data and the influence of the support on the power coefficient 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 findings can inform the design of MCTs for enhancing robustness and survivability.

1. 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 could afford more detailed hydrodynamic information about the region of interest and accurately predict the performance of the turbine. The early famously employed methods are the actuator disc theory, which ignores the shape of rotor, and blade element momentum method, which ignores the span-wise flow. These simplifications greatly reduce the computational cost. The former one can well-predict the wake field far away from the turbine [9,10,11], and the latter provides a satisfactory results of turbine performance [12,13]. However, the flow structure around the rotor, particularly the tip vortices, and the load on the blade cannot be fully described [14,15].
Complex and strong turbulence near the MCT is caused by fluid–solid interaction, including wake vortex generation, dynamic stall, or instability of hydrodynamic loads. Recently, computational fluid dynamics (CFD) methods were widely employed that use 3D cells to represent explicitly the MCT’s geometry and motion. The fluid field is computed using Reynolds-averaged Navier–Stokes (RANS) or large eddy simulation (LES) as a turbulence model. RANS models are justifiably popular, as the computational requirements are affordable though the time-averaging of the velocity field does not realistically represent the time-varying fluctuations associated with blade-generated turbulence or instantaneous fluid–blade interaction [14,16]. By contrast, LES can resolve large-scale flow structures present in the velocity field [17,18]. Furthermore, increasing computational technology efficiently compensates the known drawbacks of LES, namely its high computational cost, especially for highly turbulent flows. This makes it more possible to have highly resolved dense meshes and obtain detailed flow information.
Consequently, some researchers have used LES to investigate the fluid dynamics and performance of MCTs. McNaughton et al. [19] applied sliding mesh LES method to model the flow and tip vortices of MCT, which were influenced by the support. Kang et al. [20] utilized curvilinear immersed boundary LES method to investigate the wake of MCT.
Most of previous studies have investigated the performance of MCTs and flow field with little consideration of structural response. To reduce the maintenance cost and risk, Singh et al. [21] transferred the blade surface pressure from the steady flow field simulation to the structural analysis code. Then, the stress, strain, and deformation were calculated to verify the safety of the blades.
Due to the unstable current and induced alternating loads, MCTs face the possibility of fatigue failure. Also, it is beneficial to consider the interaction between the structure and the flow field in evaluating the performance and the details of the structure response of MCTs. Badshah et al. [22] compared the one-way and two-way FSI method applied to an MCT modeled as a structural-steel solid. The maximum blade displacement was only 1/2000 of its radius. In addition, the RANS model was used, which cannot reflect the influence of the blade variation on the flow.
Bai et al. [23] presented an immersed boundary method and two free surface methods to simulate the MCT with FSI in a relatively small computational domain without consideration of the support or discussion of the blade deformation.
The object of this paper was to investigate the interaction between the fluid field and MCT with support in a two-way FSI, using LES. LES can more accurately solve large energetic scales in the flow field better than RANS [24], thus providing more details for the FSI simulation of the MCT rotor that is typically subjected to considerable turbulence, and in reverse for the impact on the performance of MCT and the flow field. The model was validated by comparing the predicted thrust and power coefficients with experimental measurements for a 0.8-m MCT presented in the work of Bahaj et al. [6].

2. Computational Equations and Method

2.1. LES Governing Equations

In a turbulent flow, eddies of various sizes are generated. Direct numerical simulation (DNS) could, theoretically, directly solve the whole spectrum of turbulent length scale. However, its cost is proportional to the cube of Reynolds number. For most realistic engineering problems concerning high Reynolds number flows, DNS is still too expensive.
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:
Φ ¯ ( x ) = D Φ ( x ) G ( x , x ) d x ,
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:
φ ¯ ( x ) = 1 V V φ ( x ) d x ,   x V ,
where V is the volume of a computational cell. The filter function, G ( x , x ) , implied here is then
G ( x , x ) = { 1 V   ,   x V   0 ,   x   o t h e r w i s e .
For incompressible flows, filtering the continuity and momentum equations, one obtains
u ¯ i x i = 0 ,
and
u ¯ i t + ( u ¯ i u ¯ j ) x j + τ i j x j = 1 ρ p ¯ x i + f ¯ i + σ i j x j ,
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   σ i j is the stress tensor due to molecular viscosity expressed by
σ i j υ ( u ¯ i x j + u ¯ j x i 2 3 u ¯ i x i δ i j ) ,
and τ i j is the subgrid-scale stress expressed by
τ i j u i u j ¯ u ¯ i u ¯ j ,
where v is the kinematic viscosity, δ i j 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
τ i j 1 3 τ k k δ i j = 2 v t S ¯ i j ,
where τ k k is the isotropic part of the subgrid-scale stresses, which is not modeled but added to the filtered static pressure term.   S ¯ i j is the rate-of-strain tensor for the resolved scale expressed by
S ¯ i j 1 2 ( u ¯ i x ¯ j + u ¯ j x ¯ i ) ,
and v t is the subgrid-scale turbulent viscosity. Here, the Wall-Adapting Local Eddy Viscosity (WALE) Model is adopted, it is modeled by
υ t = L S 2 ( S i j d S i j d ) 3 / 2 ( S ¯ i j S ¯ i j ) 5 / 2 + ( S i j d S i j d ) 5 / 4 ,
where L s and S i j d are defined, respectively, as
L s = m i n ( κ d , C w V 1 / 3 ) ,
S i j d = 1 2 ( g ¯ i j 2 + g ¯ i j 2 ) 1 3 δ i j g ¯ k k 2 ,   g ¯ i j = u ¯ i x j ,
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].

2.2. Solid Governing Equation

The conservation equation mainly used in solid solution is as follows:
ρ s d ¨ s = f s + · σ s ,
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
σ s = H s ε s ,
where, H s is the tensor of the elastic coefficients, ε s is the Cauchy’s strain tensor expressed by
ε s = 1 2 ( i d s + ( i d s ) T ) ,
where, d s is the solid displacement [26].

2.3. 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
d s 1 = d s 0 = d s t ,
τ f 0 = τ f t ,
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 f k is obtained from
F f [ X f k , λ d d s k 1 + ( 1 λ d ) d s k 2 ] = 0 ,
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 s k is obtained from
F s [ X s k , λ τ τ f k + ( 1 λ τ ) τ f k 1 ] = 0 ,
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
d f k = λ d d s k + ( 1 λ d d s k 1 ) .
(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.

3. Computational Domain and Verification

3.1. 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.2 R ( R is the rotor radius), while the pitch angle at the blade tip is 5°.

3.2. Computational Domain

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 × 106 cells. Figure 2a shows the surface of the coarse 3 × 106 cells mesh on the rotor, hub and support. Figure 2b gives an enlargement of a part of the blade surface mesh around r = 0.4 R . Figure 2c shows the blade surface mesh around r = 0.4 R   with a mesh three-times denser in all three directions, resulting in 81 × 106 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 × 106 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.4 R .

3.3. 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 × 106. Convergence was specified as RMS residual of 1 × 10−5. The main nondimensional load parameters were
power   coefficient :   C P = Q ω 1 2 ρ π R 2 U 0 3   ,
thrust   coefficient : C T = T 1 2 ρ π R 2 U 0 2   .
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.

4. Two-way FSI Numerical Simulation

4.1. 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 × 105 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.

4.2. Results and Discussion

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.
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/27 R , 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 7 shows the displacement in flow direction at the tip of each blade versus the angle of rotation. Points 4 and 5 were at tip of the other two blades at same location on the profile as Point 3 in Figure 6. Point 3 was on the tip of the blade, which was in front of the support at 0° in Figure 7 like shown in Figure 1 and Figure 2a. Point 5 was on the tip of the blade, which was marked by the ellipse in Figure 2a. Point 4 was on the tip of the blade left in Figure 2a. There was a delay of about 15° after each blade passes in front of the support until the displacement of the blade tip decayed to a minimum value. The plotted points show a periodic trend with a period of 360° and the scatter indicates additional higher-frequency oscillations of the blades throughout the rotation.
Figure 8 shows the displacement velocity in flow direction of Point 3 against the angle of rotation. The majority of points were concentrated between −0.2 m/s to 0.2 m/s. The blade oscillating in flow direction changed the relative velocity approaching the blade and hence the local velocity triangle. 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 Figure 8 and Figure 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 R . The root can be thickened to avoid damage. However, it is worth noting that the maximum stress occurred around 0.5 R 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].
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 ( Ω i j Ω i j S i j S i j ) , 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.

5. 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/27 R 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 R 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 × 108 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.

Author Contributions

Conceptualization, J.G. and F.C.; Data curation, J.G.; Formal analysis, J.G. and Y.Z.; Funding acquisition, J.G. and Y.Z.; Investigation, J.G.; Methodology, J.G. and H.C.; Project administration, F.C. and N.M.; Resources, N.M.; Software, J.G. and H.C.; Supervision, F.C.; Validation, J.G.; Visualization, J.G.; Writing—original draft, J.G.; Writing—review & editing, F.C. and N.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by Fundamental Research Funds for the Central Universities (2016B41514), Priority Academic Program Development of Jiangsu Higher Education Institutions (YS11001), National Natural Science Foundation of China (51809083), Natural Science Foundation of Jiangsu Province (BK20180504) and China Postdoctoral Science Foundation (2019M651678), China Scholarship Council (201706710022).

Acknowledgments

This work has been performed at the Turbomachinery Laboratory at Michigan State University with the support from Research Institute of Hydraulic and Hydropower Engineering at Hohai University.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Li, W.; Zhou, H.; Liu, H.; Lin, Y.; Xu, Q. Review on the blade design technologies of tidal current turbine. Renew. Sustain. Energy Rev. 2016, 63, 414–422. [Google Scholar] [CrossRef]
  2. Sun, C.; Lam, W.H.; Lam, S.S.; Dai, M.; Hamill, G. Temporal Evolution of Seabed Scour Induced by Darrieus-Type Tidal Current Turbine. Water 2019, 11, 896. [Google Scholar] [CrossRef] [Green Version]
  3. Ma, Y.; Hu, C.; Li, Y.; Li, L.; Deng, R.; Jiang, D. Hydrodynamic Performance Analysis of the Vertical Axis Twin-Rotor Tidal Current Turbine. Water 2018, 10, 1694. [Google Scholar] [CrossRef] [Green Version]
  4. Ahmed, U.; Apsley, D.D.; Afgan, I.; Stallard, T.; Stansby, P.K. Fluctuating loads on a tidal turbine due to velocity shear and turbulence: Comparison of CFD with field data. Renew. Energy 2017, 112, 235–246. [Google Scholar] [CrossRef] [Green Version]
  5. Coiro, D.P.; Maisto, U.; Scherillo, F.; Melone, S.; Grasso, F. Horizontal Axis Tidal Current Turbine: Numerical and Experimental Investigations. Owemes 2006, 20. Available online: https://www.researchgate.net/publication/237620241_Horizontal_Axis_Tidal_Current_Turbine_Numerical_and_Experimental_Investigations (accessed on 27 December 2019).
  6. Bahaj, A.S.; Molland, A.F.; Chaplin, J.R.; Batten, W.M.J. Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank. Renew. Energy 2007, 32, 407–426. [Google Scholar] [CrossRef]
  7. Song, M.; Kim, M.; Do, I.; Rhee, S.H.; Lee, J.H.; Hyun, B. Numerical and experimental investigation on the performance of three newly designed 100 kW-class tidal current turbines. Int. J. Nav. Arch. Ocean Eng. 2012, 4, 241–255. [Google Scholar] [CrossRef] [Green Version]
  8. Gaurier, B.; Davies, P.; Deuff, A.; Germain, G. Flume tank characterization of marine current turbine blade behaviour under current and wave loading. Renew. Energy 2013, 59, 1–12. [Google Scholar] [CrossRef] [Green Version]
  9. Sun, X. Numerical and Experimental Investigation of Tidal Current Energy Extraction, Doctor of Philosophy; The University of Edinburgh: Edinburgh, Scotland, 2008. [Google Scholar]
  10. Blackmore, T.; Batten, W.M.; Bahaj, A.S. Influence of turbulence on the wake of a marine current turbine simulator. Proc. R. Soc. A 2014, 470, 20140331. [Google Scholar] [CrossRef] [Green Version]
  11. Harrison, M.E.; Batten, W.M.J.; Myers, L.E.; Bahaj, A.S. Comparison between CFD simulations and experiments for predicting the far wake of horizontal axis tidal turbines. IET Renew. Power Gen. 2010, 4, 613. [Google Scholar] [CrossRef] [Green Version]
  12. Batten, W.M.J.; Bahaj, A.S.; Molland, A.F.; Chaplin, J.R. Experimentally validated numerical method for the hydrodynamic design of horizontal axis tidal turbines. Ocean Eng. 2007, 34, 1013–1020. [Google Scholar] [CrossRef]
  13. Bahaj, A.S.; Batten, W.M.J.; McCann, G. Experimental verifications of numerical predictions for the hydrodynamic performance of horizontal axis marine current turbines. Renew. Energy 2007, 32, 2479–2490. [Google Scholar] [CrossRef]
  14. Lin, X.; Zhang, J.; Zhang, Y.; Zhang, J.; Liu, S. Comparison of Actuator Line Method and Full Rotor Geometry Simulations of the Wake Field of a Tidal Stream Turbine. Water 2019, 11, 560. [Google Scholar] [CrossRef] [Green Version]
  15. Guo, Q.; Zhou, L.; Wang, Z. Comparison of BEM-CFD and full rotor geometry simulations for the performance and flow field of a marine current turbine. Renew. Energy 2015, 75, 640–648. [Google Scholar] [CrossRef]
  16. Zhang, Y.; Zhang, J.; Zheng, Y.; Yang, C.; Zang, W.; Fernandez-Rodriguez, E. Experimental Analysis and Evaluation of the Numerical Prediction of Wake Characteristics of Tidal Stream Turbine. Energies 2017, 10, 2057. [Google Scholar] [CrossRef] [Green Version]
  17. Stoesser, T. Large-eddy simulation in hydraulics: Quo Vadis? J. Hydraul. Res. 2014, 52, 441–452. [Google Scholar] [CrossRef]
  18. Visbal, M.; Yilmaz, T.O.; Rockwell, D. Three-dimensional vortex formation on a heaving low-aspect-ratio wing: Computations and experiments. J. Fluids Struct. 2013, 38, 58–76. [Google Scholar] [CrossRef]
  19. McNaughton, J.; Afgan, I.; Apsley, D.D.; Rolfo, S.; Stallard, T.; Stansby, P.K. A simple sliding-mesh interface procedure and its application to the CFD simulation of a tidal-stream turbine. Int. J. Numer. Methods Fluids 2014, 74, 250–269. [Google Scholar] [CrossRef]
  20. Kang, S.; Yang, X.; Sotiropoulos, F. On the onset of wake meandering for an axial flow turbine in a turbulent open channel flow. J. Fluids Mech. 2014, 744, 376–403. [Google Scholar] [CrossRef]
  21. Singh, P.M.; Choi, Y. Shape design and numerical analysis on a 1 MW tidal current turbine for the south-western coast of Korea. Renew. Energy 2014, 68, 485–493. [Google Scholar] [CrossRef]
  22. Badshah, M.; Badshah, S.; Kadir, K. Fluid Structure Interaction Modelling of Tidal Turbine Performance and Structural Loads in a Velocity Shear Environment. Energies 2018, 11, 1837. [Google Scholar] [CrossRef] [Green Version]
  23. Bai, X.; Avital, E.J.; Munjiza, A.; Williams, J.J.R. Numerical simulation of a marine current turbine in free surface flow. Renew. Energy 2014, 63, 715–723. [Google Scholar] [CrossRef]
  24. Ouro, P.; Harrold, M.; Stoesser, T.; Bromley, P. Hydrodynamic loadings on a horizontal axis tidal turbine prototype. J. Fluids Struct. 2017, 71, 78–95. [Google Scholar] [CrossRef] [Green Version]
  25. ANSYS Inc. Guide, Fluent Release 19.1; ANSYS Inc.: Canonsburg, PA, USA, 2019. [Google Scholar]
  26. Sarrate, J.; Huerta, A.; Donea, J. Arbitrary Lagrangian-Eulerian formulation for fluid-rigid body interaction. Comput. Method Appl. Mech. 2001, 190, 3171–3188. [Google Scholar] [CrossRef] [Green Version]
  27. Bai, G.; Li, G.; Ye, Y.; Gao, T. Numerical analysis of the hydrodynamic performance and wake field of a horizontal axis tidal current turbine using an actuator surface model. Ocean Eng. 2015, 94, 1–9. [Google Scholar] [CrossRef]
  28. Park, S.; Park, S.; Rhee, S.H. Influence of blade deformation and yawed inflow on performance of a horizontal axis tidal stream turbine. Renew. Energy 2016, 92, 321–332. [Google Scholar] [CrossRef]
  29. Afgan, I.; McNaughton, J.; Rolfo, S.; Apsley, D.D.; Stallard, T.; Stansby, P. Turbulent flow and loading on a tidal stream turbine by LES and RANS. Int. J. Heat Fluid Flow 2013, 43, 96–108. [Google Scholar] [CrossRef]
  30. Nicholls-Lee, R.F.; Turnock, S.R.; Boyd, S.W. Application of bend-twist coupled blades for horizontal axis tidal turbines. Renew. Energy 2013, 50, 541–550. [Google Scholar] [CrossRef]
Figure 1. Computation domain discretized with 3 × 106 hexahedral cells.
Figure 1. Computation domain discretized with 3 × 106 hexahedral cells.
Water 12 00098 g001
Figure 2. (a) The surface mesh of the rotor, hub and support in coarse mesh; (b) The blade surface mesh around r = 0.4 R ; (c) The blade surface mesh around r = 0.4 R densified 3 times in all 3 directions; (d) The blade surface mesh acturally used in the simulation; (e) The section at r = 0.4 R with adapted mesh.
Figure 2. (a) The surface mesh of the rotor, hub and support in coarse mesh; (b) The blade surface mesh around r = 0.4 R ; (c) The blade surface mesh around r = 0.4 R densified 3 times in all 3 directions; (d) The blade surface mesh acturally used in the simulation; (e) The section at r = 0.4 R with adapted mesh.
Water 12 00098 g002
Figure 3. Comparisons between numerical and experimental data from Bahaj et al. [6] at different TSR. (a) power coefficient and (b) thrust coefficient
Figure 3. Comparisons between numerical and experimental data from Bahaj et al. [6] at different TSR. (a) power coefficient and (b) thrust coefficient
Water 12 00098 g003
Figure 4. The surface mesh of the solid domain of the rotor.
Figure 4. The surface mesh of the solid domain of the rotor.
Water 12 00098 g004
Figure 5. (a) Power coefficient and (b) thrust coefficient versus angle of rotation obtained from CFD only and FSI.
Figure 5. (a) Power coefficient and (b) thrust coefficient versus angle of rotation obtained from CFD only and FSI.
Water 12 00098 g005
Figure 6. Total displacement of the blade at 207°.
Figure 6. Total displacement of the blade at 207°.
Water 12 00098 g006
Figure 7. Displacement in flow direction versus the angle of rotation for the tip of each blade (Points 3–5).
Figure 7. Displacement in flow direction versus the angle of rotation for the tip of each blade (Points 3–5).
Water 12 00098 g007
Figure 8. Displacement velocity in flow direction of the Point 3 versus the angle of rotation.
Figure 8. Displacement velocity in flow direction of the Point 3 versus the angle of rotation.
Water 12 00098 g008
Figure 9. Local velocity triangle of blade tip section.
Figure 9. Local velocity triangle of blade tip section.
Water 12 00098 g009
Figure 10. Flow angle θ of the Point 3 versus the angle of rotation.
Figure 10. Flow angle θ of the Point 3 versus the angle of rotation.
Water 12 00098 g010
Figure 11. Contour of equivalent (von Mises) stress on blade suction side at 207°.
Figure 11. Contour of equivalent (von Mises) stress on blade suction side at 207°.
Water 12 00098 g011
Figure 12. Equivalent (von Mises) stress on the root of the three blades normalized by 130 MPa versus the angle of rotation.
Figure 12. Equivalent (von Mises) stress on the root of the three blades normalized by 130 MPa versus the angle of rotation.
Water 12 00098 g012
Figure 13. Iso-Q surfaces colored by velocity in flow direction.
Figure 13. Iso-Q surfaces colored by velocity in flow direction.
Water 12 00098 g013
Figure 14. The eddy viscosity in the central cut plane.
Figure 14. The eddy viscosity in the central cut plane.
Water 12 00098 g014
Table 1. Mesh details and sensitivity study for mean blade loading at TSR = 5.67.
Table 1. Mesh details and sensitivity study for mean blade loading at TSR = 5.67.
MeshNumber of CellsTime Step (s)Y PlusCPCT
Coarse3 × 1061 × 10−41.2–90.3930.707
Medium28 × 1061 × 10−40.8–60.420.761
Fine100 × 1061 × 10−40.4–30.4310.763
Fine100 × 1061 × 10−50.4–30.4310.763
Table 2. Properties of T6082-T6 aluminum alloy.
Table 2. Properties of T6082-T6 aluminum alloy.
ParametersValues
Density2700 kg/m3
Young’s Modulus70 GPa
Poisson’s Ratio0.33
Tensile yield strength270 MPa
Tensile ultimate strength330 MPa

Share and Cite

MDPI and ACS Style

Gu, J.; Cai, F.; Müller, N.; Zhang, Y.; Chen, H. Two-Way Fluid–Solid Interaction Analysis for a Horizontal Axis Marine Current Turbine with LES. Water 2020, 12, 98. https://doi.org/10.3390/w12010098

AMA Style

Gu J, Cai F, Müller N, Zhang Y, Chen H. Two-Way Fluid–Solid Interaction Analysis for a Horizontal Axis Marine Current Turbine with LES. Water. 2020; 12(1):98. https://doi.org/10.3390/w12010098

Chicago/Turabian Style

Gu, Jintong, Fulin Cai, Norbert Müller, Yuquan Zhang, and Huixiang Chen. 2020. "Two-Way Fluid–Solid Interaction Analysis for a Horizontal Axis Marine Current Turbine with LES" Water 12, no. 1: 98. https://doi.org/10.3390/w12010098

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop