Energy and Power Engineering
Vol.09 No.13(2017), Article ID:81413,15 pages
10.4236/epe.2017.913053

Effects of Terrain-Induced Turbulence on Wind Turbine Blade Fatigue Loads

Yasushi Kawashima1, Takanori Uchida2

1West Japan Engineering Consultants, Inc., Fukuoka-shi, Japan

2Research Institute for Applied Mechanics, Kyushu University, Kasuga-shi, Japan

Copyright © 2017 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

Received: November 15, 2017; Accepted: December 25, 2017; Published: December 28, 2017

ABSTRACT

Recently, the issue has surfaced that the availability factors for wind farms built on complex terrain are lower than the originally projected values. In other words, problems have occurred such as extreme decreases in generation output, failures of components inside and outside wind turbines including yaw motors and yaw gears, and cracking on wind turbine blades. As one of the causes of such issues, the effects of wind turbulence (terrain-induced turbulence) have been pointed out. In this study, we investigated the effects of terrain-induced turbulence on the structural strength of wind turbines through the measurement of strains in wind turbine blades and the analysis of wind data in order to establish a method for optimal wind turbine deployment that uses numerically simulated wind data and takes the structural strength of wind turbines into consideration. The investigation was conducted on Wind Turbine #10 of the Kushikino Reimei Wind Farm (in operation since Nov. 2012) in cooperation with Kyudenko New Energy Co., Ltd. Subsequently, we conducted numerical wind simulations (diagnoses of terrain-induced turbulence) to study the effects of the properties of airflow on the structural strength of wind turbines. For these simulations, the natural terrain version of the RIAM-COMPACT software package, which is based on large eddy simulation (LES), was used. The numerical simulations successfully reproduced the characteristics of the wind conditions and the structure of the three-dimensional airflow. These results enabled us to determine the threshold value for a turbulence index to be used for optimal wind turbine deployment planning that utilizes quantitative data from simulations with the natural terrain version of the RIAM-COMPACT software package.

Keywords:

Complex Terrain, Terrain-Induced Turbulence, LES

1. Introduction

In recent times, the availability factors of wind farms constructed over complex terrain in mountains and other areas have fallen below the initially projected values. That is, wind turbines with significantly reduced power output and damage to the interior and exterior of the wind turbines (e.g., breakdown of yaw motors and yaw gears and cracks on wind turbine blades) have surfaced as issues. As one of the causes for such issues, previous literature [1] [2] [3] has pointed out the effect of turbulence (terrain-induced turbulence) which originates from small terrain undulations in the vicinity of wind turbines.

When land-based wind turbines in Japan are constructed in the future, it can be anticipated that construction will continue to take place in complex terrain such as mountainous areas, as a result of the search for well-suited deployment sites. Therefore, highly-accurate evaluation methods for wind turbine deployment need to be established in order to reduce wind turbine accidents and failures.

Given this circumstance, our research group has conducted “highly accurate numerical wind simulations (numerical wind diagnoses) with the RIAM- COMPACT natural terrain version software package [1] ”. In the present study, in-situ wind turbine blade strain and wind data were collected for Wind Turbine (WT) #10 of the Kushikino Reimei Wind Farm, Kagoshima Prefecture, Japan. Based on the in-situ data, the wind direction is identified for which the structural strength of the wind turbine is most affected. In addition, the in-situ data are analyzed in detail. Subsequently, for the identified wind direction, numerical wind simulations (numerical wind diagnoses) based on large-eddy simulation (LES) are performed. The simulation results are examined in terms of their correlation to the in-situ data collected over the wind farm with particular attention to the characteristics as well as the three-dimensional flow structure of the wind. Finally, a planning method is proposed for optimal wind turbine deployment, a method which uses numerical simulation results and takes wind turbine structural strength criteria into consideration.

2. Overview of Kushikino Reimei Wind Farm

With the cooperation of Kyudenko New Energy Co., Ltd., the present study investigates Kushikino Reimei Wind Farm (in operation since Nov. 2012), which is located in Hatori District, Ichikikushikino City, Kagoshima Prefecture, Japan (Figure 1). At this wind farm, ten 2-MW wind turbines manufactured by Hitachi, Ltd. have been deployed. Of these turbines, the present study focuses on WT #10, for which there is a concern that it is affected by turbulence (terrain-induced turbulence) which is generated when wind flows over Mt. Benzaiten (Elevation: 519 m), a mountain located to the east of the wind farm (Figure 2 and Figure 3, Table 1).

