Assessment of the Potential and Current Limitations of Integrating CFD and Structural Analyses for the Design of Offshore Wind Turbines Master’s thesis in Applied Mechanics Julius Kandt Department of Mechanics and Maritime Sciences CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2019 Master’s thesis 2019:25 Assessment of the Potential and Current Limitations of Integrating CFD and Structural Analyses for the Design of Offshore Wind Turbines Julius Kandt Department of Mechanics and Maritime Sciences Division of Fluid Dynamics Chalmers University of Technology Gothenburg, Sweden Assessment of the Potential and Current Limitations of Integrating CFD and Structural Analyses for the Design of Offshore Wind Turbines © Julius Kandt, 2019. Supervisor: Hamidreza Abedi, Division of Fluid Dynamics Examiner: Lars Davidson, Division of Fluid Dynamics Master’s Thesis 2019:25 Department of Mechanics and Maritime Sciences Division of Fluid Dynamics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Visualization of a FAST simulation of an offshore wind turbine with OC3 Tripod configuration. Typeset in LATEX Gothenburg, Sweden 2019 iv Assessment of the Potential and Current Limitations of Integrating CFD and Structural Analyses for the Design of Offshore Wind Turbines Julius Kandt Department of Mechanics and Maritime Sciences Chalmers University of Technology Abstract While the need for electric energy grows over the years, the environmental impact is an important factor for the energy generation as well. Wind power plays a major role among the renewable energies nowadays. However, the available space onshore is restricted and complex terrain complicates the use of various areas even more. Among others this is one of the reasons why the use of offshore wind power stations becomes more attractive. Therefore, the assessment of offshore sites is an important aspect that may decide whether to place a wind plant or not. For this Computational Fluid Dynamics have become a popular approach over the past years replacing simpler spectral methods. In this study the modeling of the Atmospheric Boundary Layer for the generation of an offshore environment using a Large Eddy Simulation was investigated. Overall it was found, that this method is able to reproduce the Atmospheric Boundary Layer well, especially when it comes to the mean velocity field. Nevertheless, inaccuracies were found considering the Reynolds stresses, especially the shear components. Moreover, the flow field obtained with the Large Eddy Simulation was compared to a spectral flow field and the dynamic response of two wind turbines was assessed using the aero-elastic solver FAST. The results of both methods are in close agreement to each other. Therefore, the integration of the aero-elastic solver FAST into the process of an analysis using Computational Fluid Dynamics works well. Finally, the performance of the Actuator Disk and Actuator Line model was tested in Star-CCM+. The effect of these models on the flow field was investigated as well as the thrust and torque calculation. It was found that both models do not match the FAST results very well in terms of the torque prediction. However, the thrust calculation of the Actuator Line Model comes close to the results from FAST. Several aspects were outlined for further investigation. Keywords: Wind turbine, Offshore, Atmospheric Boundary Layer, Large Eddy Sim- ulation, Actuator Disk Model, Actuator Line Model, FAST v Acknowledgements I would like to thank all the people that made this thesis possible for all their in- put and suggestions for improvements. Most of all I want to thank my supervisor Hamidreza Abedi. He has been very supportive throughout the entire project, pro- viding valuable input and inspiring with his ideas. Furthermore, I would like to thank Lars Davidson and the Division of Fluid Dynamics, who made this project possible in the first place. Finally, I would like to thank my family and friends, who supported and encouraged me the entire time. Julius Kandt, Gothenburg, June 2019 vii Contents List of Figures xiii List of Tables xv Nomenclature xviii 1 Introduction 1 1.1 Historical Background . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Simulation Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Objectives and limitations . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Theory 5 2.1 Atmospheric Boundary Layer . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1 Onshore Atmospheric Stability . . . . . . . . . . . . . . . . . 6 2.1.2 Marine Atmospheric Boundary Layer . . . . . . . . . . . . . . 8 2.2 Computational Fluid Dynamics (CFD) . . . . . . . . . . . . . . . . . 8 2.2.1 Historical Background . . . . . . . . . . . . . . . . . . . . . . 9 2.2.2 Navier-Stokes (NS) Equations . . . . . . . . . . . . . . . . . . 9 2.2.3 Reynolds-Averaged Navier-Stokes (RANS) Equations . . . . . 10 2.2.4 Spatial Filtered Navier-Stokes Equations . . . . . . . . . . . . 11 2.2.5 Turbulence Modeling . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.6 Modeling of the Atmospheric Boundary Layer . . . . . . . . . 20 2.2.7 Turbulent Quantities . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Wind Turbine Aerodynamics and Loads . . . . . . . . . . . . . . . . 22 2.3.1 Aerodynamic Characteristics . . . . . . . . . . . . . . . . . . . 22 2.3.2 Actuator Disk Momentum Theory . . . . . . . . . . . . . . . . 23 2.3.3 Blade Element Method and Actuator Line Model . . . . . . . 27 2.3.4 Blade-Element-Momentum (BEM) Theory . . . . . . . . . . . 28 2.4 Reference Turbine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3 Methods 33 3.1 Computational Domain . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.1.1 Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2 RANS Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3 LES Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Actuator Disk and Actuator Line Model . . . . . . . . . . . . . . . . 36 ix Contents 3.5 Aero-elastic Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4 Results 39 4.1 Atmospheric Boundary Layer Modeling . . . . . . . . . . . . . . . . . 39 4.1.1 Time-Averaged Flow Field . . . . . . . . . . . . . . . . . . . . 39 4.1.2 Instantaneous Flow Field . . . . . . . . . . . . . . . . . . . . . 44 4.2 Wind Turbine Response . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.3 Actuator Disk and Actuator Line Model . . . . . . . . . . . . . . . . 49 5 Conclusion 57 Bibliography 59 A Appendix 1 I x List of Figures 2.1 Thickness of the Atmospheric Boundary Layer. Figure taken from [50] 5 2.2 Vertical structure of the Atmospheric Boundary Layer. Figure taken from [13] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3 Typical diurnal cycle of the Atmospheric Boundary Layer during fair weather. Figure taken from [50] . . . . . . . . . . . . . . . . . . . . . 6 2.4 Variation of wind speed due to atmospheric stability. Figure adapted from [50] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.5 Structure of the Marine Atmospheric Boundary Layer. Figure taken from [13] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.6 The energy density spectrum of turbulent kinetic energy. Figure taken from [11] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.7 Overview over commonly used turbulence models. Figure adapted from [39] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.8 Illustration of the Actuator Disk model. Figure taken from [8] . . . . 23 2.9 Axial velocity and static pressure curve of the flow over an Actuator Disk. Figure taken from [8] . . . . . . . . . . . . . . . . . . . . . . . 24 2.10 Cell marking within the blade element method. Figure taken from [1] 27 2.11 Steady-state behaviour of the NREL offshore 5-MW baseline wind turbine. Figure taken from [23] . . . . . . . . . . . . . . . . . . . . . 31 3.1 Computational Domain used for the CFD simulation in Star-CCM+ . 34 3.2 Sideview of the mesh for the Actuator Disk Model and Actuator Line Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.1 Comparison of the mean streamwise velocity component between the logarithmic law according to Richards & Hoxey, the RANS simulation and the LES at the location of turbine 2 . . . . . . . . . . . . . . . . 40 4.2 Comparison of the turbulent kinetic energy profiles between the RANS simulation and LES at the location of turbine 2 . . . . . . . . . . . . 40 4.3 Comparison of the turbulence intensity profile between the RANS simulation and LES at the location of turbine 2 . . . . . . . . . . . . 41 4.4 Turbulent length scale of the RANS simulation at the location of turbine 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.5 Normal stresses of the LES at the location of turbine 2 . . . . . . . . 43 4.6 Shear stresses of the LES at the location of turbine 2 . . . . . . . . . 43 4.7 Instantaneous velocity components at hub height at turbine 2 . . . . 44 xi List of Figures 4.8 Energy density spectrum of the instantaneous spanwise velocity com- ponent at hub height at turbine 2 . . . . . . . . . . . . . . . . . . . . 45 4.9 Spatial correlation of the velocity components in the z direction at hub height at turbine 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.10 Comparison of the instantaneous streamwise velocity extracted from the FAST simulation at hub height at turbine 2 of the spectral-based versus the CFD-based flow field . . . . . . . . . . . . . . . . . . . . . 46 4.11 Comparison of the rotor speed of turbine 2 of the spectral-based ver- sus the CFD-based FAST simulation . . . . . . . . . . . . . . . . . . 47 4.12 Comparison of the generated power of turbine 2 of the spectral-based versus the CFD-based FAST simulation . . . . . . . . . . . . . . . . 48 4.13 Blade deflection of turbine 2 of the spectral-based versus the CFD- based FAST simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.14 Comparison of the streamwise velocity component between the Actu- ator Disk and Actuator Line simulation near turbine 1 . . . . . . . . 50 4.15 Comparison of the turbulent kinetic energy between the Actuator Disk and Actuator Line simulation near turbine 2 . . . . . . . . . . . 51 4.16 Comparison of the turbulence intensity between the Actuator Disk and Actuator Line simulation near turbine 2 . . . . . . . . . . . . . . 52 4.17 Comparison of the normal component Rexx of the Reynolds stresses between the Actuator Disk and Actuator Line simulation near turbine 2 53 4.18 Comparison of the shear component Rexy of the Reynolds stresses between the Actuator Disk and Actuator Line simulation near turbine 2 53 4.19 Comparison of the thrust between the Actuator Disk and Actuator Line simulation at turbine 2 . . . . . . . . . . . . . . . . . . . . . . . 54 4.20 Comparison of the torque between the Actuator Disk and Actuator Line simulation at turbine 2 . . . . . . . . . . . . . . . . . . . . . . . 55 A.1 Tower-top and tower-base coordinate systems for the FAST simula- tions in the current study. The tower-base coordinate system is fixed in the support platform. It translates and rotates with the platform. The tower-top coordinate system is fixed to the top of the tower. It translates and rotates as the platform moves, but does not yaw with the nacelle [25] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I A.2 Tower-top roll moment of turbine 2. It is aligned to the x-axis of the tower-top coordinate system. . . . . . . . . . . . . . . . . . . . . . . . II A.3 Tower-top pitch moment of turbine 2. It is aligned to the y-axis of the tower-top coordinate system . . . . . . . . . . . . . . . . . . . . . II A.4 Tower-top yaw moment of turbine 2. It is aligned to the z-axis of the tower-top coordinate system . . . . . . . . . . . . . . . . . . . . . . . III A.5 Tower-based roll moment of turbine 2. It is aligned to the x-axis of the tower-base coordinate system . . . . . . . . . . . . . . . . . . . . III A.6 Tower-based pitch moment of turbine 2. It is aligned to the y-axis of the tower-base coordinate system . . . . . . . . . . . . . . . . . . . . IV A.7 Tower-based yaw moment of turbine 2. It is aligned to the z-axis of the tower-base coordinate system . . . . . . . . . . . . . . . . . . . . IV xii List of Figures A.8 Comparison of the streamwise velocity component between the Actu- ator Disk and Actuator Line simulation near turbine 2 . . . . . . . . V A.9 Comparison of the turbulent kinetic energy between the Actuator Disk and Actuator Line simulation near turbine 1 . . . . . . . . . . . V A.10 Comparison of the turbulence intensity between the Actuator Disk and Actuator Line simulation near turbine 1 . . . . . . . . . . . . . . VI A.11 Comparison of the streamwise normal component of the Reynolds stresses Rexx between the Actuator Disk and Actuator Line simula- tion near turbine 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . VI A.12 Comparison of the normal component of the Reynolds stresses in the wall-normal direction Reyy between the Actuator Disk and Actuator Line simulation near turbine 2 . . . . . . . . . . . . . . . . . . . . . . VII A.13 Comparison of the normal component of the Reynolds stresses in the spanwise direction Rezz between the Actuator Disk and Actuator Line simulation near turbine 2 . . . . . . . . . . . . . . . . . . . . . . . . . VII A.14 Comparison of the shear component of the Reynolds stresses Rexz between the Actuator Disk and Actuator Line simulation near turbine 2VIII A.15 Comparison of the shear component of the Reynolds stresses Reyz between the Actuator Disk and Actuator Line simulation near turbine 2VIII xiii List of Figures xiv List of Tables 2.1 Typical model constants for the k-ε model [10] . . . . . . . . . . . . . 18 2.2 Characteristics of the NREL offshore 5-MW baseline wind turbine [23] 30 3.1 Initial conditions, reference values and fluid characteristics for the RANS simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.1 Turbulent length scale calculated from the spatial two-point correla- tion of the velocity fluctuations in z direction . . . . . . . . . . . . . 46 4.2 Output from the FAST simulations with a flow field from spectral analysis as well as from CFD simulations . . . . . . . . . . . . . . . . 49 4.3 Output from the Actuator Disk Model, Actuator Line Model and FAST simulations for turbine 2 . . . . . . . . . . . . . . . . . . . . . 55 A.1 Different momentums from the FAST simulations with a flow field from spectral analysis as well as from CFD simulations for the tower- top coordinate system at turbine 2 . . . . . . . . . . . . . . . . . . . I A.2 Different momentums from the FAST simulations with a flow field from spectral analysis as well as from CFD simulations for the tower- based coordinate system at turbine 2 . . . . . . . . . . . . . . . . . . II xv List of Tables xvi Nomenclature Abbreviations ABL Atmospheric boundary layer. ADM Actuator Disk Model. ALM Actuator Line Model. AOA Angle Of Attack. BEM Blade-Element-Momentum. CFD Computational Fluid Dy- namics. DES Detached Eddy Simulation. DNS Direct Numerical Simula- tion. DSM Dynamic Smagorinsky Model. FAST Fatigue, Aerodynamics, Structures and Turbulence. LES Large Eddy Simulation. MABL Marine Atmospheric Boundary Layer. NREL National Renewable Energy Laboratory. NS Navier-Stokes. RANS Reynolds-Averaged Navier-Stokes. SGS Subgrid-Scales. SST Shear Stress Transport. TSR Tip-Speed Ratio. URANS Unsteady Reynolds- Averaged Navier-Stokes. WAsP Wind Atlas Analysis and Application Program. Subscripts 0 Free stream. 1 Far wake. i, ij Einstein summation notation. k Turbulent kinetic energy. t Turbulent quantity. in Inner disk. norm Normal component. out Outer disk. ref Reference value. rel Relative with respect to a blade. tang Tangential component. Greek Letters αt Thermal diffusion [m2 s ] β Inflow angle [◦] ∆ Filter width [m] δ Actuator disk thickness [m] δij Kroenecker’s δ [−] ε Turbulent dissipation [m2 s3 ] η Kolmogorov microscale [m] κ Wavenumber [m−1] xvii List of Tables κw von Kármán constant [−] Λ Turbulent length scale [m] λ Tip-speed ratio [−] µ Dynamic viscosity [ kg ms ] ν Kinematic viscosity [m2 s ] νe Subgrid-scale eddy viscosity coefficient [m2 s ] νt Turbulent viscosity [m2 s ] Ω Rotational rotor speed [ rad s ] ω Induction angular velocity [ rad s ] ρ Density [ kg m3 ] σ Solidity [−] σk Turbulent Prandtl number [−] σε Model constant for the k-ε model [−] τij Reynolds stress tensor [ kg ms2 ] τ̃ij Subgrid-scale stress tensor [m2 s2 ] Roman Letters a Axial induction factor [−] A Section area [m2] b Tangential induction fac- tor [−] c Chord length [m] cε1, cε2, cµ Model constants for the k-ε model [−] C Wall function constant [−] Cd Drag coefficient [−] Cl Lift coefficient [−] Cs Smagorinsky constant [−] Ct Thrust coefficient [−] CK Kolmogorov constant [−] D Drag force [kgm s2 ] D̃ij grid-scale rate-of-strain tensor [−] E Energy density spectrum [m2 s2 /Hz] EC Wall function constant [−] F Force [kgm s2 ] I Turbulence intensity [−] k Turbulent kinetic Energy [m2 s2 ] l Characteristic length [m] lm Mixing length [m] L Lift force [kgm s2 ] m Mass element [kg] p Hydrodynamic pressure [ kg ms2 ] Q Torque [kgm2 s2 ] r Radius [m] rh Roughness height [m] R Correlation coefficient [−] Re Reynolds number [−] Reii Normal components of the Reynolds stresses [m2 s2 ] Reij Shear components of the Reynolds stresses [m2 s2 ] t Time [s] T Thrust force [kgm s2 ] u Characteristic velocity in a steady state [m s ] u′ Fluctuation velocity [m s ] u∗ Friction velocity [m s ] vi Velocity component [m s ] V Axial velocity [m s ] xi Cartesian coordinate [−] y0 Roughness length [m] xviii 1 Introduction Due to a significant growth in the worlds population the need for electric energy increased rapidly over the last decades. While fossil fuels such as coal, oil and natural gas built the foundation of the energy supply since the beginning of the industrial evolution, there have been various concerns about the impact of these resources on the environment in recent years [18]. Global challenges like the climate change, the emission of greenhouse gases and air pollution require a change in the worldwide energy production [35]. Although nuclear power does produce only low carbon dioxide emissions, the question of what to do with radioactive waste is still not solved and a topic of huge controversies [12]. Hence, the development of sustainable and clean energies is a matter of interest in many regions all over the world. Among multiple forms of renewable energies, wind power plays an important role and is considered the most sustainable renewable energy technique [16]. Furthermore humans have used wind power for a long period of time and therefore a lot of expertise has been gathered. 1.1 Historical Background While the use of wind power is almost 3000 years old, its purpose as generation of electric energy started about 120 years ago. Since then wind power had a continuous development in terms of its technology. Although there always has been a dependency on oil prices, the electric capacity of wind power has been growing about 30 % per year over the last decade [2]. In fact there are predictions, that wind energy will generate about 5 % of the worlds energy demand in 2020 [19]. As one of its inventors Germany is one of the leading countries in terms electric energy production using wind power. The fraction of renewable energies was about 38.2 % net in 2017 and wind power provided about 18.8 % of the total net electric energy production [6][15]. In terms of the use of renewable energies Sweden is even more advanced. In this matter it leads the European Union and reached its target of producing half of its energy from renewable sources in 2012. In 2018 about 54 % came from renewable energy sources and about 11 % from wind power [46]. In order to generate more electricity, wind power stations are growing to larger scales. In the period from 2010 to 2017 the average power produced by a wind tur- bine increased from 1.77 MW to 2.7 MW [14]. However the space available onshore is limited or may have a difficult topography. 1 1. Introduction This is one of the reasons offshore wind is an attractive alternative. Moreover there are several other advantages. The sites available offshore may not only be used to increase the number of wind turbines but also to enhance their size. Furthermore the impact on human life is much less. Hence the focus on development may be more efficiency orientated without the concern of noise emission or visual distur- bances [5]. Another key fact is that offshore winds often blow stronger and more uniform than onshore. This may reduce the stresses on the turbine induced by the wind. Additionally the energy production indicator for offshore wind power stations is around 4000 full load hours per year, while it is around 2000-2500 for onshore sites [30]. Nevertheless there is still development necessary in order to reduce the high costs for offshore wind power stations, which are around 50 % higher than onshore [30]. To improve the performance and expand the lifetime of wind power stations numer- ous experiments have been conducted over the last decades. In addition compu- tational calculation and prediction methods have also improved among others due to increasing computing power. Multiple programs were invented, one of the most commonly used is Wind Atlas Analysis and Application Program (WAsP) [34], which is based on the wind-atlas model [47]. This model solves linearised fluid equations, which cannot resolve detached flow or recirculation. While it performs well on flat terrain, it tend to have problems with complex terrain or when analyzing bigger wind farms [7] which is one of the main interest in today’s research. However, since computational power became more affordable Computational Fluid Dynamics (CFD) has been applied to wind energy research to refine the results and increase accuracy. Since CFD is a non-linear method it makes it possible to capture wake effects. This is especially important for the efficiency of large-scale wind farms, as one can assess the interactions between the turbines with each other. Furthermore, it is important for the structural aspects as well. Therefore, it can affect the lifetime as well, which is a crucial point in terms of the costs. 1.2 Simulation Methods Overall, there are various types of simulations available today for assessing wind fields, the response of wind turbines as well as predicting the power outcome. Be- sides linear applications like WAsP nonlinear CFD methods are becoming more common, especially for evaluating the effect of complex terrain and wake assess- ment on single turbines and wind parks. Currently two different CFD methods are in use. On the one hand, there are the time averaged models, e.g. Reynolds-Averaged Navier-Stokes (RANS) models and on the other the time resolved approaches, e.g. Large Eddy Simulation (LES). Fur- thermore, so-called Hybrid methods are becoming more popular. They combine different methods, e.g. RANS and LES. The full resolving Direct Numerical Simu- lation (DNS) is still too power demanding to have actual use in wind engineering. This is mainly caused by the large variety of scales in the Atmospheric Boundary 2 1. Introduction Layer (ABL) and will be discussed later on. However, since atmospheric flow is a turbulent and therefore time depending phenomenon, LES is a promising approach, especially with regard to the development of the computer power. Furthermore, CFD methods have been developed to assess the influence of one or more wind turbines on the flow field without actually resolving the turbine geome- try. Two of these methods used in the present study are the so-called Actuator Disk Model (ADM) and Actuator Line Model (ALM). It is possible to study the flow field of single turbines or multiple turbine up to whole wind parks. Besides the flow assessment another important aspect is the impact of the flow on the wind turbine. Determining the response of the wind turbine to the loads and stresses induced by the wind is a crucial aspect for the design of the structure and lifetime. A popular approach for the combined use with CFD software is Blade- Element-Momentum (BEM) theory, which is a combination of two different methods itself (blade element and momentum theory). 1.3 Literature Review The modeling of the ABL has been the objective of various studies. One of the most popular approaches for modeling the mean flow with the k-ε model was proposed by Richards and Hoxey [38]. They derived the inlet conditions analytically so that they fulfill the governing equations of the k-ε model. This produces a homogeneous flow field throughout the domain. Hargreaves and Wright [17] improved this model by applying shear stresses at the top of the domain. This leads to better results, es- pecially for the turbulent kinetic energy. Other modifications were made by Parente and Benocci [36], who introduced additional source terms to the transport equations of the turbulent kinetic energy and the turbulent dissipation. For transient modeling of the ABL Porté-Agel et. al. [37] published a well-known article assessing different subgrid-scale models for LES. Comparing the introduced scale-dependent dynamic model with the standard Smagorinsky model and the orig- inal dynamic model he found improvements in terms of dissipative properties, espe- cially in near wall regions. Furthermore, Vasaturo et. al. [48] investigated different inflow methods for LES. They conclude, that the most accurate method is the pre- cursor method, which means running a precursor simulation to obtain the inflow quantities. To model especially the wake behind turbines and wind farms Stevens et. al. [44] studied the performance of the ADM and ALM and compared the results with wind tunnels experiments. In contrast to earlier results they concluded that both the ADM and ALM are capable of reproducing accurate results, especially if the na- celle and tower effects are included. Furthermore, the ALM is assumed to have its advantages over the ADM, that is the resolution of tip vortices and radial distribu- tion alongside the blades, only in high resolution grids. This is also validated by Churchfield et. al. [9] and other researchers [40][53]. Moreover, Martínez et. al. 3 1. Introduction investigated the influence of the mesh resolution on the predicted power output of both models and conclude that the power decreases with the number of cells [31]. 1.4 Objectives and limitations The main objective of this thesis is to study the potential and current limitations of integrating CFD methods and structural analysis for the design of offshore wind turbines. This main objective can be structured in three different aims. The first one is the generation of a transient flow field through LES for further investigations. To obtain a transient flow field with good quality, one needs to create a mean velocity field and describe the turbulence characteristics as input parameters for the LES. In the present study this is done with a RANS simulation. Once the transient flow field is obtained it can be transmitted to FAST for aero-elastic assessment. Then it will be compared to the original spectral input. Finally, integrated CFD methods like the ADM and ALM will be investigated and compared to the FAST simulation. The following list gives an overview over the detailed objectives of this thesis. • Set up of a RANS simulation to simulate the ABL with good accuracy and agreement with the current theoretical models and literature. • Perform a LES to produce a transient flow field. Compare this flow field to previous RANS simulation and literature. • Transmit the flow field for two different turbine locations to FAST and per- form aero-elastic simulations for each turbine. Compare the turbines to each other and the original spectral flow field as input. • Perform a LES with the ADM and ALM applied to it, respectively. Assess the influence of the respective models on the flow field. Compare the turbines and investigate the influence of the first turbine on the second. Calculate force and momentum outputs and compare them with the FAST simulation. Nevertheless, there are some limitations to the present study. The limitations for the CFD calculations are mainly the computational resources as well as the models and possibilities available in the CFD software Star-CCM+ and the aero-elastic solver FAST. Further on the project is limited by the time available to prepare and assess the simulations. 4 2 Theory In this chapter the theoretical foundation of the thesis is discussed. Therefore, the Atmospheric Boundary Layer (ABL) and its different conditions are described, as it is one of the main parts of this investigation. Fundamental basics like structure and fluid motions within the ABL as well as differences between on- and offshore wind plants are outlined. Moreover, the concept of Computational Fluid Dynamics (CFD) is presented including different models used in this thesis. The governing equations are derived and the mathematical description of turbulence and how to deal with it is shown. Furthermore the basic aerodynamics of wind turbine rotors are presented as well as the fundamental theory for aerodynamic analysis methods with the BEM approach are outlined. Finally the reference turbine used in this thesis is described. 2.1 Atmospheric Boundary Layer The ABL surrounds the earth’s surface and extends to 10-20 % of the troposphere. It is defined as the part of the troposphere that is influenced by the earth’s surface [27][50]. Its height differs dramatically from tens of meters up to 4 km [50], depend- ing on various boundary conditions, see e.g. sub-section 2.1.1. Furthermore, the thickness of the ABL varies with time and space, as shown in figure 2.1. Figure 2.1: Thickness of the Atmospheric Boundary Layer. Figure taken from [50] The ABL can be divided into three different regions, as figure 2.2 shows. One can see the roughness sublayer next to the ground, followed by the constant-flux or Prandtl sublayer and finally the Ekman sublayer. Within the roughness sublayer one can also find the laminar or viscous sublayer adjacent to the surface. It measures only a few millimeters and is therefore too thin to be shown in figure 2.2. The Prandtl sublayer is meteorologically defined as the layer, where the turbulent vertical momentum 5 2. Theory flux, heat and moisture deviate less than 10 % from their surface values and the influence of the Coriolis force is negligible. Usually it covers about 10 % of the total ABL thickness. However, modern wind turbines tend to exceed the Prandtl layer and reach into the Ekman layer. One of the main differences between these two sublayers is that the Coriolis force can not be neglected [13]. How to mathematically describe the velocities in the ABL and within the different layers is discussed in sub-section 2.2.6. Figure 2.2: Vertical structure of the Atmospheric Boundary Layer. Figure taken from [13] 2.1.1 Onshore Atmospheric Stability The atmospheric stability status has great influence on the extend of the ABL. Usually the stability changes depending on the temperature difference between the earth’s surface and the surrounding fluid [50]. A typical cycle during fair weather (which means no precipitation or clouds and moderate wind speeds) conditions over land can be seen in figure 2.3. Figure 2.3: Typical diurnal cycle of the Atmospheric Boundary Layer during fair weather. Figure taken from [50] 6 2. Theory After the sun rises, the surface heats up and the mixed layer grows due to heat con- vection. Therefore, it is also called the convective boundary layer. Strong vertical mixing takes place in this layer, thus the vertical gradients are small. Moreover, en- trainment occurs, because the mixed layer is expanding, hence inducing fluid motion in the surrounding air. This entrainment region also builds the capping inversion, which is the top end of the ABL. When the sun sets and the ground cools down and the stable boundary layer develops. Here the vertical mixing is not very distinct and therefore vertical gradients are higher than in the mixed layer. Above the stable boundary layer the residual layer occurs, as it is a remnant from the mixed layer during the daytime [50]. In total there are three different stability states for the ABL. The state is categorized as unstable if the surface is warmer than its surroundings and therefore heating up the air. This is because of the corresponding change in the air’s density and result- ing fluid motions. The air closer to the ground heats up, the density sinks and as a consequence the warmer air ascents to greater heights. During its way up the air cools down and sinks to the ground again. As a result the vertical convection in the ABL is very high, turbulence is produced and the boundary layer is well mixed. Also the ABL total thickness increases. This is typically the case during daytime. On the contrary the ground is often colder than the air during nighttime and this state is defined as a stable boundary layer, because no fluid motion is induced by cooling the air adjacent to the ground. Finally the ABL is called neutral, if the potential temperature is constant with height. Typically this is only a short period of time at late afternoon [27]. The stability has great influence on the velocity profile of the wind, as it can be seen in figure 2.4. Because of the highest vertical convection the unstable ABL has also the highest velocity gradient near the ground. After only a short vertical distance the velocity profile is fairly uniform. On the opposite side the stable ABL has the lowest velocity gradient and the neutral ABL lies somewhere in between [50]. Figure 2.4: Variation of wind speed due to atmospheric stability. Figure adapted from [50] 7 2. Theory 2.1.2 Marine Atmospheric Boundary Layer The Marine Atmospheric Boundary Layer (MABL) is somewhat different from the onshore ABL. One of the main differences is the roughness length, which is much less for the offshore, since the surface is much smoother than onshore flat terrain. As a consequence offshore wind speeds at a given height are much higher, the tur- bulence intensity is less and the thickness of the surface layer is smaller as well. Less turbulence often leads to less shear over the wind turbine. In contrast to the on- shore characteristics, the offshore roughness length is dependent on the wind speed because of the wave influence (it increases with the wind speed). Also the diurnal cycle is nearly absent due to large thermal storage capacities of the sea [13]. Figure 2.5: Structure of the Marine Atmospheric Boundary Layer. Figure taken from [13] The vertical structure of the MABL is shown in figure 2.5. Adjacent to the surface there is the wave sublayer, which is directly influenced by single waves due to the dominant pressure forces. The wave sublayer is usually about five wave amplitudes deep. Above that one can find the constant-flux sublayer or Prandtl layer. This layer is, as the entire MABL, much shallower than the ABL over land. The Ekman sublayer represents the upper 90 % of the MABL and then turns into geostrophic wind. Note, that the change in wind speed at hub height in the present study occurs within the Ekman layer and therefore the same law for the mean wind speed as for onshore wind turbines can be used [13]. This is discussed later on in section 2.2.6. 2.2 Computational Fluid Dynamics (CFD) CFD is the computer-based simulation and assessment of systems with fluid flow, heat transfer and other associated phenomena e.g. chemical reactions [49]. It is a powerful technique that can be used in a large variety of applications and industries. 8 2. Theory 2.2.1 Historical Background The industrial beginnings of CFD reach back into the 1960s, where the aerospace industry first integrated CFD into the design of aircraft and jet engines. With the further development of CFD and more computer power available over the years, it was also used to design internal combustion engines, combustion chambers for gas turbines and furnaces. Nowadays even car manufacturers use CFD on a daily basis to predict essential physical characteristics, e.g. drag forces and even the interior car environment [49]. However, due to its complex problems and relations the capa- bility of CFD is still a bit behind other computer aided engineering tools, e.g. stress analysis tools. Despite the development, that still has to be made, there are various arguments in favor of using CFD summarized in the following list [42]: • When trying to find an analytical solution for fluid mechanical problems, only highly simplified situations can be resolved. In some cases these simplifications do not represent the physical situation, leaving the initial problem unsolved. • Although empirical relations might exist and be documented for some prob- lems one can only adapt these relations for very similar situations. • An experimental assessment is not possible in many flows due to different circumstances. The assessed fluid might be very hot or chemically aggressive and therefore the measurements devices might get destroyed while in contact with the fluid. When using invasive methods, the sensors might distort the results. On the opposite site non-contact method might not be applicable for the assessed flow, e.g. due to missing visibility and other reasons. • In terms of costs CFD is getting more and more affordable due to the lowering cost for computer power, but the costs for experimental devices do not profit from such a phenomenon. • It is often easier to conduct parameter studies in CFD than in experimental methods. However, often both CFD and experimental analysis is used for obtaining best pos- sible knowledge about a flow field. On the one hand experiments validate CFD analysis in certain situations, on the other hand, CFD can help to understand the flow field. 2.2.2 Navier-Stokes (NS) Equations In the following sections the governing equations for the present study is intro- duced. For this, one presents the Navier-Stokes (NS) equations and then derive the Reynolds-Averaged Navier-Stokes (RANS) equations as well as the spatial filtered NS equations, the so-called Large Eddy Simulation (LES). The NS equations describe four equations in total. One equation represents the con- servation of mass, also known as the continuity equation. The three other equations describe the conservation of momentum in each direction. They are derived from Newton’s second law applied to an element of fluid and a given mass. 9 2. Theory The continuity equation for an incompressible fluid in Cartesian coordinates can be written as ρ ∂vi ∂xi = 0, (2.1) where ρ denotes the density, which is constant is this case [11]. The momentum equations with the Einstein summation notation1 can be formulated as ρ ∂vi ∂t + ρ ∂vivj ∂xj = − ∂p ∂xi + µ ∂2vi ∂xj∂xj . (2.2) xi is used for the Cartesian coordinate system as well, where i = 1, 2 or 3 for x, y or z, while vi denotes the velocity field with i = 1, 2 or 3 for u, v or w. The hydro- dynamic pressure is represented by p, the dynamic viscosity by µ [11]. Note that µ is constant in the present case and the external forces acting on the fluid element have been neglected in this equation. The NS equations (equations 2.1 and 2.2) build the foundation of fluid dynamics as they fully describe the flow, including all time and length scales. This leads to problems in some cases, where the range of scales is very large, e.g. in ABL simulations. The largest turbulent scales are in the order of the geometric boundary conditions and therefore around 103 m, while the smallest scales are about the size of the dissipative eddies, which is about 10−3 m [45]. To resolve all of these scales would require a massive amount of computer power, and thus it is infeasible. Moreover, in many cases in CFD it is not needed to resolve all scales, as only a part of them are in the field of interest. Sometimes one only needs to have information about the mean field and therefore does not need to resolve transient behavior at all. Thus, several approaches have been developed to make this kind of simulation feasible. The RANS simulation apply the Reynolds decomposition and time averaging to the NS equations and consequently reduce the computer time. LES apply spatial filters to the NS equations and therefore only resolve a part of the turbulent structures. Hence, they do provide transient solutions, but only calculate a certain spectrum of all turbulent scales. 2.2.3 Reynolds-Averaged Navier-Stokes (RANS) Equations For the derivation of the RANS equations one needs to perform a Reynolds decom- position. This splits the variables into a time-averaged component and a fluctuating component. After that the equation needs to be time-averaged to get to the fi- nal RANS formulation. One wants to point out, that although the equations is time-averaged, the transient term of the NS-equations can be retained. This ap- proach is called unsteady RANS (URANS) and is capable of resolving behavior on larger timescales than the averaging time [11]. However, since in this thesis only 1The Einstein summation notation is a convention, which implies summation over an indexed term in a formula. If an index appears twice in a single term it implies the summation over all values of that index. 10 2. Theory the time-averaged RANS equations are used, the derivation of the RANS equations without retaining the transient term is done. A variable v is decomposed into the time-averaged value, denoted by v, and the fluctuation value, denoted by v′. The decomposition now can be written as v = v+ v′. (2.3) By applying this decomposition to the NS equations and conducting time-averaging, the final incompressible RANS equations are obtained. Additionally a steady mean flow is assumed [11]. The continuity equation is written as ρ ∂vi ∂xi = 0. (2.4) The time-averaged momentum equations reads ρ ∂vi vj ∂xj = − ∂p ∂xi + ∂ ∂xj ( µ ∂vi ∂xj − ρv′iv′j ) (2.5) An important difference between the time-dependent NS equations and the (time- averaged) RANS equations is the new term, the so called Reynolds stress tensor τij = ρv′iv ′ j. This term is an important part of turbulence modeling and therefore a few characteristics are explained here. The tensor consists of the correlations between the different fluctuating velocities. The diagonal components represent the normal stresses, whereas the off-diagonal quantities represent the shear stresses. The basic structure of the tensor is shown in equation 2.6. It reads τij = ρ v′iv ′ j = ρ  u′2 u′v′ u′w′ v′u′ v′2 v′w′ w′u′ w′v′ w′2  . (2.6) One can see, that the tensor is symmetric and therefore u′v′ = v′u′, u′w′ = w′u′ and v′w′ = w′v′. In general the tensor represents an additional viscosity by introducing an additional stress term due to turbulent fluid motions. Additionally the Reynolds stress tensor is an unknown quantity. It introduces six new variables to the RANS equations. For the whole RANS equations there are now ten unknown variables, consisting of three velocity components, the pressure and the six Reynolds stresses. This is called a closure problem, because there are more unknown variables than number of equations. Thus several models have been developed to close the RANS equations and the main difference between these models is the way they resolve the Reynolds stresses [11][39]. 2.2.4 Spatial Filtered Navier-Stokes Equations Another approach used in the present study is the Large Eddy Simulation (LES). The main idea is to resolve the large energy-carrying eddies and to model the smaller scales. These eddies are expected to behave more homogeneous and isotropic accord- ing to the Kolmogorov theory, which is explained later in this section. Therefore, modeling is easier and closer to reality. Furthermore, the smaller the eddies the less 11 2. Theory they depend on the boundary conditions of the domain. Hence, it is easier to create universal valid models, which can be used in different applications. Compared to Direct Numerical Simulation, known as the most expensive CFD method, in LES, the mesh can be coarsened, because only larger eddies need to be captured. This makes the LES approach significantly less resource demanding than DNS [26][49]. In order to receive the governing equations for LES, one needs to filter the physi- cal quantities in the NS equations. Therefore, LES resolve transient behavior and these filtered variables are dependent on time and space. Similar to the Reynolds decomposition, a variable v is then written as v = ṽ + v′′, (2.7) where ṽ denotes a spatial filtered variable and v′′ denotes the subgrid-scale (SGS) component. Applying the filtering to the NS equations yields to the incompressible spatial filtered equations used for LES [26]. These equations have the same form as the initial NS equations. The spatial filtered continuity equation reads ρ ∂ṽi ∂xi = 0, (2.8) while spatial filtered the momentum equation reads ∂ṽi ∂t + ∂(ṽiṽj) ∂xj = −1 ρ ∂p̃ ∂xi + ∂ ∂xj (−τ̃ij + 2νD̃ij). (2.9) ν denotes the kinematic viscosity [26]. Note that D̃ij is the so called grid-scale rate-of-strain tensor, which is defined as D̃ij = 1 2 ( ∂ṽi ∂xj + ∂ṽj ∂xi ) . (2.10) The SGS stresses are given by τ̃ij = ṽivj − ṽiṽj. (2.11) Note that τ̃ij is used, to distinguish between the Reynolds stress tensor. To be precise it is incorrect to denote τ̃ij as a stress, since it is not multiplied by ρ and therefore has the wrong dimension [11]. However, it is often referred as stress and one will also refer to it as a stress in this study. An important part of every LES is the type of filter that is used. The filter defines which range of scales is resolved and which is modeled. Therefore, has a direct impact on the need of computer resources. A very common filter is the so called box filter, which is always used in finite vol- ume method [11]. With the box filter the specific variable is averaged over the entire interval, or, in the case of a volume, over all dimensions. In the present study the filter size equals the computational grid. Thus, all scales that are captured depend on the grid size. Moreover, one can denote the scales, which are not resolved with a certain grid as subgrid-scales. The SGS stresses then represent the impact of the 12 2. Theory unresolved small scale structures on the larger resolved structures. As introduced in equation 2.11, these stresses are unknown and need to be modeled, very similar to the previous presented RANS approach [11][26]. How much of all the scales should be captured can be derived from the energy density spectrum, which is shown in figure 2.6. While turbulence in general contains a wide range of different time and length scales, certain structures can be distinguished. Looking at figure 2.6, one can see three different parts. First, there are the large, energy-carrying scales, denoted by I. Second, there is the inertial subrange, labeled with II. This is also called Kolmogorov −5/3-range, which is discussed down below. Finally, the third part is the dissipation subrange, denoted by III [11]. Figure 2.6: The energy density spectrum of turbulent kinetic energy. Figure taken from [11] The different scales and the interactions between them are described by the so called energy cascade. In section I of the spectrum, the large eddies (and turbulence itself) are generated. This is denoted by the production term − 〈 v′iv ′ j 〉 ∂〈vi〉 ∂xj above the ordinate, which is discussed in section 2.2.5. The direction of motion is most likely dictated by outer conditions, i.e. the mean flow direction. Therefore, these eddies are highly anisotropic. These large eddies then transfer their energy to smaller eddies due to vortex stretching and velocity gradients. The eddies get smaller and at the same time less depending on outer parameters. Hence, at some point the eddies become isotropic. During this process the eddies also lose their kinetic energy and in relation to that viscous forces get stronger. Finally, in section III of the energy cascade, viscous forces are now strong enough to dissipate the eddies into heat energy, denoted by ε. This terminates the energy cascade as at some point small eddies can no longer exist due to viscous effects [11]. Kolmogorov first described the 13 2. Theory size of these dissipative eddies with the so-called Kolmogorov length η = ( ν3 ε )1/4 . (2.12) As equation 2.12 shows, this part of the spectrum only depends on the turbulent dissipation ε and the kinematic viscosity ν. Between these two scales one finds the inertial subrange (II), which is only influ- enced by the cascade process itself [4]. To exist this region requires a sufficient high Reynolds number, so that the flow is fully turbulent. Eddies in this region are in the mid-range. Also the turbulence is isotropic, just like region III. Looking at the whole energy cascade this region is responsible for the spectral transfer. In this region one can also find the so called cut-off, which separates the resolved structures from the SGS stresses [11]. The inertial subrange is also called the Kolmogorov −5/3-range. This is because the decay of the energy spectra in this region is described by the following law, first presented by Kolmogorov [11]: E(κ) = CKε 2 3κ − 5 3 w (2.13) In this equation CK represents the Kolmogorov constant, which is about CK = 1.5 and κw is the wavenumber. One also wants to refer to the Kolmogorov hypothesis, where Kolmogorov pointed out the characteristics stated in this section [28][29]. 2.2.5 Turbulence Modeling In the previous sections 2.2.3 and 2.2.4 one introduced the basic principles and gov- erning equations of the CFD methods used in the present study. Moreover, one presented the unknown term in RANS equations, called the Reynolds stress tensor, or, in case of the spatial filtered equations, the SGS stresses. Finally, one discussed the energy density spectrum of turbulence, explaining the different turbulent scales, behaviour and interactions. These discussions represent the physical point of view. In contrast to that, one presents the modeling approaches for these physical char- acterstics in this section. Turbulence Modeling in RANS First, turbulence models for the RANS approaches is investigated. Therefore, one can distinguish between different models based on the way they deal with the Reynolds stress tensor. The first approach is to model the unknown term τij and therefore make some kind of assumption. In most cases this is the Boussinesq as- sumption. These approaches are generally known as models of first order. The second option is to solve the respective transport equation for the Reynolds stresses. This leads to six additional differential equations. Here it is not necessary to use the Boussinesq assumption, since the Reynolds stresses are actually calculated and not modeled. This is commonly known as models of second order such as the Reynolds 14 2. Theory Stress Models (i.e. the Algebraic Stress Model). This can produce better resolved turbulence structure and therefore more accuracy in the solution. On the downside more equations have to be solved, more terms modeled and the computer power needed can increase. However, in the present study only models of first order are used which is why they are going to be in the focus here [11][39]. The models of first order can be distinguished even further by the number of differ- ential equations needed to solve the problem. First, there are models using algebraic relations to describe the Reynolds stresses. Thus, they do not need more differential equations and therefore these models are called zero-equations-models, e.g. the mix- ing length model. Second, there are one-equation-models, introducing only one more differential equation into the problem. A very popular one-equation-model is the Spalart-Allmaras model [43]. Third and one of the most commonly used methods are the two-equation-models. Very popular two-equation-models are the k-ε model and the k-ω model and also various modifications of these two [49]. A fundamentally different approach is the LES, which was introduced in section 2.2.4. Finally, the DNS do not use any turbulence modeling at all, since every scale is resolved. An overview over the most common methods for turbulence modeling is shown in figure 2.7. Figure 2.7: Overview over commonly used turbulence models. Figure adapted from [39] In the following one gives a brief discussion of the different turbulence models and the general process of turbulence modeling. Therefore, the modeling of the Reynolds stresses of the respective approach are assessed. Most of the turbulence models use the Boussinesq assumption. Boussinesq proposed the idea, that Reynolds stresses behave proportional to the mean strain [39]. The basic idea is to replace the unknown 15 2. Theory Reynolds stresses τij with a turbulent viscosity νt. For the stress tensor one can now write v′iv ′ j = −νt ( ∂vi ∂xj + ∂vj ∂xi ) + 1 3δijv ′ iv ′ i = −2νtsij + 2 3δijk. (2.14) This means, that the six unknown Reynolds stresses are replaced by only one new unknown variable, the turbulent viscosity νt (also known as eddy viscosity). This viscosity is fundamentally different than a molecular viscosity, since it is not de- pendent on the fluid nor temperature, but only on the location [11]. Moreover, δij denotes the Kroenecker’s δ, which is defined as δij = 1 i = j 0 i 6= j . (2.15) This makes sure, that the normal components of the stresses are allocated one equal third each and each Reynolds stress influences the equations in an isotropic manner [11]. Note that the Boussinesq assumption introduces a vast simplification. The turbulent viscosity now needs to be calculated. It is defined as νt = utl. (2.16) ut represents the turbulent velocity and l the turbulent length scale [39]. The tur- bulent velocity describes the typical velocity within the energy containing eddies while the turbulent length scale describes a typical length for these eddies. How to calculate these variables determines the respective turbulence model. The one-equation models introduce one additional differential equation to determine the turbulent viscosity. In the beginnings, often the transport equation for the turbulent kinetic energy k was used [26]. This equation is derived from the NS- equation [11]. Because its importance for understanding turbulence phenomena like the energy density spectrum, the different terms of the equations are discussed here. For an incompressible flow with constant viscosity the transport equation reads ∂vjk ∂xj︸ ︷︷ ︸ I = −v′iv′j ∂vi ∂xj︸ ︷︷ ︸ II − ∂ ∂xj [ 1 ρ v′jp ′+ 1 2v ′ jv ′ iv ′ i− ν ∂k ∂xj ] ︸ ︷︷ ︸ III − ν ∂v′i ∂xj ∂v′i ∂xj︸ ︷︷ ︸ IV . (2.17) In this equation term I represents the convection term. Term II describes the production of the turbulent kinetic energy. Including the minus sign it is always positive, presenting the energy that the turbulence extracts from the mean flow. The turbulent diffusion is represented by term III. One can distinguish between the pressure-, velocity fluctuation- and molecular diffusion. Finally term IV includes the turbulent dissipation ε into the equation. In contrast to the production term, it is always positive, excluding the minus sign. Therefore, it represents the conversion of the turbulent kinetic energy to thermal energy, especially taking place at smaller scales [10]. 16 2. Theory However, to actually make use of this equation a number of unknown quantities need to be determined. Using the Boussinesq assumption for the production term II the Reynolds stresses can be substituted [10]. This leads to Pk = νt ( ∂vi ∂xj + ∂vj ∂xi ) ∂vi ∂xj . (2.18) Next the diffusion term needs to be determined. Therefore the triple correlation is modeled by a gradient law. It is assumed that k is diffused from regions with higher k to a region with lower k, very similar to the transport of heat [10]. One can substitute the correlation as 1 2vjvivi = −νt σt ∂k ∂xj . (2.19) σt denotes the turbulent Prandtl number [11]. It is defined as σt = νt αt . (2.20) It represents the relation of the turbulent viscosity to the turbulent thermal diffusion αt. With DNS and experimental results it is often determined that 0.7 ≤ σt ≤ 0.9 [11]. The diffusion due to pressure is very small and therefore neglected [10]. Lastly the turbulent dissipation term is changed to ε = k 3 2 l . (2.21) The final modelled k equation now reads ∂vjk ∂xj = Pk + ∂ ∂xj (( ν + νt σt ) ∂k ∂xj ) − ε, (2.22) For convenience one refers to this equation as the transport equation for k. This equation is used to in the one-equation- and some of the two-equation-models as well [11]. Similar to the zero-equation-models, the turbulent length is needed here as well. This is why the original one-equation models did not make major improvements to turbulence modeling [10]. A more promising approach was presented with the Spalart-Allmaras model, which uses a transport equation for the eddy viscosity in- stead [26][43]. However, the most popular RANS turbulence models are the two-equations-models, especially the k-ω and the k-ε model. Since the latter of these two is used in this thesis one concentrates on that one. When using the k-ε model one calculates the eddy viscosity with replacing the characteristic length and velocity with νt = cµ k2 ε , (2.23) 17 2. Theory where cµ is one of the model constants [11]. The modelled transport equation for k has been derived above, now the modelled ε equation is presented [10]. Although the exact transport equation for ε can be derived from the NS-equations, the number of unknown variables and terms would be very large and therefore it is inconvenient to use. Thus, the transport equation often used for modelling is based on physical reasoning. It reads ∂vjε ∂xj = ∂ ∂xj (( ν + νt σε ) ∂ε ∂xj ) + ε k (cε1Pk − cε2ε). (2.24) For solving all relevant equations in total five variables now have to be determined. While there can be found certain relations for three of these variables (cµ, cε1, cε2), the two remaining variables(σk and σε) can be optimized according to the certain flow, e.g. channel flow, pipes or jets etc. The standard values are shown in table 2.1 down below. Table 2.1: Typical model constants for the k-ε model [10] cµ cε1 cε2 σk σε 0.09 1.44 1.92 1.0 1.31 Since the k-ε model is widely used in the industry, a lot of effort has been put in optimizing its performance. When using the wall law, one can significantly re- duce the number of grid points and therefore save computational time and costs. Furthermore, the model constants have been investigated in many applications and experiments providing a good accuracy for many cases. Lastly, the eddy viscosity is always positive. Thus, the numerical stability is very good [26]. However, there are cases, where the k-ε model does not provide accurate results. In general it performs poorly, when the eddy-viscosity assumption (Boussinesq) is not valid. Moreover, it has its difficulties, if the anisotropic Reynolds stresses gain in influence, since due to the Boussinesq assumption only the isotropic stresses are taken into account. Flows, which are not in local equilibrium might be resolved poorly as well, because Pk 6≈ ε. Lastly, the k-ε model might has difficulties with flows, where the approximation k/ε and k3/2/ε does not fit time and length scales. In general flows which large acceleration or deceleration, strong swirl or strongly curved streamlines are difficult to resolve correctly with the eddy viscosity approaches, like the k-ε model [49]. However, the k-ε model in its original form is typically not in use anymore. To address its weaknesses stated above, several advanced models have been developed. One of them is the Shear Stress Transport (SST) model, proposed by Menter [32]. It combines the strengths of the k-ε and k-ω model by applying the a blending func- tion. Near the wall the k-ω model is used and moving into the free flow the k-ε model is applied. This leads to improvements especially for zero pressure gradient, and adverse pressure gradient boundary layers and free shear layers [49]. 18 2. Theory Another important variation is the realizable k-ε model, which is a non-linear model. This variation is also used in the present study. It modifies the anisotropic part of the Reynolds stresses without introducing the additional transport equation, like the turbulence models of second order (Reynolds stress models) do it. Moreover, the approximation of the eddy viscosity is changed as well. Thus, the realizable model considers the relaxation time, the turbulence needs to adjust itself to the flow domain. The Reynolds stresses now depend partially on the mean strain rate. This is a major improvement to the standard k-ε model, which is unable to reproduce physical flow structures, which change very quickly in the domain. While in the original k-ε model it is assumed that Pk ≈ ε at all times, in the realizable k-ε model the turbulence is allowed to adjust itself and disturb this balance. This leads to better performance for instantaneous changes in the domain [39]. Turbulence Modeling using LES Very similar to the RANS approach, the respective LES models can be distinguished from their modelling of the SGS stresses (equation 2.11). The most common ap- proaches are discussed in this section. First one substitutes the SGS stresses with an eddy viscosity approximation, similar to a physical viscosity or the previous pre- sented turbulent viscosity [26]. It yields τ̃ij = −2νeD̃ij. (2.25) νe is denoted as the SGS eddy viscosity coefficient [26]. One now substitutes the SGS stress tensor in equation 2.9 with equation 2.25. The momentum equation then reads ∂ṽi ∂t + ∂(ṽiṽj) ∂xj = −1 ρ ∂p̃ ∂xi + ∂ ∂xj [ 2(ν + νe)D̃ij ] . (2.26) Solving the filtered continuity and momentum equation (2.8 and 2.9) and providing the SGS eddy viscosity coefficient νe builds the basis for one of the most commonly used models, the Smagorinsky model [26]. For the eddy viscosity one can substitute νe = (Cs∆)2|D̃|, (2.27) where ∆ denotes the filter width, Cs is called the Smagorinsky constant and |D̃| represents the norm of the rate-of-strain tensor |D̃| = √ 2D̃ijD̃ij. (2.28) The Smagorinsky constant is the only non-dimensional constant that needs to be given for this model. With the assumption of local equilibrium and the Kolmogorov spectra one finds the theoretical value Cs = 0.173 [26]. In general the Smagorinsky model provides good results for the energy dissipation under the condition of an appropriate Smagorinsky constant. It matches experimen- tal results well, particularly for isotropic turbulence. Nevertheless, the Smagorinsky constant needs to be adapted for each scenario. Hence, it is not universal. Moreover, 19 2. Theory it does not perform well on predicting the turbulent behaviour in the subgrid scales. Accurate data with DNS showed, that the modelled SGS stresses differ from the actual SGS stresses obtained by these simulations. Furthermore, the Smagorinsky model can only transfer energy from larger to smaller eddies, because νe is always positive. Thus, it cannot reproduce the inverse cascade (backward scatter), which can exist in reality. For practical use the high robustness of the Smagorinsky model is a very important factor [26][49]. To address the downsides of the standard Smagorinsky model, different variations have been developed. One of the most commonly used approaches are the dynamic models. These methods use the local grid-scale turbulence mechanism additionally to the grid-scale velocity gradients to resolve the SGS eddy viscosity. In the stan- dard Smagorinsky model the constant Cs was chosen invariable. In the dynamic approaches Cs changes dynamically depending on the grid-scale velocity field. This approach is known as the dynamic Smagorinsky model (DSM) [26]. Major improvements of the DSM are treatment of laminar to turbulent transition, the treatment near the wall and the inverse energy cascade, which is now possible within the model. However, the DSM suffers from numerical instability [26]. Besides the RANS and LES approaches the so called hybrid LES/RANS method is also a promising approach. This method combines RANS and LES by treating the near wall flow with the URANS approach and the free flow as the LES. The main justification for these models is the reduction of cells, especially in near wall regions. In standard LES, the grid near the walls must be very fine in order to capture the physics. This is because the scales of the turbulent eddies are very small next to the wall and getting even smaller for larger Reynolds numbers. in hybrid approaches, RANS wall functions are applied in that regions, so that the mesh can be coarsened. In the outer regions LES is used [11]. 2.2.6 Modeling of the Atmospheric Boundary Layer The modeling of the ABL is a complex task. There is a huge variety of scales and the computational domain often needs to be very large. In order to describe the physical quantities in the ABL, different approaches were developed. In the present study one focuses on the modeling approaches for the neutral ABL. A commonly used method for ABL modeling with the k-ε model was introduced by Richard and Hoxey [38] and is often referred to as the logarithmic law. The logarithmic law is derived with physical and dimensional considerations. Among others, the Prandtl mixing length is used (very similar to the zero-equation turbulence models). One obtains an equation for the vertical gradient of the wind speed. Finally, with the integration one derives the mean velocity distribution [13][50]. It reads u(y) = u∗ κ ln y0 + y y0 . (2.29) y0 denotes the roughness length, which depends on the surface characteristics and y is the distance perpendicular to the surface. For offshore environment the value is 20 2. Theory very small, because the surface is very smooth [47]. Moreover, κ represents the von Kármán constant. Finally, u∗ is the so called friction velocity. It is defined as u∗ = κuref ln y0+yref y0 (2.30) with the reference velocity uref at the reference height yref . Furthermore, to complete the boundary conditions for ABL modeling, the turbulent quantities need to be defined. The turbulent kinetic energy reads k = u2 ∗√ cµ , (2.31) while the equation for the turbulent dissipation yields ε = u3 ∗ κ(y0 + y) , (2.32) with the model constant cµ. These boundary conditions fit the k-ε turbulence model, meaning they fulfill the transport equation for k and for ε. 2.2.7 Turbulent Quantities Finally, the physical quantities, that have not been introduced yet are discussed in this section. First and foremost the Reynolds number, as it is one of the most im- portant non-dimensional parameters for viscous flows. It describes the ratio between the inertial and the viscous term and is defined as Re = UL ν . (2.33) U and L describe the characterstic velocity and length, respectively [26]. These are characteristic values for each flow domain. In the present study one assumes the characteristic velocity to be 10 m s and the characteristic length to be 1000 m, as this is the height of the Ekman sublayer investigated in this study. One calculates the approximate Reynolds number for the present study as Re = 10 m s · 1000 m 1.45 · 10−5 m2 s = 6.9 · 108. (2.34) The kinematic viscosity was calculated for an altitude of 90 m and a temperature of 287.565 K [52]. The flow over a flat plate is considered to be fully turbulent above Reynolds numbers of ReL = 5 · 105 (where L is the distance from the leading edge) [20]. Therefore, the flow investigated in the present study is fully turbulent [49]. Besides the impact of the Reynolds number on the physical behaviour of the flow, it has also a large influence on the CFD simulations. The larger the Reynolds number is, the wider the range of length scales gets in turbulent flows. Therefore, resolving all scales in high Reynolds number flows becomes infeasible at some point. This is one of the major justifications for the use of turbulence models [26]. 21 2. Theory Moreover, the energy referred to the velocity fluctuations is quantified by the tur- bulent kinetic energy k. It reads k = 1 2v ′ iv ′ i = 1 2(Rexx +Reyy +Rezz), (2.35) where Rexx, Reyy and Rezz denote the normal components of the Reynolds stresses [51]. They can be calculated according to Reii = v′2i , (2.36) as they equal the variances of the turbulent velocity fluctuations, with i = 1, 2, 3 for the streamwise, wall-normal and spanwise direction, respectively. Finally, one introduces the turbulence intensity I in the respective direction. It is defined as Ii = √ v′2i vi . (2.37) It represents the turbulent fluctuations with regard to the mean velocity [51]. The turbulence intensity is often needed to set the boundary conditions of the CFD calculation, but also for the physics of wind turbines. The higher the turbulent fluctuations the higher the induced stresses on the turbine. Therefore, it is important to predict the turbulence intensity. 2.3 Wind Turbine Aerodynamics and Loads In this section the basic principals of wind turbine rotor aerodynamics are discussed. Furthermore, one introduces the Actuator Disk momentum theory, which builds the foundation for the Actuator Disk Model (ADM). Moreover, one presents the blade element method and combines these two approaches to the Blade-Element- Momentum (BEM) theory. The BEM theory is used within the Actuator Line Model (ALM) and together with the ADM both of these models are used within the commercial CFD solver Star-CCM+ to simulate ABL and rotor aerodynamics. Finally, one presents the aero-elastic solver FAST. 2.3.1 Aerodynamic Characteristics The main intention for using a wind turbine is to extract kinetic energy from a wind flow and transform it first to mechanical and then to electrical energy. Inevitably a force acting on the rotor blades is induced by the flow. For the assessment of these forces, one looks at the individual sections of the blade, which consist of different airfoils. The resulting force acting on an airfoil can be split into the force created by pressure (normal force) and the tangential component, which is called friction force. Furthermore, the forces are often divided into a component perpendicular to the flow direction (lift) and parallel to the flow (drag) [8]. 22 2. Theory Different coefficients describe these forces and aerodynamic behaviour of an airfoil. The lift coefficient is defined as Cl = L 1 2ρV 2 0 c , (2.38) where L denotes the lift force per span, ρ is the density, V0 is the wind speed in the free stream and c is the chord length2. The lift coefficient is a function of the Reynolds number or Mach number and the angle of attack (AOA). However, the Reynolds number usually affects the state of the boundary layer and therefore the risk of flow separation, which is mainly important for the maximum lift coefficient. Higher Reynolds numbers lead to higher stall AOA and higher maximum lift co- efficients. In the usual operating range, the influence of the AOA is much more important. Therefore, the lift coefficient is usually stated as a function of the AOA. Another important aerodynamic characteristic is the airfoil drag, which is presented by the drag coefficient Cd. It is defined as Cd = D 1 2ρV 2 0 c , (2.39) where D denotes the drag force per span experienced by the airfoil. Similar to the lift coefficient the drag coefficient is often plotted against the AOA. However, it can also be shown as a function of the lift coefficient. 2.3.2 Actuator Disk Momentum Theory In the following the Actuator Disk momentum theory is derived. The momentum theory is a simplified approach for assessing the momentum balance on a rotating stream tube. Therefore, it represents the actual turbine rotor as a circular disk, where friction is neglected [8][33]. This disk is extended both upstream and down- stream, creating a stream tube as shown in figure 2.8. Figure 2.8: Illustration of the Actuator Disk model. Figure taken from [8] 2The chord length is the shortest distance between the trailing and leading edge of an airfoil. 23 2. Theory The stream tube forms a control volume with cylindrical shaped ends, for which the momentum balance is assessed. There are four important sections of the stream tube. The first one is far upstream in the free stream region, two sections are right before and after the disk (rotor) and the last one is far downstream in the far wake region. The flow enters the stream tube with the axial velocity V0 and the pressure p0, representing the free stream conditions. When passing the disk, the flow experiences a pressure drop ∆p and a reduction of the axial velocity V as well, because kinetic energy is taken from the flow by the disk. When the flow leaves the stream tube it has the axial velocity V1 and the pressure p0. The axial velocity and static pressure as a function of the axial distance x is drawn in the figure 2.9. Figure 2.9: Axial velocity and static pressure curve of the flow over an Actuator Disk. Figure taken from [8] One can see, that the velocity is reduced after the disk, but the pressure recovers to the free stream pressure p0 in the far wake region. Another important character- istic is the velocity gradually decreasing before the flow reaches the turbine. This is due to the extraction of kinetic energy by the turbine. Additionally, the rotor induces a radial velocity because of the flow expanding upstream of the turbine. In contrast to that, the pressure increases slightly just in front of the rotor because of the stagnation of the velocity and therefore transformation from dynamic pressure to static pressure. Then a sudden pressure drop is caused by the disk, which can be explained by the basic momentum theory. After the disk the pressure recovers to the free stream pressure p0 [8]. The governing equation of the Actuator Disk momentum theory are now derived on the basis of Chen et al. [8] and Schaffarczyk [41] and written in the way they are implemented in the CFD software Star-CCM+ [1]. When applying the axial momen- tum theory to the control volume mentioned above, one can state two formulations 24 2. Theory for the change in axial force (thrust), that reads dT = dm(V0 − V1) = ρV (V0 − V1)dA, (2.40) and dT = ∆pdA, (2.41) where dm represents the mass flow rate in the control volume and dA represents the change of the section area of the stream tube. Assuming incompressible and steady flow allows the application of Bernoulli’s law before the turbine and after the turbine. Because of the extraction of kinetic energy the total pressure decreases after the rotor. One obtains a relationship for the pressure drop ∆p. It reads ∆p = 1 2ρ(V 2 0 − V 2 1 ). (2.42) Combining equation 2.40, 2.41 and 2.42 yields V = V0 + V1 2 . (2.43) For convenience the axial induction factor a is introduced3. The equations for the velocities then read V = V0(1− a) and V1 = V0(1− 2a). (2.44) Inserting equation 2.44 into equation 2.40 gives dT = 4πρV 2 0 a(1− a)rdr. (2.45) Substituting the dimensionless thrust coefficient Ct = 4(1 − a)a [41] gives another formulation, which reads dT = 1 2ρV 2 0 Ct2πrdr (2.46) or for the total thrust of the wind turbine T = 1 2ρV 2 0 Ctπ(r2 out − r2 in), (2.47) where rout and rin denotes the outer and inner Radius of the disk, respectively. Although in theory the actuator disk is infinitesimally thin, the source therms are added for a actuator disk thickness δ. The thrust for a certain ∆r annulus can be written as ∆T = 1 2ρV 2 0 Ct2πr∆r. (2.48) 3The axial induction factor a, sometimes also called axial interference factor, describes the velocity induced on the profile by the free stream. According to the Betz-Joukovsky theory the power and thrust coefficient depend on a. There are maximums of both coefficients, which is called the Betz-Jokouvsky limit [41]. 25 2. Theory Using equation 2.47 yields ∆T = T 2r∆r r2 out − r2 in . (2.49) The thrust per per unit volume within a certain annulus then reads ∆T ∆V = T δπ(r2 out − r2 in) = 1 2 ρV 2 0 Ct δ , (2.50) with the unite volume of ∆V = 2πr∆rδ. For each cell that yields Tcell = ∆T ∆V Vcell. (2.51) Finally equation 2.51 needs to be scaled, in case the cell is not fully within in the disk. This leads to Tcell = T Vcell∑ allcells Vcell . (2.52) A similar approach is used for the calculation of the torque distribution. The torque acting on an element for an annular ring with the width dr and the induction angular velocity ω is obtained with the conservation of angular momentum. It reads dQ = dm(ωr)r = 2πρV ωr3dr. (2.53) Additionally to the axial induction factor a, the tangential induction factor b is defined as b = ω 2Ω , (2.54) where Ω denotes the rotational rotor speed. Inserting the induction factors and integrating equation 2.53 leads to the total torque Q = b(1− a)1 2ρV0Ωπ(r4 out − r4 in). (2.55) The torque per unite volume is defined as ∆Q ∆V = 2Wr2 δπ(r4 out − r4 in) . (2.56) Similar to the thrust, the torque for each cell then reads Qcell = ∆Q ∆V Vcell = 2Qr2 cellVcell δπ(r4 out − r4 in) . (2.57) The final equation for the torque, with applied scaling yields Qcell = Qr2 cellVcell∑ allcells r 2 cellVcell . (2.58) 26 2. Theory The tangential force for each cell is then defined as Ftangential = QrcellVcell∑ allcells r 2 cellVcell . (2.59) The thrust (axial force, equation 2.52) and the tangential force (equation 2.59) are now used to create the source term, which is added to the momentum equations [1]. 2.3.3 Blade Element Method and Actuator Line Model The blade element method is an approach to model rotor effects and especially to resolve the wake flow field without describing the exact geometry of the rotor blades. Body forces are modelled along lines, which are representing the blades [1][33]. Applying the blade element method is a two steps process. First, one needs to mark the cells that are used as a virtual disk and where the source term is added and second, the source term for the respective cell needs to be calculated [1]. The marking and allocation process is shown in figure 2.10. The red grid is the initial finite volume mesh, whereas the black grid is an interpolation grid. The marking takes place in two steps. In the first step, all cells in the volume mesh are marked, which are lying in the disk area. Then, each marked cell is assigned to an element of the interpolation grid and the location is stored with respect to the interpolation grid. One can see, that after these two steps all of the marked cells belong to one of the elements of the interpolation grid. This also means, that there has to be at least one cell in each element of the interpolation grid, which can cause problems especially in the inner area [1]. Figure 2.10: Cell marking within the blade element method. Figure taken from [1] After marking the cells the source term is calculated for each element of the inter- polation grid and afterwards transferred to the volume mesh. During this process 27 2. Theory different coordinate transformations need to be conducted. Besides the global coor- dinate system of the whole domain, a local Cartesian coordinate system (i.e. x’, y’, z’) is created in the center of the disk. The z’-direction is perpendicular to the disk, pointing in the thrust direction. Moreover, the disk rotates around the z’ axis. In case of rotor blade flapping two additional coordinate systems are created, a cylin- drical system with respect to the base of the local Cartesian coordinate system and a so called flap coordinate system, which is another Cartesian coordinate system with respect to the base of the previous created cylindrical coordinate system. This way the rotation of the disk and each blade is fully described[1]. All calculations are done within the flap coordinate system. Afterwards they are transformed to the global coordinate system and finally added to the momentum equations. For this the angle of attack and the relative velocity experienced by the blade is determined. Then, the lift and drag force is calculated according to L = 1 2ρ|vrel|2Clc dr (2.60) and D = 1 2ρ|vrel|2Cdc dr, (2.61) where ρ, vrel and c denote the density, relative velocity and chord length, respectively. The aerodynamic coefficients for lift and drag are denoted by Cl and Cd, respectively and dr represents the with of the individual element. For further improvement one can also add the tip-loss correction. However, since it was not used in the present study, it is not further explained here. Finally, the lift and drag force are transformed to the normal and tangential component [1]. It yields Fnorm = L cos β −D sin β (2.62) and Ftang = L sin β −D cos β, (2.63) where β denotes the inflow angle. Now, the force vector is assembled and trans- formed into the global coordinate system. Before adding the force vector to the momentum equations, it is time-averaged over each blade element [1]. 2.3.4 Blade-Element-Momentum (BEM) Theory The BEM theory consists of the combination of two methods to assess the dynamic response of a turbine in a flow field. The Actuator Disk momentum theory is one approach and was introduced in section 2.3.2. The other one is the Blade Element theory, which assesses the aerodynamic loads on the blade at different sections. To- gether these two approaches provide four equations for the BEM method, which is then solved iteratively [8][21]. 28 2. Theory Two equations come from the momentum theory. The axial force (thrust) is defined as in equation 2.45, whereas the torque is given as dQ = 4b(1− a)ρV0Ωr3πdr. (2.64) Note, that this formulation equals equation 2.55 in its differential form. Further- more, the tip-loss correction is not considered here. Additionally, the blade element theory delivers another two equations for the thrust and the torque, expressed by lift and drag coefficients. The thrust is defined as dT = σπρ V 2 0 (1− a)2 cos2 β (Cl sin β + Cd cos β)rdr, (2.65) whereas the equation for the torque reads dQ = σπρ V 2 0 (1− a)2 cos2 β (Cl cos β − Cd sin β)r2dr. (2.66) Note that σ denotes the solidity σ = Bc 2πr , (2.67) where B, c and r denote the number of blades, the chord length and the radius, respectively. These four equations (2.45, 2.64, 2.65 and 2.66) are now combined to receive the relationships a 1− a = σ(Cl sin β + Cd cos β) 4 cos2 β (2.68) and b 1− b = σ(Cl cos β − Cd sin β) 4λ cos2 β , (2.69) where λ denotes the tip-speed ratio (TSR) λ = Ωr V0 . (2.70) A typical iteration process goes as follows. To begin with the iteration, values for the induction factors a and b are guessed. By solving the equations 2.68 and 2.69 one calculates the inflow angel β and the TSR λ. The aerodynamic coefficients Cl and Cd are then determined using tabulated airfoil tables for the individual profile. This gives information about the aerodynamic performance of the turbine in this particular operating point. Finally, the induction factors are calculated again to start a new iteration [21]. 29 2. Theory Aero-Elastic Assessment In the previous sections 2.3.2, 2.3.3 and 2.3.4 the theory used for the aero-elastic assessment was discussed. Within Star-CCM+ the ADM and ALM are used to calculate turbine forces such as thrust and torque. Although the ALM and FAST both use the BEM theory to calculate the aerodynamic forces, there are various differences in the actual application. These are now discussed and it is outlined, why using FAST has several advantages. First, the rotation rate is fixed in Star-CCM+. In contrast to that FAST has an included control system, adapting the rotation rate and other turbine characteris- tics to the current operating point. This is especially important when simulating complex or extreme situation, where the control system plays an important role. Furthermore, the effect of the tower and nacelle can be simulated in FAST but not in Star-CCM+. However, nowadays it is possible to model these effects within the ADM and ALM in general [44]. Finally, FAST works with an instantaneous field, while the ADM and ALM are mostly designed for a steady wind field. In addition to that, FAST is a popular tool for assessing wind turbine loads and therefore it is well documented and provides reference cases and values for various applications. FAST consists of different modules, which are briefly described here. As an input one needs to provide a flow field. This is transmitted through the InflowWind mod- ule to the AeroDyn module. Here the forces and moments are calculated according to the BEM theory. An additional module called HydroDyn is available to include wave loads for offshore applications. The resulting dynamic stresses for drivetrain, rotor, nacelle and tower are then determined within the ElastoDyn module. Similar to that the substructure is calculated in the SubDyn module. Finally, the ServoDyn module calculates power outputs. 2.4 Reference Turbine For assessing the performance and accuracy of the conducted calculations in the present study one chose the NREL offshore 5-MW baseline wind turbine developed by the National Renewable Energy Laboratory (NREL) as a reference turbine. It is widely used as a standard reference for developing and improving concept studies for wind technology. The most important characteristics for this study are shown in table 2.2 according to the documentation from the NREL [23]. Table 2.2: Characteristics of the NREL offshore 5-MW baseline wind turbine [23] Rating 5 MW Rotor Orientation, Configuration Upwind, 3 Blades Control Variable Speed, Collective Pitch Rotor, Hub Diameter 126 m, 3 m Hub Height 90 m Cut-In, Rated, Cut-Out Wind Speed 3 m/s, 11.4 m/s, 25 m/s Cut-In, Rated Rotor Speed 6.9 rpm, 12.1 rpm 30 2. Theory The turbine was created as a conventional three-bladed upwind variable-speed and variable blade-pitch-to-feather-controlled turbine. In addition to the general char- acteristics the steady-state behaviour is documented as well. It is shown in figure 2.11. Figure 2.11: Steady-state behaviour of the NREL offshore 5-MW baseline wind turbine. Figure taken from [23] One can see the operating range reaching from 3 m/s wind speed (cut-in) to 25 m/s (cut-out). Outside of this range the wind is either not strong enough to sufficiently generate power, or it is too strong so the stresses on the turbine would become too high. When looking at the green curve, one can see the generated torque, which is increasing quadratically with the wind speed until it reaches its maximum and held constant at the rated wind speed of 11.4 m/s. Then the active pitch control regulates the rotor speed, which is displayed by the red curve. Furthermore, one can see that the rotor speed is held constant at velocities higher than the rated velocity due to the pitch control. The TSR is constant in Region 2 to ensure a good efficiency of the wind to power transformation and gradually decreases afterwards, because the rotor speed is constant in that region. 31 2. Theory 32 3 Methods This chapter describes the numerical methods used in the present project work. In the beginning, the computational domain is described, explaining the boundary conditions on the one hand and the computational discretization (the mesh) on the other. Furthermore, the simulation process is outlined. In total four different types of Computational Fluid Dynamics (CFD) simulations were conducted. First of all, a Reynolds-Averaged Navier-Stokes (RANS) simula- tion was implemented to generate the correct inlet conditions for a superior Large Eddy Simulation (LES). This is necessary to provide a velocity profile and turbulent quantities in order to generate correct boundary conditions for the LES. After the validation of the LES the Actuator Disk Model (ADM) and Actuator Line Model (ALM) were implemented. The ADM and ALM provide additional outputs, e.g. power and thrust of the wind turbine, which will be assessed and compared with the FAST calculations. Moreover, they deliver information about the wake of each turbine. Further on, the transient flow field generated with the LES was used for aero- elastic simulations with FAST. For comparison and validation another flow field was transmitted to FAST, which was created with the spectral method. These two simulations were investigated to outline differences between these models. 3.1 Computational Domain For all simulations the computational domain is set up as a rectangular box. It is shown in figure 3.1. In the x-direction (streamwise) the domain extends 6000 m, in the z-direction (spanwise) it measures 4000 m and the height is 1000 m (y- direction, wall-normal). Two turbines are investigated in the present study. The center of each turbine (hub height) is at the location [x, y, z] = [2000, 90, 2000] m and [x, y, z] = [4000, 90, 2000] m, respectively. This distance is chosen that the wake of the first turbine does not affect the second. Several probes are placed at these locations and at upstream and downstream positions, pointing in y and z direction. This will be discussed together with the results in chapter 4. The inlet surface is marked yellow, while the outlet region is coloured red. The bottom of the domain (gray) is set up as a flat plate, since the offshore environment is very smooth. The remaining three surfaces are symmetry planes. One discusses the boundary conditions of the individual simulations more in detail within the respective sections. 33 3. Methods Figure 3.1: Computational Domain used for the CFD simulation in Star-CCM+ 3.1.1 Mesh Overall, two different meshes are used. For the RANS simulation and LES the Trimmer and Prism Layer Mesher are selected. The base size of the cells is set to 20 m in both cases. This is also the maximum cell size in the domain. Different base sizes were tested and 20 m proved to be the most efficient option while retaining good accuracy. To refine the near wall region, a prism layer is constructed next to the plate. There are 80 prism layers, building a refined region from the plate up to an altitude of 200 m. The first cell size is very important in Atmospheric Boundary Layer (ABL) modeling. The first layer of cells is about 0.004 m thick (in y-direction), which shows good agreement to the recommended cell size of the Star-CCM+ user guide [1]. This leads to y+ values around 30. The total number of cells for the RANS simulation and LES is 7.2 million cells. Note that this mesh is not shown here due to its simplicity. Figure 3.2: Sideview of the mesh for the Actuator Disk Model and Actuator Line Model For the ADM and ALM approach the mesh is different, due to the requirements of these models. A cut through the y-z plane is shown in figure 3.2. The Surface 34 3. Methods Remesher is added as an additional mesher. Most importantly the region around the actuator disks needs to be refined. That is why two Volumetric Controls are set up for each disk. The cells adjacent to the disk are very refined cubic cells, with a cell size of 1.25 m. The Volume Growth Rate is set to very slow for a smooth growth and appropriate cell quality. A second refined region is defined to have accurate results in the wake flow field. The cubic cells in this region measure 5 m in edge length. The prism layer is significantly reduced for the ADM and ALM approaches. Since Volumetric Controls starts at a height of y = 25 m, the total thickness of the prism layer equals this value. The total number of prism layers is 30 m and the first cell height is around 0.125 m, leading to y+ values around 1000. The total number of cells is around 6.5 million. 3.2 RANS Simulation In this section one covers the setup of the conducted RANS simulation. As a tur- bulence model the realizable k-ε model is used. It is a steady, 3D simulation with a constant density. The initial conditions, reference values and fluid characteristics are defined as written in table 3.1. The density, viscosity and static pressure are chosen according to the standard atmospheric conditions at an altitude of 90 m (hub height) [52]. Table 3.1: Initial conditions, reference values and fluid characteristics for the RANS simulation Type Quantity Value Fluid characteristics Density ρ = 1.21 kg m3 Fluid characteristics Dynamic viscosity µ = 1.89 · 10−5 kg ms Reference Value Pressure p = 1 · 105 kg ms2 Initial condition Pressure p = 0 kg ms2 Initial condition Velocity [ux, uy, uz] = [0, 0, 0] m s Initial condition Turbulent kinetic energy k = 0.347 m2 s2 Initial condition Turbulent dissipation ε = 0.1 m2 s3 One now describes the boundary conditions. The inlet is set as a velocity inlet. The velocity profile is defined as introduced in chapter 2, equation 2.29, with the von Kármán constant κ = 0.42 [38] and a roughness length of y0 = 0.0002 m [47]. The friction velocity u∗ was calculated for a reference velocity of vref = 10 m s at the reference height of yref = 90 m (hub height). The turbulent conditions are specified trough the turbulent kinetic energy and the turbulent dissipation, see chapter 2 equations 2.31 and 2.32 [38]. Note that the inlet velocity and turbulent dissipation vary with the altitude, while the turbulent kinetic energy remains constant. The outlet is defined as a pressure outlet with the static pressure p = 0 kg ms2 . The definition of the plate is always crucial in ABL modeling. It was set as a wall with the no-slip condition and a rough surface. For the calculation of the wall function several constants need to be defined. They are EC = 9 and C = 0.253 [1]. 35 3. Methods Moreover, the so-called roughness height rh is calculated as rh = Ey0 C = 0.0071 m, (3.1) according to the Star-CCM+ user guide [1]. The remaining surfaces are set as symmetry planes. 3.3 LES Simulation The LES was performed on the basis of the precursor RANS simulation and is described in this section. As a turbulence model the Smagorinsky Subgrid-Scale model is chosen. Different models were tested and this one was found stable for this setup. It is an implicit unsteady, 3D simulation of second order, with a time step of ∆t = 0.05 s. This time step fulfills the CFL number on the one hand and is not too large to cause inaccuracy for the aero-elastic analysis with FAST on the other. The total simulation time is 2400 s, leading to a total number of 48000 time steps. The number of inner iterations is 5. The flow needs 600 s to pass the domain. After roughly 1.5 flow passes (1000 s or 20000 time steps), the data is expected to be converged and data is recorded for the following 1400 s. Mean values for the results are averaged over the this interval. The fluid characteristics remain similar to the RANS simulation, as well as the reference values. The boundary conditions are defined as follows. For the velocity inlet a velocity profile from the precursor RANS simulation is applied. Furthermore, the synthetic turbulence is specified with the turbulence intensity and turbulent length scale. The turbulence intensity is defined as I = 0.047, while the turbulent length scale equals l0 = 229 m. These values are obtained from the precursor RANS simulation at hub height (90 m). Similar as the RANS simulation, the outlet is defined as a pressure outlet, with a value of p = 0 kg ms2 . The bottom plate is defined as a wall with the no-slip condition. However, Star-CCM+ leaves no option for a rough wall in LES. Apart from that the other surfaces are defined as symmetry planes. For the creation of a transient flow field for the FAST simulations a presentation grid was set up to capture the flow around the two turbines. It range from y = 19.875 m to y = 161.125 m in the wall-normal direction and from z = 1929.4 m to z = 2070.6 m in the spanwise direction. Its resolution is 60x60 grid points. 3.4 Actuator Disk and Actuator Line Model For the ADM and ALM the identical setup as for the LES is used. That in- cludes the applied models, initial conditions, fluid characteristics and reference values, as well as the boundary conditions. However, two Actuator Disks are placed in the domain, with their center at [x, y, z] = [1999.375, 90, 2000] m and [x, y, z] = [3999.375, 90, 2000] m. Note that the locations differ slightly from the previously introduced locations. This is because the center of each turbine must not 36 3. Methods be between two cells. The differences in terms of the mesh were already presented, now the physical models for the AD and AL model are discussed. The ADM only requires few information. The performance of the wind turbine is specified as a function of the wind speed. Therefore, the generated power and thrust coefficients are defined within a table for the wind speed according the description from the NREL [23]. The disk geometry is specified with an inner radius of rin = 1.5 m and an outer radius of rout = 63 m, whereas the thickness is δ = 4.5 m [23]. The normal direction is perpendicular to the inlet surface, pointing in streamwise direc- tion. Moreover, the inflow for the disk needs to be specified. It was defined according to the Star-CCM+ manual [1], setting the inflow plane radius slightly larger than the actual disk (about 10 %), which leads to a radius of r = 70 m. The offset of the plane is 12.6 m, which is 10 % of the diameter as well. Finally, the rotation rate is defined as the operating point specified as 12.1 rpm [23]. The ALM requires more details of the rotor blades including airfoil properties. Therefore, the blade is divided into 18 nodes. For each node a table is provided, delivering the lift coefficient Cl and the drag coefficient Cd as a function of the angle of attack (AOA) and Mach number. Moreover, the chord distribution is tabulated as the chord length for each node over the normalized span (i.e. radius divided by the current location). Similarly, the aerodynamic twist is given for the normalized span. Similar to the ADM, all tables can be found in the description of the NREL [23]. As the ALM is originally developed for helicopter applications, a sweep angle can be specified, but for the current study it is a constant value of 0 radian. The same goes for the disk stick and disk flap specification. For the geometry the thick- ness is defined as 1.25 m, since at least 3 cells should cover the actuator disk [1]. Despite that, the geometry is defined identical to the ADM, as well as the rotation rate. The azimuthal disk resolution is set to 7 and the radial resolution to 49 [1]. No tip-loss correction is applied for the present study. 3.5 Aero-elastic Simulation Finally, the aero-elastic simulations using FAST are discussed. After performing the LES, the transient velocity field