The Fifth International Symposium on Computational Wind Engineering (CWE2010) Chapel Hill, North Carolina, USA May 23-27, 2010

Development of a numerical code for 3D LES simulation of Thunderstorm Downburst

Kalyan Kumar Dasa, A.K.Ghosh b, K.P.Sinhamahapatra c

aLecturer(Selection Grade) Mechanical Engineering, Assam Engg. college, Guwahati, Assam, India,kkdas_aec@yahoo.co.in

bProfessor,Department of Aerospace Engg., IIT Kharagpur, Kharagpur, West Bengal, India, akghosh@aero.iitkgp.ernet.in

cAssociate Professor,Department of Aerospace Engineering,IIT Kharagpur,Kharagpur, West Bengal, India, Kalyanps@aero.iitkgp.ernet.in

ABSTRACT: Severe thunderstorms are important weather phenomena which impact on various facets of national activity like civil and defense operation, particularly aviation, space vehicle launching, agriculture in addition to its damage potential to life and properties. Experimental and numerical simulation studies on thunderstorm downburst have been reported by many researchers during the past two decades. Most of the numerical studies are based on commercially available codes with k-epsilon and 2D LES models for turbulence. In the present work a 3D numerical code has been developed using vorticity-vector potential formulation, with LES model for the turbulence. In addition a microburst simulator has been fabricated and used to generate experi-mental data for validation of predicted results from the code.

1 INTRODUCTION

The famous atmospheric scientist Fujita[1981] has observed and studied the flow due to downburst impacting on the ground and spreading outward in the different directions. He classi-fied downburst as either microburst or macroburst depending on their horizontal extent of dam-age. For the complexity of the full scale phenomenon, the physical simulation of the downburst is confined to the generic experiments of density currents impinging on a wall. Alahyari and Long-mire[1995], Lundgren et al.[1992], Cooper et al. [1993], Didden and Ho[1985], Knowles and Myszko[1998] have studied experimental simulation of the downburst. Letchford and Chay[2002], Chay and Letchford[2002] and Sengupta and Sarkar [2007] performed physical modelling to study the flow field characteristics and pressure distribution the stationary and translational downburst.

Numerical simulation of the downburst is performed by Proctor[1988],Craft et al. [1993] and Selvam and Homes[1992]. Kim and Hangan[2006] and Sengupta and Sarkar [2007] simulated the downburst flow field with different turbulence model using FLUENT software. The primary ob-jective of the present work is (a) to develop a 3D code based on the impinging jet model to study the wind field characteristics of the microburst flow field (b) to study the flow dynamics due to the variation of the H/Djet ratio and the Reynold’s number. (c) to validate the 3D CFD code with results from the physical simulation. (d) to study the radial velocity profile under partial slip con-dition imposed on the solid surfaces. 2.0 NUMERICAL MODELLING :The impinging jet model is considered for the numerical si-mulation of the microburst wind. The 3D N-S equations have been solved using vorticity-vector

The Fifth International Symposium on Computational Wind Engineering (CWE2010) Chapel Hill, North Carolina, USA May 23-27, 2010

potential formulation and the LES is used for the turbulence. The 3D incompressible Navier-Stokes equations in Vorticity-Vector potential formulation are as follows neglecting the second order terms[1987].

The three vortcity components are related to the velocity components by u . vector po-tential( ) is related to the velocity fields by u=curl .

………………………..(3)

………………………………….(4) ……………………(5)

The vector potential(ψ) is related to the vorticity(ξ) by equation 3. The vorticity transport equ-ation 1 is normalized with the jet exit parameters..The sub-grid scale model is given by equation 4 and . Vorticity transport equations (1) are solved by implicit ADI technique and the Poisson’s equations(3) for the vector potential are solved by SOR technique. Following Sakamota et al. [22], the Smagorinsky constant (Cs) is taken as 0.15. The density and molecular viscosity of air are taken as 1.225 kg/m3 and 1.7894×10-5 N-sec/m2 for the present numerical simulations. It is further assumed that the flow enters the computational domain where the fluid is stationary at t = 0 with the jet exit velocity. Square jet is considered for this work.