3. Analysis of Wind Turbine Blade Strain Data (Analysis of In-Situ Data)

Since wind energy enters the wind turbine system via the turbine blades, assessment

Figure 1. Map of the Kushikino Reimei Wind Farm and the surrounding area (Google Earth).

Figure 2. Photo of WT #10.

Figure 3. Relative locations of Mt. Benzaiten (Elevation: 519 m) and WT #10 (Google Earth).

Table 1. Elevation information for WT #10 and distance between Mt. Benzaiten (elevation: 519 m) and WT #10.

Figure 4. Blade strain measurement positions on WT #10.

of the strength of the blade roots (assessment of blade bending load) is highly important when evaluating the structural strength of the wind turbine. Accordingly, electrical strain sensors were installed at the roots (root section: the section of the blade, approximately 1.3 m long, which extends outward from the blade-hub junction) of the three blades of WT #10 in the present study (Figure 4).

In-situ blade strain data from WT #10 are analyzed for the time period between 00:00 JST, Nov. 3, 2015 and 07:00 JST, Mar. 17, 2016, and the results of the analysis reveal that the magnitude of the fluctuations of blade strain was the largest on Nov. 13, 2015. Given this finding, with the cooperation of the wind turbine manufacturer, an additional analysis is performed with the wind turbine operation data and the in-situ blade strain data to investigate the effect of the fluctuations of blade strain on the fatigue load on the wind turbine blades. From the analysis, it is determined that for the case of the easterly wind in the 10-min period starting at 09:40 JST (10-min average wind speed: 9.1 m/s), the value of the damage equivalent load (DEL) for blade bending is 2.03 (equivalent to reaching the design load in 5.88 years under a continuously blowing wind with the same flow characteristics as those from the 10-min period) [4] (Figure 5).

In the present study, the in-situ blade strain and DEL data under easterly wind, in which the value of DEL becomes larger than that in wind flowing from any other direction, are compared to those under wind from other wind directions. Here, the results are shown of this comparison between a time period with northerly wind, which is wind from the most frequently-occurring wind direction during the entire in-situ data collection period (Figure 6, Table 2), and a time period with easterly wind. In both time periods, the average wind speed was approximately 9 m/s. The results reveal that DEL is 0.95 under northerly wind in contrast to 2.03 under easterly wind, indicating that there exists a distinct difference in DEL and the magnitude of the fluctuations of blade strain between the time periods with northerly and easterly wind (Figure 7, Table 4).

Figure 5. Time series of wind speed and normalized blade bending moment. Plotted values were evaluated from instantaneous data values within 10-min periods. (Analysis period: 08:00 - 13:00 JST, Nov. 13, 2015): (a) Wind speed and (b) Damage equivalent load for flapwise blade bending.

Figure 6. Frequency distribution of the direction of the 10-min average wind (%): (a) Wind rose (Frequency distribution) and the average of the 10-min average wind speed observed for 16 wind directions (m/s): (b) Wind speed by direction (Wind measurement height: hub height (60 m), analysis period: 00:00 JST, Nov. 3, 2015 - 07:00 JST, Mar. 17, 2016).

Figure 7. Blade strain data (Blade flapwise bending data) in the cases of (a) easterly and (b) northerly wind (Interval: 0.02 second, average wind speed: 9 m/s).

Table 2. Frequency distribution of the direction of the 10-min average wind (%) and the average of the 10-min average wind speed observed for 16 directions (Wind measurement height: hub height (60 m), analysis period: 00:00 JST, Nov. 3, 2015 - 07:00 JST, Mar. 17, 2016).

Table 3. Wind direction range and total number of data values.

Note: includes only data from 10-min time periods with an average wind speed of 4 m/s (cut-in wind speed) or higher.

4. Airflow Field Analysis Using Data from the Nacelle Prop-Vane Anemometer (In-Situ Data Analyses)

Subsequently, in order to evaluate the relationship between the results of the DEL analysis (Section 3) and three-dimensional airflow characteristics with the use of the in-situ data, the in-situ wind data from WT #10 are analyzed for the time period between 00:00 JST, Nov. 3, 2015 and 07:00 JST, Mar. 17, 2016, which is the period for which both blade strain and wind data were collected for this wind turbine. Specifically, 10-min average wind speed, standard deviation, and turbulence intensity are analyzed for the 16 wind directions. Due to limited space, only the results from the analyses of the standard deviation of the wind speed and turbulence intensity in northerly and easterly wind (Table 3) are shown in Figure 8 and Figure 9. For these analyses, only the data from 10-min time periods with an average wind speed of 4 m/s (cut-in wind speed) or higher are included.

