Near-wake analysis of actuator line method immersed in turbulent flow using large-eddy simulations

The interaction between wind turbines through their wakes is an important aspect of the conception and operation of a wind farm. Wakes are characterized by an elevated turbulence level and a noticeable velocity deficit, which causes a decrease in energy output and fatigue on downstream turbines. In order to gain a better understanding of this phenomenon this work uses large-eddy simulations together with an actuator line model and different ambient turbulence imposed as boundary conditions. This is achieved by using the Simulator fOr Wind Farm Applications (SOWFA) framework from the National Renewable Energy Laboratory (NREL) (USA), which is first validated against another popular Computational Fluid Dynamics (CFD) framework for wind energy, EllipSys3D, and then verified against the experimental results from the Model Experiment in Controlled Conditions (MEXICO) and New Model Experiment in Controlled Conditions (NEW MEXICO) wind tunnel experiments. By using the predicted torque as a global indicator, the optimal width of the distribution kernel for the actuator line is determined for different grid resolutions. Then, the rotor is immersed in homogeneous isotropic turbulence and a shear layer turbulence with different turbulence intensities, allowing us to determine how far downstream the effect of the distinct blades is discernible. This can be used as an indicator of the extents of the near wake for different flow conditions.


Introduction
An important aspect for the conception of wind farms is the turbine spacing which depends on the interaction of wind turbines through their wakes. This phenomenon can decrease the wind park energy output by up to 20% due to the velocity deficit propagated by the wakes (Manwell et al., 2010). Additionally, it can increase the turbine fatigue due to the increased turbulence 15 intensity. In order to study wake interactions, the flow around the rotor has to be modelled correctly. Hence the model should account for the apparition of turbulent structures of different magnitudes. For instance the vortices created by the blade tips and its interaction with the ambient turbulence.
As opposed to the far wake region (Olivares Espinosa, 2017), the near wake representation in a computational fluid dynamics simulation depends heavily on the applied rotor model. Approaches range from an actuator force representation inserted as 20 momentum sink in the Navier-Stokes equations to full rotor modelling where the attached boundary layers on the blades are simulated (Sanderse et al., 2011). This work will apply the actuator line method (ALM) in order to model the transient behaviour of the rotor by representing distinctly the rotating blades as presented by Troldborg (2009). Each blade is represented by a force line allowing to reproduce the helicoidal vortical structure in the near wake allowing to to assess its interaction with the flow.
1 In order to evaluate the soundness of the present method a comparative study of the SOWFA framework, from NREL, and EllipSys3D, from DTU, was conducted as initally ::::::: initially presented in Nathan et al. (2017). Based on this study, the method used throughout this work will be evaluated before proceeding to establish the base case for the non-turbulent inflow.
For establishing the base case, the optimal width of the distribution kernel of the forces of the actuator line is determined.
While previous work often focused on numerical stability as in Troldborg (2009) or Ivanell et al. (2010), when choosing the 5 distribution width, Martínez-Tossas et al. (2015a) states that with decreasing distribution width, the line forces are getting too concentrated resulting in a wrong prediction of the rotor torque. Hence this work tries to evaluate the optimal width for each mesh resolution by using the predicted torque as a global indicator.
For the introduction of a turbulent inflow different methods exist for imposing a stastically ::::::::: statistically : generated velocity field, such as inserting it via a momentum sink as done in Troldborg et al. (2011) or as boundary conditions as done in 10 Olivares Espinosa (2017). This work adheres to the latter approach, as it was seen as more straightforward than the conversion of the velocity field to a force which then is translated back to the velocity field by the numerical solveras done the former approach.
Then shear layer turbulence (Muller et al., 2014) is introduced exposing the rotor model to a more realistic wind flow situation bearing more resemblance to applied wind energy. This novel approach takes into consideration the temporal evolution of the 15 sheared velocity field, hence allowing it to be imposed as boundary condition as well. Also this method was not yet applied to the actuator line method.
Finally, the numerical results are used to examine the spatial extents of the near-wake region. While in previous work such as Krogstad and Eriksen (2013) or Sarmast et al. (2016) often the profiles of velocity deficits or turbulent kinetic energy are taken into consideration for evaluating the near-wake, this work uses the energy spectra to determine how far downstream the 20 discernible effects of the distinct blades are noticable :::::::: noticeable. While the analysis of the turbulent inflow in previous work (Olivares Espinosa, 2017) has often been conducted using energy spectra, seldom energy spectra including the rotor effects are included in the analysis of the near-wake. As with increasing turbulence intensity the statistical convergence tends to be longer, the energy spectra approach in this work permits the analysis of the spatial extensions of the near-wake region even without fully convergence of second-order statistics. 25

