Axisymmetric Tornado Simulations with a Semi-Slip Boundary

Abstract

The structure of natural tornadoes and simulated analogs are sensitive to the lower boundary condition for friction. Three-dimensional numerical simulations of storms require a choice for turbulence parameterizations and resolution of wind near the lower boundary. This article explores some of the consequences of choices of a surface drag coefficient on the structure of a mature simulated tornado, using a conventional axisymmetric model. The surface drag parameterization is explored over the range of the semi-slip condition, including the extremes of no-slip and free-slip. A moderate semi-slip condition allows for an extreme pressure deficit, but without the unrealistic vortex breakdown of the no-slip condition.

1. Introduction

In [1], we find this statement about including surface drag in a storm model: “Our philosophy is that the inclusion of modest surface drag represents a more physically consistent bottom boundary condition for tornadogenesis than does the habitually employed free-slip assumption.” Several other recent storm modeling efforts have investigated surface drag in storm evolution and tornadogenesis [2,3,4]. This article attempts something simpler: documenting the effect of surface roughness on the mature tornadoes that can occur in storm models.

In [5] we read: “Thus, surface friction paradoxically instigates the extreme pressure deficits and upward and rotary winds in end-wall tornadoes.” This paradoxical effect has long been recognized and studied in models. In models, “extreme” instances are transient events, for example the three-dimensional suction vortices of [6], or the axisymmetric super-critical events of [7]. Axisymmetric transient events did occur in the model simulations here, but are not documented. Also, presumably, repetitive and longer-lasting non-axisymmetric suction vortex events could occur in similar simulations if a non-axisymmetric model were employed, as in [8]. However, the purpose of this article is to help diagnose the effect of semi-slip boundary conditions on larger, average features of a mature tornado, rather than extreme events. The axisymmetric simulations are similar to the two in [9], though here we investigate a range of semi-slip conditions. We also present the model results in a dimensionless form, to allow for wider generalization.

As used here, semi-slip refers to the simplest application that we see atmospheric modeling, such as in [1,3,10]. It could also be called bulk drag. There is a more sophisticated method referred to as partial-slip, as in [11], discussion of which is beyond the scope of this work. The intent here is not to offer the best model for nature, but to offer predictions of what might be produced in current three-dimensional storm models that use semi-slip, or equivalently bulk drag.

The model is the axisymmetric Fiedler Chamber [12]. This article documents and extends the simulations that were published in [12]. The Python source code (which uses numpy) is freely available to completely reproduce the simulations shown here, as well as easily investigate transient events. Similarly, ref. [13] has extended the axisymmetric simulations in [12], going to more than 10 times higher Reynolds number (which is not done here) but not investigating semi-slip boundary conditions.

2. The Axisymmetric Model

The model domain extends between $0\le r\le 2$ and $0\le z\le 1$. Acting by itself, the stationary buoyancy field $b(r,z)$ could accelerate fluid to a speed of 1 on the axis $r=0$ (with parcels starting from rest). The chamber, including all boundaries, rotates about $r=0$ at rate $\Omega$.

The radial, azimuthal, and vertical momentum equations are:

$$\frac{\mathrm{d}u}{\mathrm{d}t}=-\frac{\partial \varphi }{\partial r}+\frac{v^2}{r}+2\Omega r+\frac{1}{R{e}_F}\left[\frac{\partial }{\partial r}\left(\frac{1}{r}\frac{\partial ru}{\partial r}\right)+\frac{{\partial }^{2}u}{\partial z^2}\right]$$

(1)

$$\frac{\mathrm{d}v}{\mathrm{d}t}=-\frac{uv}{r}-2\Omega r+\frac{1}{R{e}_F}\left[\frac{\partial }{\partial r}\left(\frac{1}{r}\frac{\partial rv}{\partial r}\right)+\frac{{\partial }^{2}v}{\partial z^2}\right]$$

(2)

$$\frac{\mathrm{d}w}{\mathrm{d}t}=-\frac{\partial \varphi }{\partial z}+b+\frac{1}{R{e}_F}\left[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial w}{\partial r}\right)+\frac{{\partial }^{2}w}{\partial z^2}\right]$$

