Numerical Study of the Flow Fields in Downburst with Consideration of the Rough Condition on the Ground

The downburst is an extreme weather condition whose resulting load will affect the stability of the building structures. The characteristic of downbursts is required during the design of buildings. In order to achieve the characteristics of the downbursts, large eddy simulation (LES) is adopted. The method impinging jet is used to reproduce the downbursts, meanwhile smooth and rough ground conditions are examined. The setting of the rough layer of the ground is done by using the immersion boundary method (IBM). After the simulation, the wind field is decomposed into the mean component and the turbulence component. In this paper, the radial wind speed profile and the wind speed time diagram of the simulation experiment are compared with the previous measured data and the simulation results. This paper reveals that the radial wind speed is the key part of the downburst wind speed, and it rapidly increases with height. It is also found that the rough ground will cause the peak radial velocity to move up, which is consistent with the change of the main annular vorticity in vortex cloud image analysis. Finally, the turbulence intensity is found to be relatively small at the position where the radial wind speed is relatively large.


Introduction
The definition of downburst was first proposed by Fujita [1] that in a thunderstorm, a very destructive strong wind will be formed by the divergence of the air from the impact center to the circumferential direction along the surface after the rapid descent of the airflow and impact on the ground.The formation of downburst is closely related to thunderstorm weather.Its source is also one of the manifestations of the sudden release of unstable energy in the atmosphere, but it is not the same as the common wind phenomenon.The downburst is a great threat to the takeoff and landing of all kinds of aircrafts in the field of aerospace.And it is also one of the important factors of thunderstorm wind load in the field of designing building structures.Its impact on human production and life in many areas has gradually become a hot issue of research.
Tornadoes and downburst are extreme weather events in thunderstorms.The former is well-known because of the horrible phenomena and great destructiveness.However, the power of the latter one also cannot be neglected.Many disasters caused by it are not uncommon at home and abroad.Hawers [2] et al. selected and analyzed 94 transmission towers damage cases in Australia.It was found that 90 percent of the damages were caused by downbursts and tornadoes and other extreme weather phenomena.In 1996, a downburst occurred during a thunderstorm in Ontario, Canada.This led to the collapse of 19 transmission towers and the economic losses amounted to more than 1000 million dollars.Due to geographical environment and other factors, China is also a country where thunderstorms occur frequently and disasters caused by downburst also occur frequently.In June 22, 2000, a passenger plane crash occurred when it was passing Wuhan, Hubei and there are 42 victims in the crash.The meteorological bureau investigated the local weather conditions, combined with the accident scene investigation, and finally concluded that the downburst was the main culprit of the disaster.June 14, 2005, due to the occurrence of downburst weather, near Siyang, Jiangsu, the transmission towers were destroyed and the power consumption was affected in a large area.June 2009, in Zhenjiang, Jiangsu, five base transmission tower collapsed.And the research results showed that it was related to downbursts.It is obvious that the extreme weather such as downbursts may cause great damage to the infrastructure, and it also has a huge impact on transportation, threating people's life and economic interests and other aspects.
At present, scholars mainly use three research methods to study it.It is the wind field measurement research based on historical records of meteorological data, analysis of meteorological theory based on the generation and development principle of downburst and the laboratory simulation based on the special computer.These three methods have distinctive features, but they complement each other, which play vital roles for researchers to understand this kind of special weather phenomenon.
The direct simulation of thunderstorm weather is a relatively close approximation to the actual situation, but its calculation is huge and it is difficult to implement.From the past decade, scientists mainly use IJ (Impinging Jet) and CS (Cooling Source) methods to do the simulation.The advantage is that it can reduce the huge amount of computation.In addition, it has good consistency with the real cases.The IJ method referring to the numerical simulation of impinging jet is in the paper of Fujita [3].CS method is a physical simulation method of using cooling source proposed by Anderson [4].Kim and Hangan [5] et al. used IJ method to establish axisymmetric two-dimensional computational domain to simulate the downburst, and the average wind speed was compared with the real case.The results were good.Sengupta and Sarkar [6] used k-epsilon, k-omega turbulence model, Shear stress transfer (SST) model and LES turbulence model to simulate the downburst.The simulation results were compared with the actual measurement results, which confirmed that the simulation method is consistent with the actual downbursts.Mason [7][8][9] used Shear adaptive simulation (SAS) to do the simulation of downbursts.However, it was concluded that the turbulent viscosity of the jet is too high when dealing with the simulation results.Then they also studied the effect of impinging jet inclination and ground roughness on the wind field simulation.At the same time, the impinging jet device was set up in the laboratory, and the data were recorded as a comparison.The reliability of simulation is proved.Vermeire [10] et al. used CS method to simulate the influence of different terrain conditions on the characteristics of the downbursts wind field.Haitham [11] et al. simulated the downbursts in actual space size.At the same time, the roughness of the ground corresponding to four different realistic ground structures was added.Honghai Li [12] used the method of wall jet, choosing k-Epsilon turbulent model to simulate the downbursts.Meanwhile, the scaled building model was added into the calculation domain.The wind pressure distribution on the surface of the building is obtained, and the wind load value of the building is analyzed.Wenfu Zhang [13] et al. combined the DESH model, Wood vertical wind profile model and Holmes horizontal radial wind profile model to simulate the downbursts.The results turned out to be good, but the time variation was not considered in the calculation of correlation function.Weilian Qu [14] et al. used the CFD to simulate the evolution of downbursts and divided it into three stages: formation, development and steady states.The radial wind velocity profiles at different positions were analyzed, and the effects of different initial parameters were taken into account.Zhiwei Peng [15] also used CFD method to simulate the downburst characteristic, and the influence of motion of downbursts was considered.
Looking for ways to prevent and cope with downburst from structures design angle, and to reduce or avoid its harm and loss, it is necessary to take the downburst directly as an entry point to understand its characteristics.Therefore, this paper will reproduce the process of evolution and development in a certain space by means of the simulation methods which is mentioned above in the three research ways.
Different researchers have different initial parameters such as the size of the calculation domain and the velocity of the impinging jet.In this paper, the modeling method of Haitham [11] is taken as reference.In the research of Weilian Qu [14], the actual time is used to select the radial velocity profile, but this paper will be based on dimensionless time for analysis, being more universal.At the same time, many scholars did not consider the influence of ground roughness.The difference between previous research and study of this paper is shown in Table 1.Based on the immersed boundary method (IBM) adopted by Zhenqing Liu [16] et al., this paper investigates the influence of two kinds of ground conditions on the wind field, which is smooth and rough respectively, and further investigates the characteristics of the wind field.