Numerical methodology
The numerical simulations are based on the incompressible Navier-Stokes equations ∂U ∂t + ∇ · (UU) − ∇ · (ν∇U) = −∇p/ρ + F (1) with F representing the actuator force inserted as a momentum sink, U as the velocity, ∇p/ρ as the modified pressure and ν as an overview is given on the methods generating the two different ambient turbulences :::::::: turbulence : and how they are imposed as boundary conditions. Finallythere is : , a summary of the numerical framework and its setup : is ::::: given : in section 2.3.

Rotor model
The force term F in Eqn.
(1) is obtained using with the lift and drag forces shown in Figure (1) and defined as with c l and c d as the lift and drag coefficient, U mag the sampled velocity magnitude in the blade reference frame, c the chord width, l s the length of the actuator segment, Gaussian kernel G :: the :::::::: Gaussian :::::: kernel and tip correction f tip ::: f tip .

Homogeneous isotropic turbulence
A synthetic velocity field representing homogeneous isotropic turbulence based on the von-Kármán energy spectrum (Pope, 2000) 15 is obtained by using the algorithm proposed by Mann (1998). The technical details can be found in the article of Mann (1998) or more recently in Olivares Espinosa (2017). The main parameters for this approach are the integral length-scale L and the coefficient α 2/3 which can be used as a scaling factor to obtain the desired amplitude of the turbulent structures. The range of wavenumber :::: wave ::::::: number κ depends on the grid resolution and dimension extents. Hence, these parameters determine the 20 ability of the numerical mesh to resolve a certain range of turbulent scales.
In Figure (2) the midplane of a generated turbulent field is shown for different turbulence intensities. The flow structures are identical apart from the different scaling of the velocity fluctuations. This results from using the same seed for the random number generator in the Mann algorithm and by scaling the obtained velocity field with α 2/3 to obtain the desired synthetic 5 turbulence intensity T I syn .
Contrary to Troldborg (2009) where the synthetic turbulence is imposed as a momentum source, in this work it is imposed as a boundary condition as done in Muller et al. (2014) and Olivares Espinosa (2017). The velocities are imposed by convecting the velocity field of the synthetic turbulence through the computational domain by the mean velocity U ∞ at each simulation time step. They are then projected by trilinear interpolation onto the computational points. In order to speed up the statistical 10 convergence, the simulation is also initialized with the synthetic turbulence field.

Numerical framework
This work is realized within the open-source framework OpenFOAM 1 (version 2.2.2) together with the SOWFA 2 project, which contains a similar implementation of the ALM as presented by Troldborg (2009). A more detailed explanation of the 5 implementation can be found in Martínez-Tossas et al. (2016). OpenFOAM is a set of libraries and executables entirely written in C++. While the first released scientific article about the framework was by Weller et al. (1998) its inner workings are described more in-depth by Jasak (1996).
The large eddy simulations use the dynamic Lagrangian sub-grid scale model (Meneveau et al., 1996). For the discretization of the convective term a linear combination of 75% central differencing and 25% of a second-order upwind scheme is applied as presented by Warming and Beam (1976). In OpenFOAM terminology this scheme is called "Linear-upwind stabilized 5 transport" (LUST). The choice of the scheme is made as a trade-off between the accuracy of a linear discretization and the stability of an up-winding scheme. This scheme proved to preserve well the turbulent structures (Nathan, 2018). The remaining spatial terms are discretized by central differencing and for the time discretization the Crank-Nicolson method is used.
The pressure is resolved using a geometric agglomerated algebraic multi-grid solver and the remaining variables are solved for with a bi-conjugate gradient method using a diagonal-based incomplete LU preconditioner. The total simulation run-time 10 comprises 60 rotor revolutions ::::::: (≈ 8.5 s) : and the time step has to be small enough to avoid the actuator point representing the blade tip skipping a computational cell during rotation. It is also an integer fraction of the rotor revolution time. This ::: set :: at :::::::::::: 0.327 · 10 −3 s. ::: The :::: total : run-time is chosen as first and second order statistics are deemed to be converged.
For the parametrization of the ALM, different distribution widths are chosen in order to obtain the optimum for the examined case and 40 actuator points are used to represent one blade in accordance with what was found in Nathan (2018). 15 3 Results