(3)

As in [12], $R{e}_F^{-1}$ is the dimensionless viscosity and $\varphi$ is pressure divided by a constant density. The model is Boussinesq.

Steady solutions are facilitated by increased damping in the top half of the domain:

$$\frac{1}{R{e}_F}=max\left({10}^{-4},2\times {10}^{-3}(z-0.5)\right)$$

(4)

For example, if we interpret the model as representing a strong thunderstorm with domain height $h={10}^4\mathrm{m}$ and buoyancy that could potentially produce an updraft of $W=100\mathrm{m}·{\mathrm{s}}^{-1}$, then using a value of $R{e}_F^{-1}={10}^{-4}$ implies an eddy diffusivity of $K=100{\mathrm{m}}^2·{\mathrm{s}}^{-1}$.

The continuity equation is:

$$\frac{1}{r}\frac{\partial ru}{\partial r}+\frac{\partial w}{\partial z}=0$$

(5)

The boundaries at $z=1$ and $r=2$ are no-slip. A semi-slip boundary is applied at $z=0$:

$$\frac{1}{R{e}_F}\frac{\partial u}{\partial z}=\gamma u\frac{1}{R{e}_F}\frac{\partial v}{\partial z}=\gamma v$$

(6)

These semi-slip conditions arise from the dimensional equation for a drag force , such as in [3,11]:

$$K\frac{\partial u^{†}}{\partial z^{†}}=C_D{V}_H{u}^{†}$$

(7)

where $u^{†}$ and $z^{†}$ are the dimensional horizontal velocity and vertical coordinate. In terms of our dimensionless u and z, (7) is

$$\frac{KW}{h}\frac{\partial u}{\partial z}=C_D{V}_{H}Wu$$

(8)

Dividing by $W^2$ gives:

$$\frac{K}{Wh}\frac{\partial u}{\partial z}=C_D\frac{V_H}{W}u$$

(9)

and (6) results with

$$\gamma =C_D\frac{V_H}{W}$$

(10)

A surface with identical roughness everywhere presents a constant value of $C_D$. Motivated by simplicity, various constant values of $\gamma$ are specified in the simulations here. In hindsight, a more applicable set of experiments would have allowed $\gamma$ to vary with $V_H$. Nevertheless, even in laminar flow, (6) with constant $\gamma$ has been applied when a thin film of reduced viscosity is adjacent to a wall [14]. So it is not impossible that somebody could reproduce the constant-$\gamma$ simulations here in a laboratory apparatus.

If we have a storm model employing (7) with a grassland value $C_D=0.01$ [3], a surface wind speed of order $V_H=30\mathrm{m}·{\mathrm{s}}^{-1}$, and h and K as above, then we would employ $\gamma =0.003$ in this axisymmetric model to anticipate the structures we would see in the storm model.

For diagnostic purposes, we plot a quantity related to head H of [12]:

$$H=\frac{1}{2}\left(u^2+v^2+w^2\right)+\varphi$$

(11)

In steady flow and in the absence of significant viscous and buoyancy effects, H is conserved along streamlines and is useful in predicting properties of the vortex corner region [12] .

In this article, we plot a modified head $H^{\prime }$:

$$H^{\prime }=H-q$$

(12)

where q is the kinetic energy contributed by buoyancy (but not pressure), which is easily prognosticated with:

$$\frac{\mathrm{d}q}{\mathrm{d}t}=wb$$

(13)

In principle, $H^{\prime }$ is conserved along streamlines of steady flow, even in the presence of b, but again without viscous effects. Changes in $H^{\prime }$ along streamlines are attributable to viscous effects. In the model calculations, the diagnostic q is reset to zero as fluid comes in through $r=0.5$, so any q that is accumulated, and used in the plot of $H^{\prime }$, is from one cycle through the central region of the chamber. By visualizing $H^{\prime }$, instead of H, dissipation in the core and updraft can be discerned. Near the lower boundary, H and $H^{\prime }$ will be nearly identical, and either can be used in the analysis of the corner flow region as in [12].

3. Results

