A Mathematical Model Governing Tornado Dynamics: An Exact Solution of a Generalized Model

1 Introduction

A tornado is a fast whirling columnar vortex wind system hanging as a pendent from a cumuliform cloud in contact with the surface of the earth. It is witnessed as a funnel merging into clouds with circulating dust at the foot. This low-pressure visual core vortex rotating with terrific energy also travels along the ground surface. It is the most violent and destructive atmospheric vortex on the surface of the Earth. It is observed all over the world such as in Japan, Bangladesh, Britain, and Australia, but in large number in Tornado Alley in the USA. It is generally believed that the circulation of a tornado vortex is maintained by the rotating mother cloud and the sense of the vortex rotation coincides with that of the associated mother cloud [1].

Numerous theoretical and empirical models are available in the literature. The tangential velocity in tornadoes are generally approximated by continuous functions that are zero at the center of the tornado, increase to a maximum at some radial distance, and then decrease asymptotically to zero at points infinitely distant from the center. Different forms of the tangential component of velocity have been proposed using the idealized Rankine combined vortex model [2] as a first approximation. However, because of the absence of the radial and vertical components of velocity, the Rankine combined vortex model is doubtfully a good approximation.

In the literature, a number of analytical, empirical, and numerical models are available, discussing the flow field of mature tornado-like vortices for single cells and two cells above the tornado boundary layer. Burgers [3] and Rott [4] both individually presented the exact solution of full Navier-Stokes equations for viscous incompressible steady flow of a vortex with the radial velocity as u = −ar (a being the constant of proportionality), the azimuthal velocity v=Γ(1−exp(−ar2/2ν))/2πr (Γ being the circulation and ν the kinematic viscosity), and the axial velocity w = 2az. While the azimuthal velocity has been found to fit well to some observed and experimental data by several researchers (Wood and Brown [5], Kim and Matsui [6], Gillmeier et al. [7], etc.), the model may not be useful because other components of velocity are unbounded. Other models available in the literature are either empirical or are unable to model a real tornado.

Sullivan [8] gave an exact solution with some similarity to the Burgers-Rott vortex model. He discussed both one-celled and two-celled vortices. The two-celled vortex possesses an inner cell in which wind descends from above and flows outward to meet a separate wind that converges radially. Both the winds rise at the meeting point. The Sullivan vortex is a very simple vortex but can describe the flow in an intense tornado having a central downdraft, and also its updraft is localized to a place meant for the thunderstorm. The axial pressure gradient however increases vertically without bound. The two-cell analytical Sullivan [8] vortex model for steady incompressible viscous flow is quite complex.

Kuo [9] analytically modeled the three-dimensional flow in the boundary layer of a tornado-like vortex and alternatively solved the two nonlinear boundary-layer equations for the radial and vertical velocities. The Bloor and Ingham vortex model [10] and the Vyas-Majdalani vortex model [11] are exact solutions for inviscid flows using the Euler’s equations, respectively, in a conical and a cylindrical domain. Xu and Hangan [12] used a free narrow jet solution combined with a modified Rankine vortex to analytically model an inviscid tornado-like vortex. However, the combined model is not an exact solution to the Navier-Stokes equations. Wood and White [13] reported a new parametric model of vortex tangential-wind profiles, which is based on the Vatistas et al. [14] model and is mainly designed to represent realistic-looking tangential wind profiles observed in atmospheric vortices.

Tornado-like flow field has been studied experimentally and/or numerically in numerous reports. A few of them are as follows: Ward [15] simulated in a laboratory system the three features of tornadoes, viz., characteristic surface pressure profile, bulging deformation on the vortex core, and multiple vortices in a single convergence system; Church et al. [16] discussed the dynamics of natural tornadoes based on laboratory simulations; Lewellen et al. [17] simulated tornado’s interaction with the surface; Natarajan [18] discussed large eddy simulations of translation and surface roughness effects on tornado-like vortices; Sabareesh et al. [19] also discussed surface pressure and surface roughness, whereas Liu and Ishihara [20] studied translation and roughness on tornado-like vortices; Haan et al. [21] designed a large tornado simulator for wind engineering applications; Mishra et al. [22] physically simulated the single-cell tornado-like vortex; Refan et al. [23] tried to reproduce tornadoes in a laboratory using proper scaling; Gillmeier et al. [24] analyzed the influence of tornado generator’s geometry on the flow field; Nolan et al. [25] worked on tornado vortex structure, intensity, and surface wind gusts in large eddy simulations; and Tang et al. [26] worked on the characteristics of tornado-like vortex simulation.

