LES : Unsteady Atmospheric Turbulent Layer Inlet . A Precursor Method Application and Its Quality Check †

The motivation of this work is to bridge the gap between experimental approaches in wind tunnel testing and numerical computations, in the field of structural design against strong winds. This paper focuses on the generation of an unsteady flow field, representative of a natural wind field, but still compatible with Computational Fluid Dynamics inlet requirements. A simple and “naive” procedure is explained, and the results are in good agreement with some international standards.


Introduction
Inlet conditions for Large Eddy Simulations (LES) are of the utmost importance, as it determines a large part of the fluid behavior within the computational domain.Tabor and Baba-Ahmadi [1] gave a good overview of the different techniques, classifying them as synthesized turbulence methods and precursor simulation methods.
The first category presents the critical limitation of not achieving a flow field inlet 100% compatible with Computational Fluid Dynamics (CFD) requirements (temporal and spatial fluctuations, divergence-free, spectrum).For instance, the method developed by Jin et al. [2] is not divergence-free.

OPEN ACCESS
These requirements are of great importance when the purpose of the computations is to generate pressure fields on obstacles.Some regularization methods have to be introduced, as for instance the Synthetic Eddy Method [3].
The second category of techniques is due to Spalart and Leonard [4].The idea of these methods is to generate inflow condition using CFD itself, as a genuine simulation of turbulence.The inflow of this precursor domain is obtained through a rescaling of the flow field extracted from a downstream location.Nevertheless, even with different stages of simplification [5,6], those techniques are still complicated and not easy to implement.Additionally, they seem more adequate to small turbulence level (Reynolds number Re about 2 × 10 3 ), which does not correspond to structural design against strong winds (Re about 4 × 10 7 ).
This paper deals with a more "naive" precursor model.The idea is to draw inspiration from the wind tunnel testing community, with the use of roughness blocks lying on the floor.The large x-axis dimension of the wind tunnel is replaced by a short cyclic domain, in order to convert the large domain issues into a duration matter.Similarly to atmospheric turbulent layer which is driven by geostrophic wind (instead of a horizontal pressure gradient), our flow is naturally driven by the upper boundary condition, acting as a conveyor belt.Re-introducing a flow field extracted downstream becomes straightforward.
The final aim of the method being building dimensioning to high winds, our reference will be an international code concerning these issues: Eurocode I (EN-1991-1-4): "Actions on structures-Wind actions" [7].This code proposes a systematic approach to describe some of the wind characteristics; we will focus on these characteristics in order to validate-or not-the approach.

Aim of the Simulation
In order to compute the extreme loads on a building, CFD calculations have to immerse it in a turbulent atmospheric air flow.Eurocode I (EN-1991-1-4) proposes a systematic and simple approach to describe this air flow: a log-law profile for the mean wind speed, associated with an inverse log-law profile for turbulence intensity.
where is the 10-min averaged wind speed, is the height above ground, is the ground roughness length, is the regional base velocity, and is a parameter depending only on .% is the turbulence intensity, σ is the standard deviation of the instant velocity (computed on a 10-min window), and is a parameter depending only on .
Wind speeds distributions are Gaussian and the spectrum of turbulence is driven by a modified Kaimal model: where is the single sided power density function,  is the frequency, and is the reduced frequency: The turbulence scale is also given by Eurocode: There is no mention of correlation lengths within the Eurocode.Another regulation, the English code ESDU 75001 [8], gives the correlation lengths in the streamwise direction: where , and are the streamwise correlation lengths, respectively for the components (stream wise), (lateral) and (vertical).
Our target is to generate an unsteady flow corresponding to this wind model, with a roughness length equal to = 5 cm and a base velocity equal to = 32 m/s.In this case, we have the fallowing values: At 5 m high: (5) = 44 m,

Originality of the Method
The main idea of the method consists of two parts: explicit roughness-cubes on the ground, and short cyclic domain.
In true life, roughness is a succession of obstacles that "sculpt" the air flow and add turbulence.In classic RANS computations, these can be modeled by a statistical approach using some wall functions that describe the mean profiles close to the rough surface.This generates a smooth and uniform flow, with a profile corresponding to the roughness length.However, for LES inlet, it is necessary to describes explicitly the changes in time and space of the wind field: this needs horizontal and vertical strong gradients, which can be generated by some blocks on the floor.The wakes of each element will act as a force opposite to the air flow, and a local turbulence generator.Large scale turbulence will be then explicitly resolved.
A fully developed turbulent layer takes a long span to establish (width δ is (√ ) and also depends on the roughness length).Thus, to obtain about 20 m-high boundary layer, the domain needs to be several hundred meters streamwise.This would result in a very heavy computation.Nevertheless, from a geometrical point of view, the domain would be similar to it-self; the only difference between a virtual segment of 10 m with its neighbor, would be the inlet.Thus, instead of computing a very long domain, a short domain is preferred with cyclic boundary conditions.However, the segment has to be long enough to integrate several rows of ground cubes.Finally, the initial large domain issues have been converted into a duration matter: the beginning of the simulation will be dedicated to the development of the turbulent boundary layer.