Five values of $\gamma$ are considered, ranging from free-slip to no-slip: 0, $0.003$, $0.01$, $0.03$, ∞. The solutions are initialized as motionless. The solutions are integrated out to time $t=225$. Ultimately, the grid is $181\times 181$ grid points, with the grid doubled at $t=200$ and $t=220$. The solutions with $181\times 181$ grid points are nearly identical to those with $91\times 91$ grid points.

For each value of $\gamma$, at least three values of $\Omega$ are specified. The range of values of $\Omega$ are chosen so that a lower one allows for a vortex with updraft in the center of vortex to a larger one with a downdraft, and at least one central value showing a transition from an updraft at lower levels to a downdraft at upper levels. The exception to this is the free-slip case, $\gamma =0$, which does not display a solution with an updraft in the core. The free-slip solutions in Figure 1 and the no-slip solutions in Figure 2 have been analyzed in [12].

The solutions with $\gamma =0.03$ (Figure 3) are similar to the no-slip case $\gamma =\infty$ (Figure 2). The transitions for mid-range $\Omega =0.025$ and $\Omega =0.03$ are vortex breakdown events: the lower portion of the vortex has an updraft that prevents the centrifugal wave of enlargement from traveling down to the surface. This transition to the wider vortex aloft is analogous to a hydraulic jump. As in the no-slip case, the core of vortex, below the transition, is fluid that has lost circulation by friction but has lost little $H^{\prime }$ (or H) in the process of traveling into the core of the vortex. The size of the simulated vortices, particularly those with an updraft in the core, are much too large, as compared with the buoyant “storm” to be credible analogs to tornadoes. Giant axisymmetric vortex breakdown events are not reported to have been observed in tornadoes. More likely, the supercritical end-wall vortices exist as non-axisymmetric, and possibly multiple, suction-vortex events near the surface [6,8]. Though the core radius could be diminished by working with higher values of $R{e}_F$ and greater resolution, such simulations stray from useful application to storm models, where a value of $R{e}_F$ based on molecular viscosity is not relevant.

Moving on to less friction, we see that for $\gamma =0.01$ (Figure 4) that a supercritical vortex is still possible, with an updraft in the core. The size of the vortex, relative to the buoyant “storm” is smaller and more realistic than those for $\gamma =\infty$ or $\gamma =0.03$. Note, unlike the simulations with larger values of $\gamma$, at $r=0$ the updraft is slightly diminished from that in the surrounding core. This may seem paradoxical that this is happening with less friction, but the paradox perhaps could be resolved with an understanding of the boundary layer of a potential vortex [12]. As $\Omega$ is increased, a downdraft descends to the surface. Again, the core size is more realistic, as compared with those for larger $\gamma$.

For $\gamma =0.003$ (Figure 5), solutions are naturally tending towards those of $\gamma =0$. Some differences are: greater amplification of v near the surface, a pinch in the extent of the pressure deficit at the surface, and a stronger downdraft and updraft near the surface. In these, and in all simulations, we note the strongest winds near the surface are just outside a core of reduced $H^{\prime }$. In comparison with the free-slip case of Figure 1, we can see that some, but not all, of the reduction in $H^{\prime }$ is due to friction at the surface.

4. Conclusions

We anticipate that a storm model with $C_D=0.01$ as in [3] will produce tornadoes like those here with $\gamma =0.003$, which would have properties close to being free-slip. A downdraft will extend to the surface, and the cross-section in [3] shows exactly that. A value of $C_D=0.0014$ as in [1] would tend even more so to free-slip properties. However, the judgements being express here are based on the diffusion coefficient $R{e}_F^{-1}$ used here.

The laminar, axisymmetric solutions we present with $0.003\le \gamma \le 0.03$ are akin to the azimuthally averaged solutions of a three-dimensional large-eddy model [8]. The model here is computationally cheaper and simpler to analyze than that of [8], but can investigate only a more limited suite of questions. Insofar as this model with a semi-slip boundary is consistent with features in [8] and real tornadoes, it could be used to anticipate the structure and intensity of the lower region of the tornado and its breakup into multiple suction vortices.