Vatistas [27] experimentally observed that in the concentrated vortex, the azimuthal velocity component does not vary strongly in the axial direction. Therefore, with those assumptions, the radial velocity component can be obtained from the θ-momentum equation. This approaches Rankine model as a parameter denoting the sharpness of the velocity profile near the radius of the maximum wind increases infinitely (the details are given the main text). The Vatistas et al. [28] model is a generalization of a few well-known vortex tangential-velocity models. The Vatistas et al. [14] proposed the tangential velocity profiles for vortices with continuous distributions of flow quantities. Recently, Gillmeier et al. [7] have reviewed some classical analytical tornado-like vortex flow field models.

The full-scale structure of tornado is highly complex, and therefore, several issues such as instabilities, singularities, and nonlinearity pop up to be addressed (Lewellen [29], Alexander and Wurman [30], and Karstens et al. [31]). Hence, to understand the complete physical processes adhering to the tornado flow field, simplified mathematical models are required, which minimize the error existing in the full-scale data observations and allow to explain original velocity and pressure fields. Collecting real wind field data from tornadoes in nature has been a difficult task for observers because of their destructive nature. Nowadays, researchers use Doppler radars for enabling full-scale tornado data from a safe distance. However, the data collected from such measurements are largely limited to the genesis of tornadoes (Bluestein et al. [32], [33], [34], Lund and Snow [35], Wurman et al. [36], and Jones et al. [37]).

Yin and Chang [1] opine, “Tornado is a huge vortex column with a low pressure visual core,” and with this consideration, Pandey and Maurya [38] floated a mathematical model for atmospheric vortices by assuming an annular vortex. However, despite the fact that the characteristics discovered hold for all whirlwinds, it was discussed in detail with regard to dust devils.

Baker and Sterling [39] have recently published a paper in which for inviscid vortex flows they assumed the dimensionless radial velocity as u¯=−4δr¯z¯(1+r¯2)(1+z¯2), where δ is the ratio between the vertical and horizontal length scales and the rest bear the usual meaning. The other two components of velocity have been deduced as v¯=Kr¯ln(1+z¯2)(1+r¯2) (K being a constant) and w¯=4δln(1+z¯2)(1+r¯2)2 from the Euler equation. The viscous effects remained untouched.

With the objective of getting new exact solutions to the equations governing the dynamics of tornadoes duly considering viscous effects, we target to model single-cell tornadoes by considering radial velocity whose variation in vertical direction is consistent with the flow field in the laboratory vortex simulator (Ward [15]). We aim to deduce vertical and azimuthal velocities and pressure as well with due consideration to the inferences made by Makarieva et al. [40], who concluded, “The decrease of pressure along the vertical axis sustains the ascending air motion with vertical velocity w and induces a compensating horizontal air inflow with radial velocity u. The converging radial flow has maximal velocity at the surface, where the magnitude of the condensation-induced pressure drop is the largest. Radial velocity approaches zero at a certain height z = h, which approximately coincides with the cloud height.”

2 Physical Model and Mathematical Formulation

Equations governing the motion of a steady incompressible Newtonian viscous fluid with axial symmetry are given by

where u, v, and w are, respectively, the radial, azimuthal, and vertical components of fluid velocity and r and z are the radial and axial coordinates; p stands for the pressure, ν for the kinematic viscosity, and ρ for the density. The body forces (buoyancy) in the vertical direction are denoted by F_z.

Variables in (1)–(4) are made dimensionless with primed notations as follows:

where r_m is the core radius and v_m is the maximum azimuthal velocity.

In view of (5), (1)–(4), on dropping the primes, are transformed to

where Re = r_mv_m/ν denotes the Reynolds number.