Wind Field Model of Downburst
Only by choosing the appropriate model of the downburst wind field, can the simulation result be similar to the real case, which will be more convincing in comparison with the wind field measurement results.At the same time, different wind models need to select different wind speed profiles to analyze the wind field characteristics, so as to make the results clearer.The wind field model of the normal wind in the atmospheric boundary layer is relatively simple.It can be simplified into uniform unidirectional horizontal flow.The wind speed increases with the increase of height from the ground.Besides, the wind profile is easy to select and the characteristic is more obvious.Due to the development process of downburst is special, the wind field become more complex.In the process of simplifying the wind field model, because it is non-uniform flow, the magnitude and direction of wind speed vary greatly in each point of wind field.It needs to consider the horizontal speed and vertical speed characteristics and do radial and vertical wind profile analysis.
Recently, there are two types of numerical simulation of downburst.One is annular eddy model as shown in Figure 1(a) [14].Another is impinging jet model as shown in Figure 1(b) [14].Impinging jet model is to jet the flow vertically to the ground surface.The annular vortex is also produced in the process of flow sinking.The results show that the synergy effect is good for the simulation of the wind field characteristics of downburst [18,23], making up for the unsatisfactory numerical effects in the core and peripheral areas of wind field simulated by annular eddy model.Therefore, the impinging jet model is chosen in this paper.