Geometry and Boundary Conditions
The domain used is a 10 × 20 × 20 m 3 box.On the ground, the roughness is modeled by random cubic elements: their size and density are given by some "rule of the thumb" from the experimental wind-tunnel community [9,10].We start from several cubes sided = 20 (here = 1 m), with an areal density around 10% (the gap between two consecutive elements is 2 ).Then some randomness is included: the eight of each block is randomized ( between 0 and 30 ), as well as its position (each block is allowed to move on its 3 3 individual panel).Special attention is focused on the lateral homogeneity (the largest cubes should not be grouped on the left or on the right).Figure 1 shows the roughness cubes at the bottom of the computational domain for this study.No slip Boundary Conditions (BC) are applied on the ground and on the roughness-cubes.Because the cubes are randomly located, a lateral force can appear on the ground.In order to compensate this force and to guide the flow within the x-axis, the lateral boundaries have to block the flow: symmetry boundary condition is chosen on these boundaries.Concerning the inlet and outlet, cyclic BC are used to re-introduce the flow field as it is.The upper-boundary condition is also no-slip on the "roof", but with a driven wall at 67.5 m/s (67.5 m/s instead of 0 m/s on the ground).This value of 67.5 m/s has been obtained through a try-and-error process: a first calculation, using 40 m/s on this Boundary Condition, gave a base velocity (evaluated as in the Section 3.3) equal to 19 m/s.As the flow is Reynolds quite independent, the results of the simulation may be proportional to this moving roof speed.The final calculation used the roof velocity of 40 × 32 ÷ 19 = 67.5 m/s.
Figure 2 summarizes these boundary conditions.

Running OpenFOAM
The mesh is a uniform grid (as recommended by A. Passalacqua [11]), with dimensions of 25 cm × 25 cm × 25 cm.In order to slightly reduce the cells number, a pinch of grading is introduced in the upper half space (last upper cell is 50 cm tall, see Figure 3a).This paragraph describes the parametrization used to operate OpenFOAM [12].The aim is to make people able to reproduce the computation.Pimplefoam is run in order to optimize the time step: maxCo (limit of CFL number) is 10, with a time step close to 4 × 10 −2 s.The time scheme is Crank-Nicholson, the spatial scheme is Gauss Linear (and Gauss LinearUpWindV GradU).The LES is performed with implicit filtering (cubeRootVol with deltaCoeff=1), and subgrid turbulence model is oneEqEddy.The solver for pressure is GAMG, and PBiCG DILU for U and k.The initial condition is to impose a uniform field U equal to 67.5 m/s.Computation has been parallelized on six CPUs using scotch decomposition method.Computation longs about 52 s per time step, which is equivalent to 23 min for one simulated second.The total duration for 600 s is about 10 days.

Spatial and Temporal Extraction Method
Many probes are placed onto an x-constant plane (1 m close to the outlet), in order to evaluate the success of the method to generate an atmospheric boundary layer: three vertical lines ( = 0, center of the domain, and = ± 5 m, each quarter of domain), to estimate the flow homogeneity and the vertical profiles, and a uniform horizontal line ( = 10 m, middle of domain) to assess the horizontal homogeneity and correlations (see Figure 3b).
A time span of 300 s is first run.During this first period, the velocity of flow rapidly decreases (from uniform 67.5 m/s), starting to brake from the ground.300 s turns out to be a good time span in order to achieve a global equilibrium (stationarity).Then the statistics are evaluated over a 300 s window, which is half the meteorological international standard [13].

First Check
Visual inspection (Figure 4) shows a first validation: the flow field looks like one could expect (heterogeneous, globally the higher the faster, and with some large "structures".There are no small structures because of the LES filter).