3 Radial Velocity

The winds blowing toward the center of the tornado in the radial direction below few meters (or kilometers) of height play a crucial role in the formation of tornadoes. We consider the radial velocity of the form u(r, z) = −U(r) F(z), where U(r) is a function of only r and F(z) is a function of only z and the negative sign indicates that the radial wind blows toward the center.

It is mostly considered that the strong inflow weakens along the radial direction and reduces to zero at the center of the tornado. When the wind starts to rotate about the axis of rotation, it is physically accepted that the radial velocity diminishes with height and approaches zero at a certain height z = h, which approximately coincides with the cloud height (Makarieva et al. [40]).

We therefore assume that the radial velocity decreases linearly with height and vanishes at height z = h, i.e. F(z) = 1 − az/h, where the parameter a controls the shape of F(z) and a must be 1 when z tends to h so that F(z) vanishes. This is consistent with the flow field in the laboratory vortex simulator [15], [16]. Hence, the radial velocity may be taken as

For one-cell tornado model, Burger [3]-Rott [4] considered the radial velocity proportional to radial distance, which is physically unacceptable because it has no upper bound in the radial direction.

In the opinion of Vatistas [28], concentrated vortices produced in air and water follow a certain relationship. That relationship inspired Vatistas et al. [14] to assume the tangential velocity (i.e. azimuthal velocity) of the form v = r/(1 + r^2β)^1/β, where β governs the sharpness of the velocity profile near the radius of the maximum wind. β→∞ leads to Rankine’s model. Smaller integral values of β give smoother turns.

Then with the assumption that the radial velocity does not depend strongly on the axial coordinates, Vatistas et al. [14] derived from momentum conservation equation that the dimensionless radial velocity is U(r) = −2(β + 1)r^2β−1/(1 + r^2β). In that model, U is well behaved except for the numerical value of β < 1, where U has a singularity near the tornado center. For simple one-cell vortex, we accept it for β = 1, and so the dimensionless form of U(r) is U(r) = 4r/(1 + r²). Thus, the final form of the radial velocity may be considered as

where the second factor is based on the conclusion made by Makarieva et al. [40]; i.e. the pressure drop along the vertical axis sustains the ascending air motion with vertical velocity and induces a compensating horizontal air inflow with radial velocity. The converging radial flow has maximal velocity at the surface, where the magnitude of the condensation-induced pressure drop is the largest. u→0, at a certain height z = h, which is approximately the cloud height. Further, the coefficient 4 has been dropped as it does not modify the model qualitatively, and further, we shall be deriving other components from mass conservation and momentum conservation equations.

4 Vertical Velocity

Substituting (11) into the continuity equation (4) with the consideration that w(r, z) = w₁(r) × w₂(z), we obtain the vertical velocity as

where K is an integrating constant and is determined by using the boundary condition that vertical velocity w₂(z) = w₀ at z = 0. This yields

5 Azimuthal Velocity

6 Pressure

6.3 Determination of Pressure for Viscous Flow

Substituting u from (11), v from (28), and w from (13) into (6), the radial pressure gradient may be given by

(35)−∂p∂r=r(1+r2)3[(1−r2)(1−azh)2+2ah(z−az22h)−8Re(1−azh)]−v2r,

which, on performing integration from 0 to r, gives

(36)p(r,z)−p(0,z)=−r22(1+r2)2(1−azh)2+12{1−1(1+r2)2}{ah(z−az22h)−4Re(1−azh)}+∫0rv2(s,z)sds.

Similarly, from (11), (13), and (8), the axial pressure gradient may be obtained as

−∂p∂z+Fz=(z−az22h)4(1+r2)4{(1−azh)(1+2r2)+4(1−2r2)Re}+2aReh(1+r2)2,

integration of which from 0 to z yields

(37)p(r,z)−p(r,0)=−4(1+r2)4{(1+2r2)2(z−az22h)2+2(1−2r2)Re(z2−az33h)}−2azReh(1+r2)2+∫0zFsds.

From (37) at r=0, we get

(38)p(0,z)=p(0,0)−2{(z−az22h)2+4Re(z2−az33h)}−2azReh+∫0zFs(0,s)ds.