Governing Equation
This paper uses the large eddy simulation to model the wind field.As discussed above, in the research, only the large scale eddy will be directly simulated.The small scale vortices which below the grid resolution will be considered by subgrid model to simulate their effects on large eddy.Boussinesq hypothesis and the standard Smagorinsky-Lilly model are used to calculate Subgrid stress (SGS).
The time after which the subgrid vortex are filtered depends on Navier-Stokes Equation (Cartesian coordinate system), as follow: Where u ̃ is velocity of fluid that filtered subgrid vortex, p ̃ is the pressure and μ is air viscosity coefficient.The standard Smagorinsky-Lilly model is used to calculate the subgrid stress as follows: Where Ls is the mixing length of subgrid, κ is karman constant, set as 0.42, d is the distance from wall grid to wall surface, V is the volume of control body, Cs in this model is taken as 0.1, which is commonly used in the atmospheric boundary layer simulation.
In order to compare the wind field with the smooth ground condition, this paper chose the immersed boundary method (IBM) to realize the simulation of the ground area covered by buildings.In this method, the resistance term caused by ground roughness is directly added in the Navier-Stokes Equation.So steps of establishing geometric shape of rough element is effectively avoided.Formula is as follow: , () The calculation formula of resistance term is from the studies of Enoki and Ishihara, as follow: , 1 2 Where Cd is the resistance coefficient Cd=Cfiγ0/l0, γ0 is the rough element density, l0 is the roughness length, u ̃mag represents the incoming wind speed.In this paper, Cd is taken as 0.4 according to the previous research by Kaimal [24].
The finite volume method is used for spatial discretization.Two order central difference scheme is applied to convection term and viscous term.The two order implicit scheme is used for the unsteady term time advancing.SIMPLE algorithm is used for pressure and velocity decoupling, see Ferziger [25].Fluent 6.3.26[26] is used as the solver.

Geometry Model
A complete 3D cylinder is chosen as the computational domain, as shown in Figure 2.This makes the evolution of flow field not to be limited, and the simulation results more accurate.The geometric scale is selected as 1:1000.The diameter of impinging jet Dj is set as 1m.The size of the selected cylinder domain is shown in Figure 2(a).The computational domain diameter is 8Dj.Vertical height is 4Dj.The impinging jet inlet is located at the vertical height of 2Dj.The impinging jet channel is also a cylinder, extending to the upper circular surface of the cylinder.According to the wind speed characteristics of downburst and the selection of computational domain size, the speed scale 1:3 recommended by Mason [8] is used in this paper.The velocity of impinging jet is taken as 20m/s.For rough surface conditions, the resistance term in the momentum balance equation is added to the range within 10cm to the ground surface.

Mesh Distribution and Boundary Conditions
The 3D structural grids are meshed by GAMBIT.The total amount of meshes is 594000.In the process of meshing, to ensure the good quality of grid structure, firstly, the grid is divided into squares at the center of the horizontal section; then, the grid is divided along the circumferential direction; finally, the volume meshing is done.The height of the grids which are closest to the ground is 0.1 mm.The maximum mesh height is 10 cm in the whole calculation domain.The minimum horizontal grid size is 1 mm while the maximum horizontal grid size is 10 cm.The growth rates of vertical and horizontal grids are less than 1.05.The cylindrical surface intersected by the impinging jet channel with computational domain is the sliding wall without shear stress to avoid the impinging jet's influence before entering the computational domain and forming downburst.In the cylinder computation domain, the bottom boundary is wall condition, and the cylindrical surface as well as the upper surface are pressure outlets.The calculation time step is 0.0001s, the total computation time is 50s.

Data Comparison
Normalization processing and wind speed profile are made to compare the present numerical data with previous wind field measurement as well as numerical results, shown in Figure 3.It can be clearly seen that the numerical simulation is in good agreement with the previous research results [5,10,11,27,28].Thus, the reliability of the turbulence model, the mesh distribution, the boundary condition and the time step are verified.