Results-Comparison with Eurocode
The flow is quite homogeneous: at 10 m high, the mean wind speed and turbulence intensity do not depend on the lateral position (see Figure 5a), at least in the central part of the flow (one half centered span).At 5 m high, densities are found symmetric, close to Gaussian (Figure 6, bottom).At 16 m high, this is not the case, especially for vertical component ( , purple points).This might be linked with the proximity of the "roof"; this has to be confirmed in a future work.Vertical profiles (mean wind speed, Figure 7, and turbulence intensity, Figure 8) agree very well with the Eurocode model.In those figures, the results from LES are compared with the Eurocode profiles (Equations ( 1) and ( 2)), with a roughness length of to = 5 cm and a based velocity equal to = 31.6m/s.This based velocity is found to be directly proportional to the top boundary condition speed.
Once again, close to the roof (17 m and above), results suffer discrepancy with the free atmospheric boundary layer modeled by Eurocode [7].The validity of our turbulent wind field is limited to the three fourth of the vertical domain (3/4 H), from the ground.3)).The peak of energy is slightly shifted to low frequencies, corresponding to a turbulence length scale slightly larger than the one proposed by Eurocode (Equation ( 4)).Nevertheless, the global balance between high and low frequencies is not accurately reproduced: low frequencies are over estimated whereas high frequencies are underestimated.This behavior might be linked with the LES filtering of high frequencies; this has to be confirmed in a future work.
Close to the ground (5 m high), there is a second peak of energy at a reduced frequency around = 4.This is the "signature" of the domain size: 10 m long, about 28 m/s, makes a cycling frequency equal to 2.8 Hz, corresponding to a reduced frequency (Equations ( 3) and ( 4 This peak of energy is also visible at 10 m high, at a slightly higher frequency, linked with a slightly higher turbulence length (  = = = = 6.3).
One possibility to push at higher frequency this annoying peak could be to use a shorter domain (e.g., 5 m); this has to be confirmed in a future work.

Correlations-Comparison with ESDU
Streamwise correlation length is evaluated thanks to Taylor's hypothesis: the temporal correlation (of fluctuating component u, v, and w) is calculated and compared with a negative exponential (Figure 10).The parameter of this negative exponential is given by Formulae (6a) (135 m at 10 m high, 105 m at 5 meter high), divided by the local mean wind speed (given on Figure 7).The agreement is very good for axial component , especially at the middle of the domain.For the other two components (v and w), the situation is not so satisfactory: for ESDU, correlation lengths are different for the three axes (Equation (6b) and (6c)), whereas in our LES calculation they appear to be the same.
Through the use of the horizontal array of probes, the lateral correlation can be directly measured.The central vertical array is used for vertical correlation.Results are plotted Figure 11.Lateral and vertical correlation lengths are much smaller than streamwise: we can estimate from Figure 11 around 5 m, which is a quarter of the width and height of the domain.These correlation lengths seem limited by the domain size much more than any other aspect; this assumption has to be confirmed in a future work.

Conclusions
A new solution for generating unsteady inlet conditions for Large Eddy Simulation is presented in this paper.This solution, which draws its inspiration from the wind tunnel testing community [14], is as simple as possible and straightforward to implement in OpenFoam.The flow field generated can be directly used as inlet conditions for Large Eddy Simulations.
In the central part of the flow (half centered span laterally, three fourth from ground vertically), the results are compared with international standards concerning wind engineering (Eurocode 1 and ESDU 75001).Densities are found symmetric, close to Gaussian, and the flow is homogeneous (no lateral fluctuations).Vertical profiles of mean wind speed and turbulence intensity fit well with the Eurocode model.The stream-wise correlation length is about 120 m which agrees with ESDU; lateral and vertical correlations are much smaller (about 4~5 m, limited by the computation domain size).The main discrepancies concern the power spectral density: the low frequency domain is a bit over energetic while there is a lack of high frequencies.This misbalance might be linked with the LES filtering of high frequencies.
Future work has to be conducted to verify some assumptions and check the reliability of the method with other characteristics.Authors do not plan to conduct a classic mesh sensitivity-study, because in LES, the physics change with the grid size (implicit filtering [15,16]).What seems important is the ratio between the grid size and the roughness cubes size.Thus, future work will address the problem of changing the roughness length and the whole domain size.

Figure 1 .
Figure 1.Geometry of the roughness-cubes lying on the floor.

Figure 4 .
Figure 4. Instant wind field slices in the computational domain.

Figure 7 .
Figure 7. Vertical profile of mean wind speed.Blue stared: LES results.Black line: Eurocode.

Figure 11 .
Figure 11.Lateral and vertical correlation for the three components.Blue: u-component, red: v-component; purple: w-component.