Figure 1:

The radial profile of the radial component u(r) of velocity based on (11). The parameters assumed here for the diagram are h = 50 and a = 1.

Figure 2:

(a) The radial profile of the vertical velocity w and (b) the dependence of w on z. The plots are based on (13), and the parameters assumed here for the diagram are h = 50 and a = 1.

Substituting (38) into (36), we get the pressure given by

(39)p(r,z)−p(0,0)=−r22(1+r2)2(1−azh)2+12{1−1(1+r2)2}{ah(z−az22h)−4Re(1−azh)}+∫0rv2(s,z)sds−{2(z−az22h)2+8Re(z2−az33h)}−2azReh+∫0zFs(0,s)ds.

7 Results and Discussion

The present model is more general in the sense that all the three components are also dependent on the axial coordinate. The derived velocity components are plotted to study the radial and axial profiles and the impact of viscosity thereon.

7.2 Vertical Component of Velocity

The vertical component of velocity has been derived from the continuity equation with substitution of the radial component designed in Section 2. The plots are given in Figure 2.

Figure 3:

Plots for (a) the radial profile of the azimuthal velocity v and (b) the dependence of v on z. The plots are based on (28), and the parameters assumed here for the diagram are S = 0.50, h = 50, and a = 1.

Figure 4:

Plots for vertical profiles of the azimuthal velocity v based on (28). The parameters assumed here for the diagram are S = 0.98, h = 50, and a = 1.

The radial profile of the axial component of velocity reveals that it is maximum at the centerline, weakens as we move radially outward, reduces at large distances from the axis, and then vanishes gradually (Fig. 2a). Baker and Sterling [39] showed very similar results with their derived dimensionless vertical velocity.

Vertical velocity increases along the axis (Fig. 2b). However, it diminishes as we move outward from the center. Similar patterns with similar radial profiles were observed by Liu and Ishihara [24] in their numerical study for tornadoes with swirl ratio of 0.02 for weak vortices. As swirl ratio is not a parameter in our formula, it is not possible to compare with other observations or results given in that paper. Therefore, a further improvement is required considering some other aspects of the vortex motion.

7.3 Azimuthal Component of Velocity

The azimuthal velocity is the most significant component in a whirling motion. In this investigation, it has been derived based on the assumption made for the radial component and also the axial component that was derived from the continuity equation by the substitution of the radial component of velocity. Unlike this, most of the previous theoretical investigations showed its dependence on the radial coordinate only ignoring its obvious variation in the axial direction. However, during numerical simulation, this fact could not be ruled out. Duly considering this aspect, Tang et al. [41] simulated azimuthal velocity for various fixed values of axial coordinates and found it to match relatively less general models of Rankine [2] and Burgers [3]-Rott [4]. As an exception to most of the previous results, Baker and Sterling [39] formulated the radial as well as axial dependence of the azimuthal velocity but for inviscid flows (equation (9) and Figure 1b in Baker and Sterling [39]). The results they obtained are similar to those of our model for the inviscid part. However, the impact of viscosity cannot be compared. As both the models are based on the model in [14] in terms of the radial component, the observations are alike. Graphs given in the following paragraph reveal more.

Figure 5:

(a) The radial distribution of pressure difference p(r, z) − p(0, z) based on (36). (b) The axial distribution of pressure difference p(r, z) − p(r, 0) based on (37) with S = 0.98 and (c) the axial distribution of pressure difference p(r, z) − p(r, 0) based on (37) with S = 0.15. The other parameters assumed here for the diagrams are h = 50 and a = 1.

Taking care of axial dependence of the azimuthal component of velocity, we draw plots showing its radial and axial profiles (Fig. 3). This component is strongest at the core, and inside and outside the core, it diminishes. This is revealed in the plots drawn. The radial profile has resemblance with simulated models of Tang et al. [41], and the vertical profile is parabolic even for viscous flows, obvious from the constructed model given by (28) for high Reynolds numbers.