The Evolution of Wind Field with Time
Because of the symmetry of the wind field, a vertical section is selected in this paper.Figure 4 is a vorticity image of a downburst wind field at selected characteristic moments.It can clearly show the evolution process of the flow field.In Figure 4(a), the main annular vortex produced by the impinging jet begins to impact the ground, and secondary vortices have already been generated at the ground boundary.In Figure 4(b), the annular eddy is impacting the ground, and the secondary vortex is generated at the ground boundary.Figure 4(c), the main annular vortex sheds, and the displacement is in the radial direction.In Figure 4(d), it is the main annular vortex.In Figure 4(e), the main annular vortex has larger radial movement and moves vertically.

Wind Field Analysis
In order to reduce the influence of coordinate system on data, the velocity is divided into two parts: axial velocity and radial velocity.And this paper only focus on the analysis of complex radial velocity characteristics as well as the wind speed profile characteristics on four particular positions.The wind speed profile was drawn on four particular positions in the radius from 1Dj to 2.5Dj with the most complex radial wind speed variation, as shown in Figure 5.The radius of R were 1Dj, 1.5Dj, 2Dj and 2.5Dj, respectively.The time point for drawing the wind profile is that the maximum radial velocity at the four positions are reached.Tn is 7.8, 8.8, 11 and 12.8, respectively, corresponding to the vorticity maps at the same time points mentioned above.It can be clearly seen that when the R=1Dj, the maximum radial wind speed is reached, and the main annular vortex center has passed this position, but the vortex does not completely leave this position.At this time, the upper part of the profile is close to the vortex center, so the wind speed is lower.The lower part of the profile is influenced by the following vortex of the main annular eddy impacting the ground, so the wind velocity is higher.When R=1.5Dj, the maximum radial velocity is reached and the profile just passes through the vortex center, which is the reason why the wind velocity is negative at the upper end of the profile.Secondary vortex exists in the lower end of the profile, and the nearby wind speed increases to 1.6Vj.After the downburst hitting the ground, the flow moving along the ground boundary is generated.This flow is influenced by the negative pressure of the annular vortex.After violent change, the secondary vortex is formed which is opposite to the direction of main vortex.The wind speed near the ground is greatly enhanced under the synergistic action of the main annular vortex and the secondary vortex.When R=2Dj, the maximum radial velocity is reached.There is a marked upward movement of annular vortex.At the same time, the intensity of secondary vortex decreases gradually.Thus, the maximum velocity is lower than that in R=1.5Dj.Meanwhile, the upward movement of the annular vortex also leads to the wind velocity of some place, which is not enhanced by secondary eddy, become lower.It is shown in the figure as the wind speed curve concave part.However, at the R=1Dj, the increase of wind speed is the result of that the secondary annular vortex impacts the ground.Because R=2Dj is far from wind field center, the energy dissipation and the maximum wind speed are also lower.The secondary vortex is obviously dissipated, so the near surface wind speed is not obviously enhanced.