Non-turbulent flow
When refining the grid using the actuator line method the distribution parameter has to be adjusted to obtain a global torque T close to the reference value T ref . In the following only the case for U ∞ = 15 m/s will be examined as the other cases in Nathan et al. (2017) served as extreme cases for determining how the model behaves at its limits.
Instead of relying on a constant /∆x for different grid resolutions, this work adapts /∆x depending on the grid resolution, 5 or number of cells across the rotor diameter N = 2R/∆x. This can be seen as a first step towards an actuator surface method, where the force is distributed with respect to the blades chord. The results are shown in Figure (6). A confidence interval of ±1 % was established around the reference torque value T ref . Through iterations, an optimal distribution parameter is found to fall in this range.
The lower bound for the distribution parameter here is = 1.7∆x for the sake of numerical stability of the here chosen  (2010). By doing so it can be seen that the best solution in terms of global torque for a resolution of N = D/∆x = 32 is off by around 4 % in Figure (6).
As a general trend it can be seen that /∆x has to be increased with increasing resolution. This stems from the fact that by refining the mesh with a constant /∆x the punctual induction caused by the blade would be too high and eventually the torque would be below the reference value, e.g. ::: e.g. for = 2∆x for N ≥ 64. On the contrary, when having a very low resolution 5 a constant distributes the force too widely, causing a lower induction around the rotor resulting in an overestimation of the torque, e.g. ::: e.g. for = 2∆x for N ≤ 48.
The optimal distribution parameter /∆x found in Figure (6)  It can be seen that the method seems to converge towards a solution when refining the mesh. As shown in Figure (6) the lowest resolution at N = 32 over-predicts the torque by distributing the force to widely which also reflects in the low axial induction downstream at x/R = 0.13. Despite following well the trend of the experimental values the method seems to converge towards radial profiles which are especially off in the tip and hub region where the strongest vortices are shed.

Homogeneous isotropic turbulence
In Figure (11) the longitudinal evolution of the turbulence intensities can be seen. There is a stronger decay for higher turbulence intensities which was also found in Olivares Espinosa (2017). In that work EllipSys3D was compared to a solution based on OpenFOAM and it was found that over the same longitudinal distance of 10R a decay :: an ::::::: absolute ::::::::: difference :: in ::: the ::::::::: turbulence ::::::: intensity of 48 % and 44 % occurred for each framework respectively. This stands in a stark contrast to the 4 % in this case for 5 the high turbulence intensity case. This huge decay, which is even more significant for EllipSys3D, necessitates to approach the introduction of the turbulence close to the turbine for high turbulence intensity cases (Olivares Espinosa, 2017).
An important aspect when imposing a synthetic turbulence as boundary conditions of a CFD simulation is respecting the Nyquist-Shannon sampling theorem (Shannon, 1949) as also mentioned by Muller et al. (2014). Hence a study considering different ratios between the grid resolution of the synthetic turbulence and the simulation was undertaken. It is found that the 10 higher the computational resolution is compared to the one of the synthetic turbulence, the less the turbulence intensity decays in longitudinal direction. While the criterion of Nyquist-Shannon states that the resolution of the computational domain should be at least twice as big ::::::::: dx/∆x > 2, this work uses the ratio of dx/∆x = 2.5 with dx as the cell width of the synthetical ::::::: synthetic field and ∆x as the cell width of the computational mesh.
It is interesting to notice the distinct peaks in the spectra occur at the wavenumber relating to the frequency of the blade passage and its harmonics.

5
Although flow properties were maintained at a similar level, the sheared flow has a clear impact on the power extraction and also on vortex properties of the structures emitted by the blade. A very interesting observation is the fact that for higher turbulence intensities the effects of the distinct blades using the ALM seems to vanish at relatively short downstream distances.
This poses the question of the usability of the ALM when arguing for its capabilities of representing the transient behaviour and its impact on downstream turbines.