Figure 8 and Figure 9 show that, for the class of winds with a 10-min average wind speed of 10 m/s or less, the values of the standard deviation of the wind speed and the turbulence intensity from the corresponding time period under easterly wind are much larger than those under northerly wind. Furthermore, in easterly wind, for the class of winds with a 10-min average wind speed of 10 m/s or less, many of the values of turbulence intensity exceed the value of turbulence intensity for IEC Category A turbulence [5] (Figure 9). As already shown in Section 3, the value of DEL from a period with easterly wind with an average wind speed of approximately 9 m/s is roughly twice as large as that from a period with northerly wind with an average wind speed of approximately 9 m/s (Table 4). The combined results of the in-situ data analyses from above reveal for the first time that there exists a very high correlation between DEL and the standard deviation of the wind speed.

It is speculated that the differences in the magnitude of the fluctuations of blade strain, the standard deviation of the wind speed, and the turbulence intensity that were found for WT #10 between the time periods with northerly and

Figure 8. Relationship between the standard deviation and the average of the wind speed in 10-min periods for two wind directions (Wind measurement height: hub-height (60 m), analysis period: 00:00 JST, Nov. 3, 2015 - 07:00 JST, Mar. 17, 2016): (a) Easterly wind and (b) Northerly wind.

Figure 9. Relationship between turbulence intensity and the average of the wind speed in 10-min periods for two wind directions (Wind measurement height: hub-height (60 m), analysis period: 00:00 JST, Nov. 3, 2015 - 07:00 JST, Mar. 17, 2016): (a) Easterly wind and (b) Northerly wind.

Table 4. Measured data analysis and LES simulation results (Summary).

Note: Turbulence evaluation index [=σu (standard deviation of the streamwise wind velocity at the hub height)/U (inlet boundary streamwise wind velocity (=10 m/s))].

easterly wind are attributable to the effect of Mt. Benzaiten (elevation: 519 m), which is located approximately 300 m to the east (bearing: 78 degrees) of WT #10 (Figure 2 and Figure 3).

5. Overview of Numerical Simulation Method

To quantitatively evaluate the characteristics of wind which affects the durability of wind turbine blades (blade bending), two numerical wind simulations are performed. Specifically, one is for a case of easterly wind, for which the effect of terrain-induced turbulence is speculated to be large, and the other is for a case of northerly wind, for which the magnitude of the fluctuation of blade strain and the values of DEL, turbulence intensity, and the standard deviation of the wind speed differed significantly from those in easterly wind (Sections 3 and 4).

For the numerical simulations, the RIAM-COMPACT® natural terrain version software package is used, for which a collocated grid in a general curvilinear coordinate system is adopted [1] [2] [3] . In this collocated grid, the velocity components and pressure are defined at the grid cell centers, and variables that result from multiplying the contravariant velocity components by the Jacobian are defined at the cell faces. For the numerical technique, the finite difference method (FDM) is adopted, and a large-eddy simulation (LES) model is used for the turbulence model. In the LES model, a spatial filter is applied to the flow field to separate eddies of various scales into grid-scale (GS) components, which are larger than the computational grid cells, and sub-grid scale (SGS) components, which are smaller than the computational grid cells. Large-scale eddies, i.e., the GS components of turbulence eddies, are directly numerically simulated without the use of a physically simplified model. In contrast, dissipation of energy, which is the main effect of small-scale eddies, i.e., the SGS components, is modeled according to a physics-based analysis of the SGS stress.