2.1 Computational domain and boundary conditions: Dimensions of the computational domain are 12Djet X12Djet X8Djet in the x(radial), y(spanwise) and z(axial) directions respectively as shown in figure 1. Uniform Cartesian cells are Considered for the computation. The computational grid consist of 96X96X64 nodes along the x,y and z directions. No-slip boundary condition is imposed for the solid surfaces (b) Gradient of all parameters are equated to zero at the outflow boundaries (c) For the jet walls free slip condi-tion is imposed (d) For the inlet boundary velocity is assumed to be axial and in the negative z di-rection(e) symmetric boundary condition is considered along the jet axis. In the numerical investi-gation three surfaces are considered having different slip conditions. Partial slip is the physical state between the no-slip and free slip conditions. Mathematically, slipping can be determined from the velocity of the surface layer fluid. In this numerical work three surfaces are considered where fluid velocity on the surface are 0% (no slip), 2.5% and 5% of the velocity of its neigh-bouring fluid layer. Henceforth, these surfaces will be represented by S1, S2 and S3 respectively. 3.0 PHYSICAL MODELLING A nozzle of diameter (D jet) 165 mm(6.5 in) is used to pro-duce the impinging jet. The size of the impinging platform is 1.5 m X 1.5 m. The distance of the

The Fifth International Symposium on Computational Wind Engineering (CWE2010) Chapel Hill, North Carolina, USA May 23-27, 2010

jet from the impinging platform (H) can be varied from 165 mm to 600 mm using an adjustable frame. The H/D jet ratios considered for this work are 1.0, 1.25 and 1.5. Henceforth, these ratios are termed as R1, R2 and R3 respectively. Three jet velocities(Vjet) 20 m/s, 25 m/s and 30 m/s are used in the experimentation. Henceforth, these velocities are termed as V1, V2 and V3 respec-tively. For the velocity measurements in the flow field hot wire anemometer system (DANTEC 56C17) is used. To determine the velocity profiles in the radial and axial directions hot wire ane-mometer probe is placed in the DANTEC traverser system. For the pressure measurement scani-valve PDCR23 pressure transducer system is used. Experimental setup is shown in figure 2. 4.0 RESULTS AND DISCUSSION

A grid dependency study is made initially to verify the algorithm and the code. The computed results using two separate grids for the flow with the jet velocity V1 at height/diameter ratio R1 are presented in figure 3. No significant change in the radial velocity profile u(y) is observed for the two grids. Similar comparisons are found in other cases as well. Consequently the grid with the resolution 96X96X64 is adopted for all subsequent computations.

Figure 4 shows radial velocity contour for the jet velocity V1 at R1. It show high radial veloci-ty closer to the impinging plateform. Similar results are also obtained for other jet velocities at dif-ferent height/diameter ratio. The maximum radial velocity is seen in the region of x = 1.1-1.5Djet. Axial velocity contour for V3 at R1, R2 and R3 are presented in figure 7, where a region of high axial velocity is seen closer to the axis of the downdraft and opposing axial velocities are seen in the region of primary vortex ring. It can be concluded from figure 7 that the primary vortex is closer to the impinging plate for lower value of the height/diameter ratio. Hence the region of high radial velocity will be closer to the impinging plate for lower value of height/diameter ratio. Figure 5 and figure 6 show the velocity vector plot showing the formation of the primary and the sec-ondary vortices due to initial Kelvin-Helmholtz instability. It is also seen from figure 6 that the primary and secondary vortices have opposite rotational velocity. Formation of the primary and secondary vortices can be seen in figure 13 also which shows the vorticity contour for jet velocity V3 at R1 and R2. In the free jet region the initial Kelvin-Helmholtz instability generates the vortex rings because of stronger shear generated at the start-up as surrounding fluid is at rest. In the im-pinging jet region the vortex rings touch down the surface and start moving radially triggering an unsteady separation of the newly formed boundary layer. The counter rotating primary and sec-ondary vortices moving radially outward produce maximum velocity closer to the surface. Figures 8 and figure 9 show the variation of radial velocity with jet velocity and height/diameter ratio at different locations from the point of impact. Clearly the maximum radial velocity decreases with an increase of distance from the point of impact. It can be seen from the these figures that the maximum radial velocity increases with increase in jet velocity and decrease in height/diameter ra-tio. Also the location of the maximum radial velocity moves closer to the impinging plateform for higher jet velocity and lower height/diameter ratio. Figure 11 shows the effect of variation of the radial velocity profile for partial slip condition in the range of 0-5%. No, significant variation in radial velocity profile is seen jet velocity of V1 at height/diameter ratio R1 and R2. Similar results are also obtained for different jet velocities in the same range of partial slip. Figure 10 show the results from the microburst simulator fabricated by the authors and computed radial velocity for various Reynolds numbers and height/diameter ratio. The figures show good agreement between CFD solutions and experimental results. Figure 12 shows the measured and computed radial ve-locity at a height of 0.1Djet from the plate for the configuration having = R1 at Reynolds