Conclusions
By using a validated actuator line implementation (Nathan et al., 2017), it was shown that the distribution width /∆x has a non-linear dependence on the grid resolution and converges probably towards values suggested in Martínez-Tossas et al.
(2016). The rotor torque is used as a global indicator for determining the distribution width, but the rotor thrust followed the same trend. Hence it is interesting to see that while the rotor induction is predicted well, the velocity deficit agrees well only for x/R > 5 but not in the ultimate rotor vicinity.
It is also shown that with increasing grid resolution the spatial profiles seem to converge. This would be one aspect of a grid independent solution, but it is still very far away from resolving correctly the shed tip vortices. Although it seems to converge 5 towards a value of /∆x ≈ 4−5 for which the dimensions of the experimental vortices would be attained, this causes excessive computational costs due to the large mesh.
When looking at the turbulent inflow, : a ::::::: synthetic ::::::::: turbulence ::::::::: generated :: by ::: the ::::: Mann ::::::::: algorithm ::::::::::: (Mann, 1998) : , it was shown that the decay of the turbulence intensity in longitudinal direction is much less pronounced than in previous work. As shown for the axial decay of the turbulence intensity a significant part of the difference between the resolved turbulence intensity and 10 the imposed one from the synthetic field, resides within the sub-grid scales. Hence there is very little loss due to numerical dissipation which also reflects in the energy spectra which are the better the higher the turbulent content is.
As expected the wake does recover at a faster pace for a higher turbulence intensity. It is very interesting to notice that the turbulent structures of the ambient flow eventually catch up with the amplitude of the structures emitted by the rotor. This is already noticable ::::::::: noticeable in the instantaneous velocity fields but becomes even clearer when evaluating the spectra. When can be observed that in this case for T I syn ≥ 10 % the near wake already ends at x/R = 0.4. This is particular interesting as a turbulence intensity of 15 % at hub height is still considered to be low turbulence intensity according to ISO 61400 and many real sites exhibit even higher turbulence intensities. Hence for some cases the limit of the near wake would be x/R = 0.4 and even lower.

20
Code and data availability. The SOWFA framework on which this work is based is made available by NREL https://github.com/NREL/ SOWFA/ and the turbulence generator for the homogeneous isotropic turbulence can be obtained via http://vbn.aau.dk/en/publications/ tugen(3e097a90-b3d8-11de-a179-000ea68e967b).html. The results for the NEW MEXICO experiments were provided upon request by Gerard Schepers.
Competing interests. Christian Masson is a member of the editorial board of the journal.