For the governing equations of the flow, a filtered continuity equation for incompressible fluid (Equation (1)) and a filtered Navier-Stokes equation (Equation (2)) are used. Since the present study investigates wind fields with a hub-height mean wind speed of 5.0 - 6.0 m/s, the effects of the vertical thermal stratification of the atmosphere (atmospheric stability) are neglected. In addition, the effects of wind turbine wakes are not taken into consideration. Regarding the effects of the surface roughness, they are also neglected in the present study for the following reason. When the surface roughness is distributed nearly evenly, terrain undulations affect the local wind more significantly than the surface roughness as discussed in the previous literature [1] [2] [3] . For the computational algorithm, a method similar to a fractional step (FS) method [6] is used, and a time marching method based on the Euler explicit method is adopted. The Poisson’s equation for pressure is solved by the successive over-relaxation (SOR) method. For discretization of all the spatial terms except for the convective term in Equation (2), a second-order central difference scheme is applied. For the convective term, a third-order upwind difference scheme is applied. The interpolation technique by Kajishima [7] is used for the fourth-order central differencing that appears in the discretized form of the convective term. For the weighting of the numerical diffusion term in the convective term discretized by third-order upwind differencing, α = 0.5 is used as opposed to α = 3.0 from the Kawamura-Kuwahara scheme [8] in order to minimize the influence of numerical diffusion. For LES subgrid-scale modeling, the standard Smagorinsky model [9] is adopted with a model coefficient of 0.1 in conjunction with a wall-damping function (Equations (3)-(8)).

u ¯ i x i = 0 (1)

u ¯ i t + u ¯ j u ¯ i x j = p ¯ x i + 1 R e 2 u ¯ i x j x j τ i j x j (2)

τ i j u i u j ¯ 1 3 u k u k ¯ δ i j 2 ν S G S S ¯ i j (3)

ν S G S = ( C s f s Δ ) 2 | S ¯ | (4)

| S ¯ | = ( 2 S ¯ i j S ¯ i j ) 1 / 2 (5)

S ¯ i j = 1 2 ( u ¯ i x j + u ¯ j x i ) (6)

f s = 1 exp ( z + / 25 ) (7)

Δ = ( h x h y h z ) 1 / 3 (8)

6. Overview of Numerical Wind Simulation Set-Up

In this section, the numerical wind simulation set-up for the case of easterly wind is described. The computational domain used in the present study extends over a space of 12.0 km (x) × 2.0 km (y) × 2.6 km (z), where x, y, and z are the streamwise, spanwise, and vertical directions, respectively (Figure 10). A buffer zone is established in the upstream end of the computational domain, in which the terrain irregularities are reduced by 97% to form flat terrain. Similarly, a buffer zone is added to the downstream end of the domain. The maximum and minimum elevations of the surface in the computational domain are 523 m and 0 m, respectively. For the terrain elevation data, 10 m resolution data from the Geospatial Information Authority of Japan are used. The total number of computational grid points including those in the buffer zones in the upstream and downstream ends of the computational domain are 2401 (x) × 401 (y) × 101 (z), that is, approximately 9.8 million. The grid spacing is uniformly 5 m in both the x- and y-directions. The vertical grid spacing decreases smoothly down to 1.5 m at the ground surface.

As for the boundary conditions that are adopted in the simulation, a vertical profile of wind speed is applied at the inflow boundary according to Roughness Classification III [10] . Turbulence fluctuations in the inflow wind are neglected in order to focus our investigation on the effect of terrain-induced turbulence which originates from terrain undulations. For the side boundaries and upper boundary, the free slip boundary condition is imposed, and a convective outflow condition is applied to the outflow boundary. On the ground surface, the

Figure 10. Computational domain.

Figure 11. Characteristic two scales (U and h).

no-slip boundary condition (the viscous boundary condition) is imposed. The non-dimensional parameter Re in Equation (2) is the Reynolds number (= Uh/ν). For the simulation, Re = 104 is used based on Kato [11] . Figure 11 illustrates the characteristic scales used in the present simulation: h is the difference between the minimum and maximum terrain elevations within the computational domain, U is the streamwise wind velocity at the inflow boundary at the height of the maximum terrain elevation within the computational domain, and ν is the kinematic viscosity of air. The time increment is set to Δt = 2 × 10−3 h/U. The numerical wind simulation set-up for the case of northerly wind is identical to that for the case of easterly wind.

7. Simulation Results and Discussions

Here, discussions proceed with a focus on WT #10, which was speculated to be affected significantly by terrain-induced turbulence in easterly wind. In addition, the simulation results for easterly wind are compared to those for northerly wind, which was, in the earlier in-situ data analyses, associated with values of DEL and standard deviations of the wind speed that were smaller than those in easterly wind (hub-height wind speeds of approximately 9 m/s (Section 3)).

Figure 12 shows the distribution of the instantaneous streamwise wind velocity component in the vicinity of WT #10. This figure provides a visual indication that, when easterly flow is present, separated flow (terrain-induced turbulence) forms over Mt. Benzaiten, which is located upstream of WT #10, and that this separated flow strongly affects the wind turbine.