The Fifth International Symposium on Computational Wind Engineering (CWE2010) Chapel Hill, North Carolina, USA May 23-27, 2010

number of 85573 (Re1). The results are compared with the well-known correlation proposed by Holmes and Oliver [2000]. Holmes and Oliver [2000] reported that the radial velocity profile us-ing the empirical formula agrees well with the radar observation of the full scale phenomena by Hjelmfelt [1988] when radial length scale (R) is 50% of . Following the observation R is tak-en as 50% of for the present comparison. Figure 17 shows fairly good agreement between the current CFD and experimental results and the results from the empirical formula. However, the experimental and computational data indicate that the radial velocity in the stagnation region is not truly linear as given by potential flow solution and used in Holmes and Oliver [2000]. The subsequent decay follows the exponential law to a certain distance, nearly 1.5 for the cases studied here. However, beyond this distance the velocity decay rate reduces considerably due to the influence of the secondary vortex generated just outside and below the main vortex ring. Thus, the secondary vortex helps to increase the radial extent of the downburst wind. 5.0 CONCLUSION: An explicit second order accurate finite difference based LES method em-ploying 3D vorticity-vector potential approach (using FTCS scheme) has been developed. A mi-croburst simulator is also fabricated to validate the CFD code. Current CFD results from the axi-symmetric microburst model , 3D experimental results and results from the famous empirical relation of Holmes and Oliver [13] are compared and good agreement is observed between these results. The present study reveals the following facts regarding the microburst flow field,

(a) The maximum radial velocity of the microburst near the ground increases with the in-crease in Reynolds number and decrease in .

(b) The region of high radial velocity moves closer to the impinging plate for an increase in Reynolds number and decrease in . c) Microburst flow field characteristics depend on Reynolds number of the flow,

ratio. (d) The secondary vortex reduces the velocity decay rate beyond 1.5 and hence in-creases the radial extent of the downburst. (e) The maximum radial velocity and its location is insignificantly affected by slip condi-tion considered in the range of 0-5% for this work.

Figure captions: Figure 1: Computational domain Figure 2: Experimental set up for the physical simulation Figure 3: Grid independence test for the two grids Figure 4: High radial velocity near the impinging platform Figure 5 Velocity vector plot showing the formation of the vortex Figure 6: Velocity vector plot showing the formation of the primary and the secondary vortic es near the impinging plate at Vjet=V1 at R1

Figure 7 Showing the variation of the axial velocity(w) for different H/Djet values Figure 8: Radial velocity profiles at different locations for different H/Djet at fixed Vjet Figure 9: Radial velocity profile for different H/Djet values Figure 10: Comparison of experimental and computational velocity profiles Figure 11: Radial velocity profiles for different partial slip condition

The Fifth International Symposium on Computational Wind Engineering (CWE2010) Chapel Hill, North Carolina, USA May 23-27, 2010

Figure 12: Comparison of radial velocity at an altitude of 0.1 Djet Figure 13: Vorticity contour