General remarks:
This paper is of overall very good quality. A detailed discussion on the impact of the actuator line parameters is given. Even so no generic solutions are given, some metrics can be extracted from the work, regarding the optimal values of the Gaussian width parameter. Based on this validated actuatorline model, the authors propose a study on the impact of a homogeneous, isotropic turbulent inlet and shear layer turbulence. While the impact of the homogeneous isotropic turbulence is clear and well exposed thanks to the provided spectra, the additional impact of the shear layer is less obvious and should be discussed in more details. Done, see P3L11.
-P3 L12. It is disappointing to see that purely 2D airfoil polars are used, while 3D effects are discussed.
It was shown in Nathan 2017, that it is more than sufficient for the 15 m/s case.
-P3 L16 → P4 L3. The authors argue the Glauert tip correction should be used due to the low resolution. According to Churchfield et al. (2017), this is due to the isotropic kernel that is used, and the virtual projection of forces outside of the "blade domain". It could be interesting to see which phenomena is dominating, i.e. the lower resolution or the isotropic projection. Furthermore, I was not able to find reference (Nathan, 2018). Is it already published? Otherwise, please mention it in the references.
It would definitely be interesting, but it would be the matter of future work. Nathan, 2018 references to my PhD thesis, which should already be available since a while in a digital form at the library of ETS. I am in contact with the school in order to find out, why this has yet happened. I handed it in more than 6 months ago.
-P6 L7. The domain is rather small in length compared with standard recommendations. As a comparison, N. Troldborg (2009) uses a domain length of 18R, MartinezTossas (2015) use a domain length of 21D. . . It could be useful to provide some proof of convergence.
A complete sensitivity study was conducted in my PhD thesis. Including the numerous graphics and tables would blow up this article unnecessarily.
-P6 L25. If possible, provide some orders of magnitude for the time step.
-P7 L6. The under-prediction of the axial induction is not clear to me; results are almost super-imposed to the NewMexico measurements. C2 True, I rephrased that passage.
-P8 L5 → L7. I do not understand the link with the actuator surface method. Even so epsilon over dx is adapted, depending on the Cartesian grid refinement, it is still an actuator line, and no chord-wise meshing is used.
I agree, it is still an ALM, passage adapted.
-P8 L14 → 18. Discussion regarding the impact of the epsilon parameter is very instructive. Giving a look at the results, it seems to me there is an almost linear relation between the optimal epsilon value (leading to T/T_{ref} = 1) and the mesh refinement parameter N in the range 50 < N < 150. It could be interesting to derive an empirical law from it. In case the results presented in Figure 6 are not "rotor specific", this could lead to a very simple law to derive the value of epsilon "on the fly".
I share your enthusiasm on this. Probably this method should also be applied to e.g. NTNU blind comparison test or other experiments and see whether an empiricial relation could be derived. In my opinion it is slightly non-linear as seen in Fig.7.
-P9 L13→L16. The "bad" resolution near the tip is, according to the authors, attributed to the actuatorline representation of the blades. In Blondel et al. 2017, a lifting-line model is used together with a vortex model of the wake, and better correlation with experimental data was obtained (compared to SOWFA simulations). Thus, from my point of view, discrepancies should be attributed to the isotropic kernel in use, which is unrealistic near the tip, or the potential excessive diffusion of the finite volume scheme. The effect of the isotropic Gaussian kernel and the mesh effects, as discussed P10., should be further analyzed (not necessarily in this publication, as this is not the main topic).
I agree, I think it goes in hand with your point mentioned above and a possibility for future work.
-P12. In this part, synthetic turbulence is imposed at the inlet of the domain. I guess results presented in Figure 11 are based on simulations without the wind turbines, considering only the evolution of the TI in a channel. A net decrease of the turbulence intensity is observed along the channel. The simulations are not really described here. It could be interesting to provide some proof of convergence. Is the simulation time long enough to transport the characteristics defined at the inlet? This could be a basic explanation to the decrease of the TI that is observed. Also, it seems something is happening at the outlet, with some kind of sharp TI recovery. If this belongs to ghost cells, it should not be present in this plot. The effect of the ratio between the CFD grid and the SEM grid is discussed. One can wonder about the turbulent length scales used C3 in the SEM algorithm, and their relation with the size of the computational grid. Is the CFD grid fine enough to catch the vortices given at the inlet?
Yes, the numerical grid cell size was chosen to respect the Nyquist-Shannon criterium, hence it was more than twice as refined as the synthetical grid.
-P14. Figure 15 is not useful and could be removed. Same remark holds for Figure 17, results are similar to the homogeneous isotropic turbulence case.
-P15. In 2.2.2., the evolution of mean velocity with height is presented. However, is seems to me, based on the contours of vorticity, that the TI is constant with height. Is that correct? From a physical point of view, higher TI is expected near the ground. This should lead to higher vorticity. Can the author provide some insight? In figure, y label should be z/R.
Passage adapted, TI changes with height, but you are right, it is rather hard to notice with the vorticity plots.
-P16. L7. I do not understand the point here. Please reformulate.
-P16. L11. Are the authors talking about the global rotor power? No metrics are given. The impact on vortex structures is not clear to me. Differences in the spectra between figure 19 and 14 are rather small. I would suggest including an additional metrics here to clarify the impact of the shear on the near wake.
Passage removed. It is discussed more in detail in my thesis, but it does not necessarily has to be in the article.
-P17. Conclusions. Based on the observation of the turbulence decay ( fig. 11) with axial position, it seems that at high TI, a large part of the TI is included in the subgrid scales. Therefore, I am not totally comfortable with the conclusions that are given: the length of the near-wake is determined based on the observation of the vorticity. However, the subgrid-scales are not included in the vorticity. Therefore, it seems difficult to draw definitive conclusions. As a more general remark, this work emphasize on the impact of the TI on the near-wake. However, blade loads may also be impacted, even at the airfoil polar level. One might expect a delay in the stall at high TI, which could impact the blade root loads. Also, it could have been interesting to use the NTNU Blind Comparison experiments as a complementary validation case.
True, I think it is a shortcoming of my work, that I did not have the time to conduct this experiment with another rotor. It is true, that some of the turbulent motion falls within the subgrid, but as shown in my PhD thesis, more than 90% fall within the resolved scales. Hence I only take the resolved scales into consideration when looking at the limits of the near wake. Done.
-P3 L8 Gaussian Kernel G → G the Gaussian Kernel Done.