In order to examine the impact of viscosity, we further draw some plots for different values of Reynolds number. We fix S = 0.98; i.e. the azimuthal and radial velocities are almost the same at their maxima. The radial length is varied in the range r = 0.5–3.0. We take the following three different cases and plot the axial profiles of the azimuthal component of velocity:

1. r = 1.0: This is the core where the velocity is maximum. Plots are drawn for Re = 5, 10, 100, and 10,000. The plots, drawn in Figure 4, show that the larger the Reynolds number, the lesser is the velocity. However, once Re = 100, further increase in Reynolds number has insignificant impact as we observe that curves corresponding to Re = 100–10,000 almost coincide.

2. r = 0.5: This is a point inside the core and plots are drawn for Re = 10, 100, and 10,000. The trends are reversed. The larger the Reynolds number, the lesser is the velocity.

3. r = 3.0: This is a point outside the core, and plots are drawn for Re = 10, 100, and 10,000. Beyond the core, the azimuthal velocity decreases. The trends for the Reynolds number are similar to that inside the core; i.e. the larger the Reynolds number, the lesser is the velocity.

7.4 Pressure Distribution

Radial pressure distributions for different axial positions are shown in Figure 5a. We set z = h/2a, h/3a, and h/4a and Re = 10,000. Profiles are similar to the theoretical observation of Arsen’yev et al. [42] and the numerically simulated observation of Tang et al. [41]. That is, as we move outward from the axis, the pressure increases, and also that pressure decreases as height increases. The drop from the circumferential pressure to the pressure at the axis, based on (36), increases with height for swirl ratio S = 0.98, or, in other words, the pressure at the axis, i.e. p(0, z), is minimum. Plots are differently drawn in Figure 5b, in which the axial distribution of pressure, based on (37), has been plotted for different radial distances. The pressure is observed to fall as the height is scaled; in other words, i.e. p(r, 0) is maximum. If we fix z, trends are similar to what is observed in Figure 5a. Further, pressure drop becomes more in magnitude with height. This is a combined effect of the two observations.

However, for S = 0.15 (i.e. the value Wang et al. [43] talked about), the trends for pressure drop, shown in Figure 5c, are similar to what Wang et al. [43] observed in an experimental investigation; i.e. for some height, pressure decreases (i.e. pressure drop increases in magnitude), while above that the trend gets reversed. However, it needs to be noted that neither we talk of roughness nor has it anything to do with roughness of the surface; this is just based on the model we have formulated. This is an indication that the swirl ratio plays a big role. Tang et al. [41] have discussed this in detail.

In order to examine the effect of viscosity, we study pressure against different Reynolds numbers. The plots are given in Figure 6. Pressure decreases with rising Reynolds number uniformly for all radial distances. This is an indication that quantitative difference in pressure is large between viscous and inviscid flows. Thus, an approximation, from viscous to inviscid flow, does not look fair; it will be difficult to fit experimental data with theoretical inviscid models.

Figure 6:

(a) The radial distribution of pressure difference p(r, z) − p(0, z) based on (36) with z = h/4a. (b) The axial distribution of pressure difference p(r, z) − p(r, 0) based on (37). The other parameters assumed here for the diagrams are S = 0.50, h = 50, a = 1, and r = 1.

8 Conclusions

New results are mainly related to the axial dependence of the various components of velocity. Although all the three components are functions of radial and axial coordinates, viscosity affects the azimuthal velocity and the pressure.

It is observed that the magnitude of the radial velocity increases to the maximum at the core but reverses the trend beyond and vanishes as it reaches the centerline. The magnitude reduces linearly with axial distance as per the supposition.

At the core, the larger the Reynolds number, the lesser is the velocity for moderate Reynolds number. For larger Reynolds number, insignificant impact is observed. However, inside and outside the core, the trends are reversed; i.e. the larger the Reynolds number, the lesser the velocity.

Radial pressure distributions for different axial positions are similar to theoretical, experimental, and numerical observations as given in the discussion. As we move outward from the axis, pressure increases, but pressure decreases with height. Further, drop of pressure from the circumference to the axis increases in magnitude with height.

Pressure decreases with rising Reynolds number uniformly for all radial distances. This is an indication that quantitative difference in pressure is large between viscous and inviscid flows.