Figure 13 illustrates a vertical profile of the simulated wind velocity vectors (instantaneous flow field) at the site of WT #10 in easterly wind. Figure 14 shows

Figure 12. Distribution of the streamwise wind velocity component on a vertical cross-section which includes WT #10, instantaneous flow field, easterly wind.

Figure 13. Wind velocity vectors (instantaneous flow field) at WT #10, easterly wind.

Figure 14. Vertical profiles of the streamwise wind velocity and related variables at WT #10, easterly wind.

the vertical profile of the streamwise wind velocity (instantaneous flow field) for the same time as Figure 13. Figure 14 shows the vertical profiles of the streamwise wind velocity averaged over 10 min (actual time) (purple line) and the range of the instantaneous streamwise wind velocity (blue line). A close examination of this figure makes it clear that the magnitude of wind velocity fluctuations is large within the swept area of WT #10 in the case of easterly wind.

Figure 15 shows the time series data of the simulated streamwise (x) wind velocity component (u) at the wind turbine hub height (60 m above the terrain surface) for 10 min in actual time. In the easterly wind case shown in Figure 15(a), the fluctuation of the streamwise velocity component is larger than that in the northerly wind case shown in Figure 15(b). Further analysis (not shown) also reveals that there exists a periodicity of about 6 to 7 seconds in actual time in the streamwise wind velocity fluctuation in the easterly wind case (Figure 15(a)).

On the other hand, in the case of northerly wind, in which the value of DEL was 0.95 in the earlier in-situ data analyses (Section 3), the magnitude of the fluctuation of the streamwise wind velocity is smaller than that in the case of easterly wind, yielding 0.99 m/s for the standard deviation of the streamwise wind velocity (Figure 15(b)). Here, the present study proposes a planning method for optimal wind turbine deployment which utilizes numerical simulation results. For this purpose, an original index (= σu (the standard deviation of the wind turbine hub-height streamwise wind velocity)/U (inlet boundary streamwise wind velocity = 10 m/s)) called the turbulence evaluation index is defined. Table 4 shows the values of the turbulence evaluation index calculated with the

Figure 15. Time-series data of streamwise wind velocity from the numerical simulations: (a) Easterly wind and (b) Northerly wind.

results from Figure 15 and also shows the in-situ data that were collected over the wind farm. This table shows a strong correlation between the values of DEL that were evaluated from the in-situ data and the values of the turbulence evaluation index that were calculated from the simulation results. This result is consistent to the strong correlation between the values of DEL and the standard deviation of the wind speed that was found from the in-situ data analyses in Section 4. Finally, it is speculated that terrain-induced turbulence which results from complex terrain and is shown in Figures 12-15 causes metal fatigue of the wind turbine components to accumulate at a faster rate than originally estimated.

8. Conclusions

In the present study, in-situ data analyses and high-resolution numerical wind simulations that are based on large-eddy simulation (LES) were performed for Wind Turbine (WT) #10 at the Kushikino Reimei Wind Farm (in operation since Nov. 2012), which is owned by Kyudenko New Energy Co., Ltd. and located in Hatori District, Ichikikushikino City, Kagoshima Prefecture, Japan. The results of the analyses and simulations revealed that in easterly wind, WT #10 was directly affected by terrain-induced turbulence which originated at and was generated by Mt. Benzaiten (elevation: 519 m), a mountain upstream of (to the east of) WT #10.

In the hours selected for detailed in-situ data analysis, in which the mean hub-height wind speed at WT#10 was approximately 9 m/s, the value of DEL (blade bending) for WT #10 in easterly wind was approximately twice as large as that in northerly wind. In addition, the values of the newly-proposed turbulence evaluation index that were evaluated from the numerical wind simulations for easterly and northerly wind were 0.236 and 0.099, respectively (Table 4). The value of this index in easterly wind is approximately twice that in northerly wind, thus, the ratio of the value of the index in easterly wind to that in northerly wind was similar to the ratio of DEL from the same wind directions in the in-situ analysis. These results revealed that the values of DEL (blade bending) that were based on the in-situ blade strain data and the turbulence evaluation index that were calculated from the numerical simulation results were positively correlated. Furthermore, the results from above confirm for the first time that the use of 0.2 or less for the value of the turbulence evaluation index, which our research group has frequently adopted as a guideline for wind turbine deployment to reduce DEL, is appropriate.