Figure 1 Figure 2

Figure3 Figure 4

Figure 5 Figure 6

The Fifth International Symposium on Computational Wind Engineering (CWE2010) Chapel Hill, North Carolina, USA May 23-27, 2010

Figure 7

Figure 8

Figure 9

The Fifth International Symposium on Computational Wind Engineering (CWE2010) Chapel Hill, North Carolina, USA May 23-27, 2010

Figure 11

Figure 12 Figure 13

Figure 10

The Fifth International Symposium on Computational Wind Engineering (CWE2010) Chapel Hill, North Carolina, USA May 23-27, 2010

Nomenclature:

Re Reynolds number downburst diameter H height Si j strain tensor t time downdraft velocity ψx ψy ψz vector potential in x,y and z Δx, Δy and Δz mesh size in x, y and z direction

ξx , ξy and ξz vorticity in x,y and z directions ω relaxation factor

sub-grid scale eddy viscosity Cs Smagorinsky constant

References: 1. Alahyari, A., Longmire, E.K., 1995. Dynamics of experimentally simulated microbursts. AIAA J. 33

(11), 2128-2136. 2. Chay, M.T., Letchford, C.W., 2002. Pressure distribution on a cube in a simulated thunderstorm downburst—Part A: stationary downburst simulation. J. Wind Eng. Ind. Aerodyn. 90, 711-732.

3 .Cooper, D., Jackson, D.C., Launder, B.E., Liao, G.X., 1993. Impinging jet studies for turbulence model

Assessment-I. Flow-field experiments. Int. J. Heat Mass Transfer 36 (10), 2675–2684. 4. Craft, T.J., Graham, L.J.W., Launder, B.E., 1993. Impinging jet studies for turbulence model As-

sessment-II: an examination of the performance of four turbulence models. Int. J. Heat Mass Transfer 36 (10), 2685–2697. 5. Didden, N., Ho, C.M., 1985. Unsteady separation in a boundary layer produced by an impinging jet. J. Fluid

Mech. 160, 235–256. 6 .Fujita, T.T., 1981. Tornadoes and downbursts in the context of generalized planetary scales. J. At-

mos. Sci. 38, 1511–1534. 7. Holmes, J.D., Oliver, S.E., 2000. An empirical model of a downburst. Eng. Struct. 22, 1167–1172. 8 .Kim, J., Hangan, H., 2007. Numerical simulation of impinging jets with application to downbursts. J.

Wind Eng. Ind. Aerodyn. 95, 279–298. 9. Knowles, K., Myszko, M., 1998. Turbulence measurement in radial wall-jets. Exp. Thermal Fluid

Sci. 17, 71–78. 10.Letchford, C.W., Chay, M.T., 2002. Pressure distributions on a cube in a simulated thunderstorm

downburst, Part B: moving downburst observations. J. Wind Eng. Ind. Aerodyn. 90, 733–753. 11.Lundgren, T.S., Yao, J., Mansour, N.N., 1992. Microburst modeling and scaling. J. Fluid Mech. 239, 461–488.

.12.Proctor, F.H., 1988. Numerical simulations of an isolated microburst. Part I: Dynamics and Struc-ture. J. Atmos. Sci. 45, 3137–3160 13.Sakamota, S, Murakami S., Mochida A., 1993. Numerical study on flow past 2D square cylinder by Large Eddy Simulation Comparison between 2D and 3D computations, J. Wind Eng. Ind. Aerodynam-ics 50 (1993) 61-68. 14.Selvam, R.P., Holmes, J.D., 1992. Numerical simulation of thunderstorm downdrafts. J. Wind Eng. Ind. Aerodyn. 41–44, 2817–2825

15.Sengupta, A., Sarskar, P. P., 2007. Experimental measurement and numerical simulation of an impinging jet with application to thunderstorm microburst winds. J. Wind Eng. Ind. Aerodyn. (2007)

16. Hjelmfelt, M.R., 1988. Structure and life cycle of micoburst outflows observed in Colorado. J. Appl. Met. 27, 1988, 900-927