Anonymous Referee #2
The manuscript entitled "Near wake analysis of actuator line method immersed in turbulent flow using large-eddy simulations" deals with numerical computations of a modelled wind turbine embedded in LES computations. Several types of inflow turbulence are implemented in order to study its effect on the near wake development. The spectral content of the wake flow is then used as indicator of modifications. The present study is of interest for the wind energy community ad gives valuable insights on the near wake extent depending on turbulence intensities and types. On the other hand, the manuscript needs several improvements before to be accepted for publications.
-In general, the figure captions are poor and lack of essential information. Added, see P7L10-15.
-What are the properties of the generated turbulence in term s of scales: integral length scale, ratio between this scale and the rotor radius? Since the authors present some turbulence spectra without wind turbines, it would be worth to develop a better description of the generated turbulence.
Added length scale of incoming turbulence P5L9, R/4.
-The PhD thesis of the 1st author is cited regularly, whereas the reference is not precise enough to ensure that the reader will find it easily on the web. Additionally, this reference is sometimes cited for results which are not specifically a new outcome of this thesis. So please cite more relevant sources when possible.
Nathan, 2018 references to my PhD thesis, which should already be available since a while in a digital form at the library of ETS. I am in contact with the school in order to find out, why this has yet happened. I handed it in more than 6 months ago.

Major comments:
-P2, lines 22-24: the spectra are also based on statistics. I do not see why it would be less sensitive than second-order statistics to the convergence issue, since the spectra is a frequency distribution of the variance.
In wavenumber space it is easier to detect at which point the structures emitted from the blade merge with the surrounding turbulence. Also time series for energy spectra were taken at different points at the circonference of the rotor. Hence the spatial aspect can be neglected and it the energy spectra for each point can be averaged and the result is a smoother energy spectrum. This would not be possible with a radial profile of TI due to the asymmetric nature of the helicoidal vortical structures.
-P3, line 10-12 : this data was obtained from wind tunnel experiments and therefore it does not include the stall delay due to boundary layer stabilizing effects such a Coriolis and centrifugal forcing which enhance the lift of the airfoil". These types of effects can be reproduced in a wind tunnel. Do you mean 2D experiment? Without rotation?
Yes, passage was adapted.
-Page 5, lines 4-6. Give some details about the synthetic turbulence fields: what are the turbulent length scales ? it seems that the turbulence does not dissipate with the axial distance. That is not so common, as explained later on in the manuscript. Some comments should be already mentioned here.
Added, P5L9. It does dissipate, see Fig. 11, but at a lower pace than in other work.
- Figure 5 : improve the caption and give details about the experimental configuration used as reference here : power coefficient, thrust coefficient, TSR, etc Done, P7L10-15.
-Page 10, line 15-17 : this part is confusing : it seems that a relative decrease of turbulence is given (48% and 44%), whereas the 4% stands for a decrease of turbulence from 15% to 11%. This would correspond to a relative decrease of 25%... please rephrase this part.
-Page 16, lines 11-15. Please elaborate more on the discrepancies between sheared and un-sheared conditions Removed passage.
-Conclusion: there is not discussion about the relative size of the inflow turbulent structures compared to the wake turbulent structures (rotor size, blade size, tip vortex size, shear layer size?). It is indeed a very important parameter to justify the observations mentioned in page 17, lines 16-23. Minor comments: Done, P7L10-15.
- Figure 2 : if it the midplane at y/R = 0, the plot should be dependant of z/R and x/R Done.
-P4, line 3: please give the definition of sigma and Delta x. Done.
- Figure 8 is too small. Additionally, it is difficult to differentiate the experimental and numerical results Done.
-Page 10, line 2: The Reynolds number is based on the circulation: Please explain why you use this definition and not another one.
- Figure 10 : the caption is wrong Done.
- Figure 11: the caption is poor and the authors should also better explain in the body text what this figure is for. Which computation solver is used here? Done.
- Figure 18 : Y axis is not consistent with the caption Done.
-Conclusion: remind in the conclusion the used method to generate the turbulent inflow Done.