The series of qualitative and quantitative examinations in the present study has revealed that wind turbine deployment planning which takes into account the durability of blades (blade bending) is possible with the use of the turbulence evaluation index, which is output from the RIAM-COMPACT software package at the planning stage of wind turbine deployment. In the future, in order to assess and predict the fatigue strength and lifetime of wind turbine components, our research group will attempt to identify the relationship between wind turbine operation data and three-dimensional wind data which are output from RIAM-COMPACT and also establish an advanced method for analyzing wind turbine failures.

Acknowledgements

For conducting this research, the authors were provided with in-situ data by Kyudenko New Energy Co., Ltd. For the analyses of the in-situ data, the authors received cooperation from E-WIND, Co., Ltd. Furthermore, the present research was supported by joint research between Kyushu University and West Japan Engineering Consultants, Inc. (Research title: Establishment of a method for optimal wind turbine deployment which uses numerical simulation results and takes wind turbine structural strength criteria into consideration; Research period: Oct. 13, 2015 - Mar. 19, 2016; Principle investigator: Takanori Uchida) and joint research between Kyushu University and Hitachi, Ltd. (Research title: Joint research on the effect of terrain-induced turbulence on wind turbine structural strength; Research period: Jan. 8, 2016 - Mar. 31, 2016; Principle Investigator: Takanori Uchida). The authors wish to express their gratitude to all the organizations and individuals who have contributed to the present research.

Cite this paper

Kawashima, Y. and Uchida, T. (2017) Effects of Terrain-Induced Turbulence on Wind Turbine Blade Fatigue Loads. Energy and Power Engineering, 9, 843-857. https://doi.org/10.4236/epe.2017.913053

References

  1. 1. Uchida, T. (2017) CFD Prediction of the Airflow at a Large-scale Wind Farm above a Steep, Three-dimensional Escarpment. Energy and Power Engineering, 9, 829-842. https://doi.org/10.4236/epe.2017.913052

  2. 2. Uchida, T. (2017) Large-Eddy Simulation and Wind Tunnel Experiment of Airflow over Bolund Hill. Open Journal of Fluid Dynamics, Accepted.

  3. 3. Uchida, T. and Ohya, Y. (2011) Latest Developments in Numerical Wind Synopsis Prediction Using the RIAM-COMPACT® CFD Model Design Wind Speed Evaluation and Wind Risk (Terrain-Induced Turbulence) Diagnostics in Japan. Energies, 4, 458-474. https://doi.org/10.3390/en4030458

  4. 4. Kiyoki, S., Uchida, T., Kawashima, Y., Kondo, K. and Nishida, T. (2016) A Research of the Effects of Terrain-Induced Turbulence on the Wind Turbine Fatigue Load (Part 3: Actual Measurement Evaluation and Prediction Technoloty Development for Fatigue Load in Complex Terrain). Japan Wind Energy Association, In Proceedings of the Japan Wind Energy Symposium, 38, 459-462.

  5. 5. (2014) IEC 61400-1 Ed 3.1: Wind Turbines - Part1: Design Requirements. International Electrotechnical Commission, ISBN: 978-2-8322-1525-8.

  6. 6. (2010) Wind Power Generation Equipment Support Structure Design Guideline. Japan Society of Civil Engineers, ISBN: 978-4-8106-0705-5. (in Japanese)

  7. 7. Kato, M. (1994) Wind Tunnel Simulation of Atmospheric Turbulence over Complex Terrains. Japan Association for Wind Engineering Magazine, 59, 89-92.

  8. 8. Kim, J. and Moin, P. (1985) Application of a Fractional-Step Method to Incompressible Navier-Stokes Equations. Journal of Computational Physics, 59, 308-323. https://doi.org/10.1016/0021-9991(85)90148-2

  9. 9. Kajishima, T. and Taira, K. (2016) Computational Fluid Dynamics: Incompressible Turbulent Flows. Springer, ISBN: 978-3319453026.

  10. 10. Kawamura, T., Takami, H. and Kuwahara, K. (1986) Computation of High Reynolds Number Flow around a Circular Cylinder with Surface Roughness. Fluid Dynamics Research, 1, 145-162. https://doi.org/10.1016/0169-5983(86)90014-6

  11. 11. Smagorinsky, J. (1963) General Circulation Experiments with the Primitive Equations: I. The Basic Equations. Monthly Weather Review, 91, 99-164. https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2