Influence of Surface Roughness
Figure 6 is the maximum wind speed envelope graph which contains all of the time points from the height 0 to 0.3Dj.It is clear to see the effect of ground roughness on maximum radial wind speed.When the roughness term is not added, the maximum radial wind speed reaches 1.8Vj and at the height of 0.025Dj.However, when the roughness term is added, the near surface wind speed deceases to 0.2Vj, and steep slope appears in the near ground range.Then the wind speed increases rapidly with the height to 1.4Vj, distributing uniformly from 0.1Dj to 0.18Dj.Combined with the established part of the model, it can be analyzed that the speed steep slope is caused by adding the surface roughness layer.Its height is consistent with that of the rough layer.After the maximum velocity was observed, from the slope of velocity decreasing with the increase of height, it can be found that the slopes of the two cases are similar, indicating that the general characteristics of the attenuation of the downburst with altitude.In summary, after increasing the roughness of the ground, the maximum radial wind velocity decreases and the location of occurrence moves upward.This is consistent with the experimental results of CS physical simulation by Mason [27] et al.
Figure 7 is the contour map drawn after the statistics of the maximum radial velocity of all points from vertical range in 0≤Z≤0.2Djand radial range in 0.5Dj≤R≤3Dj.There are two peaks in Figure 7(a), which are located radially within the 0.9Dj to 1.0Dj and 1.2Dj to 1.6Dj, respectively.And the distribution range of the former peak is less than the latter.Combined with the vorticity map above, it can be explained that the place between the two peaks is the center of annular vortex impacting the ground.Therefore, larger radial wind speed appears at both sides of the location.The distribution of the latter one is bigger due to the enhancement of the secondary vortex and the vertical movement of the vortex.Under the rough condition, the maximum radial wind speed is lower than the smooth situation in general.This is consistent with the previous conclusions, and peaks are located within 1.1Dj to 1.4Dj and 0.07Dj to 0.18Dj in a long narrow shape.The different distribution results are caused by the different moving paths of the main annular vortices in two kinds of ground conditions.In order to study the turbulent, the turbulence intensity is chosen as the basis.The calculation method of turbulence intensity is the same as that of Holmes [22] et al.The formula is as Equation 9.This method is equivalent to regard the turbulent part of the downburst as a piecewise stationary process in a short period of time before and after the maximum speed is reached.
in which Ur maxis the maximum radial velocity at the position, σur max is the mean square deviation of the wind speed fluctuation in time interval from tmax-(1/2fcut) to tmax+(1/2fcut), tmax is the time that the radial wind speed is reached in the location, in which fcut is the cut-off frequency.The calculation formula is fcut=2×0.3×Vj/Dj=0.6×Vj/Dj.After the calculation of turbulence intensity is done, the contour distribution map is shown in Figure 8.It is known that the turbulence intensity near the ground is large while the turbulence intensity is low when it is far away from the ground.For the two different ground conditions, smooth and rough, the turbulence intensity at the place near the ground where maximum radial mean velocity is reached, is roughly between 0.16 and 0.20, and between 0.16 and 0.24, which is slightly higher than the results of Aboshosha [13].The turbulence intensity peak is between radial 1.6Dj to 1.7Dj.It can be contrasted in Figure 7 and 8 that in the higher radial wind speed position, the turbulence intensity is smaller, indicating that the maximum radial wind velocity in the wind field is mainly related to the laminar flow.

Conclusions
In this paper, the CFD numerical simulation of the 1:1000 scaled downburst wind field is completed by the LES method.The evolution process and the wind field characteristics of the downburst are studied, and the following conclusions are obtained:  In the smooth ground condition, the peak, maximum and instantaneous radial velocity wind profiles in the wind field are consistent with the previous research by wind field measurement, laboratory wind tunnel experiment and numerical simulations.Thus, the reliability of the simulation results is proved.
 The main reason for the complexity of downburst wind field near the ground surface is because of the dramatic changes in the flow field after the annular vortex impacting the ground.The roughness of the ground will aggravate the evolution intensity and speed of the wind field.
 The axial velocity is mainly influenced by flow path of impinging jet.The axial velocity is larger from the outflow region of impinging jet to the ground.Its velocity direction is the same as the impinging jet, in the negative direction, and the wind speed is also similar.The positive axial velocity mainly appears in the outer part of the vortex ring axis, and is relatively small.
 The radial velocity changes very sharply.It is greatly influenced by annular vortex and secondary vortex.The main feature is the rapid increase near the ground.The maximum radial wind speed can reach about 1.6V j .The locations where the maximum wind speed appear move up slightly with the distance away from center point increasing.The increase of the roughness of the ground can also because the locations where maximum wind speed appear have a larger upward motion.

Figure 2 .
Figure 2. LES computational domain (a) and Mesh distribution (b)

Figure 3 .
Figure 3. Contrast chart of wind speed

Figure 4 (
f) is the final stage of the wind field.The vortex gradually dissipates.