Numerical analysis of fluid-body interaction considering vortex and tornado-like flows

Abstract

A numerical analysis of fluid-body interaction is performed in this work in order to evaluate the influence of vortex and tornado-like flows on immersed objects. Velocity profile models are adopted to generate vortical flow fields based on time-dependent boundary conditions and a finite element formulation is used for spatial discretization, where eight-node hexahedral elements with one-point integration are adopted. In addition, an arbitrary Lagrangian-Eulerian (ALE) approach is proposed to describe the relative motion between vortex flow and immersed objects. The flow governing equations are discretized using an explicit two-step Taylor-Galerkin scheme and tornado flow fields are simulated using the Rankine Combined Vortex Model (RCVM) and the Vatistas Model. Turbulence modeling is performed using Large Eddy Simulation (LES) with the Smagorinsky’s sub-grid scale model. Problems involving moving and stationary tornadoes interacting with fixed and moving objects are analyzed, where significant aerodynamic forces are observed on the immersed bodies, producing also significant changes in the vortex flow characteristics.

Keywords:
Computational Wind Engineering (CWE); Finite Element Method (FEM); Large Eddy Simulation (LES); Tornado Flows; Vortex Profile Models

Graphical Abstract

1 INTRODUCTION

The problem of body-vortex interaction is very important for Wind Engineering applications, where structures are subjected to vortical flows induced by upstream obstacles or special atmospheric events such as tornadoes and downbursts. Considering the flow complexity and many other aspects involved in this problem, additional investigations are required in order to further understand the effects of an impinging vortex impacting on immersed objects. In this sense, advances observed in computer technology and numerical techniques have contributed to the development of reliable numerical models for flow simulation and fluid-body interaction.

Numerical models based on vortex velocity profiles take into account vortex motion and immersed bodies in a simple manner, although all the possible flow patterns observed in a tornado cannot be represented completely. Mathematical expressions for vortex models can be obtained by solving the Navier-Stokes equations in cylindrical coordinates for particular situations (Lewellen, 1976). The Rankine Combined Vortex Model (RCVM) (Rankine, 1891), for instance, is deduced by considering an ideal fluid with a velocity field composed of two regions: the forced vortex region in the inner core and the external free vortex region, where only the tangential velocity is not null. In the Burgers-Rott Model, the three velocity components and pressure are considered to form a unicellular vortex (Burgers, 1948; Rott, 1958). A similar approach is also considered in the Sullivan Model (Sullivan, 1959), but forming a two-cell vortex. In recent years, algebraic velocity profiles have been developed, which try to approximate all these models simultaneously using a single model (Wood and White, 2011).

A numerical scheme based on the RCVM was first implemented by Wilson (1977) in a finite difference code using a system of two-dimensional Euler equations to determine loads induced by tornado flows on bodies with square and rectangular shapes. The RCVM was first utilized in a finite element simulation by Selvam (1985), where a two-dimensional tornado model was adopted. Later, Selvam (1993) applied his RCVM formulation to a three-dimensional analysis considering a logarithmic vertical profile, where the effect of viscosity was considered and turbulence was simulated using RANS and the k-ε turbulence model (see also, Kareem and Kijewski, 1996). No difficulty with respect to the imposition of boundary conditions was reported. The same profile model was employed by Selvam and Millett (2003) and Selvam and Millett (2005) using a three-dimensional model and LES for turbulence modeling. A sufficiently refined mesh was utilized and different angles of attack of the tornado trajectory line with respect to the position of a building model immersed in the flow were considered. Strasser and Selvam (2015b) used the Algebraic Model proposed by Vatistas et al. (1991) to simulate the vortex-cylinder interaction in a two-dimensional model by varying their relative size and the impact time. A detailed comparison of vortex velocity models may be found in Alrasheedi (2012) and Strasser and Selvam (2015a).

A numerical investigation considering vortex and tornado-like flows based on velocity profile models is performed in this work in order to evaluate the aerodynamic loads produced by tornado flows on immersed structures. In addition, an arbitrary Lagrangian-Eulerian (ALE) approach is proposed to simulate the effects of relative movement between vortices and objects. The numerical model adopted here is built using the explicit two-step Taylor-Galerkin scheme in the context of the Finite Element Method (FEM), where hexahedral elements with one-point integration and hourglass control are utilized for spatial discretization. Turbulence is considered using Large Eddy Simulation (LES) and the Smagorinsky’s sub-grid scale model. The flow is supposed to be incompressible, although the pseudo-compressibility hypothesis is adopted, and an isothermal process is assumed. A mesh movement scheme is utilized to simulate relative motions between vortices and immersed structures, independently if the moving object is identified with the immersed structure or the flow vortex. Vortex and tornado-like flows are simulated using the Modified Rankine Combined Vortex Model (MRCVM) formulation and the algebraic model proposed by Vatistas et al. (1991), in which a translation velocity is incorporated in order to allow the vortex structure to follow any direction in the computational space. The different methodologies utilized here are evaluated and tested comparatively using a classical application involving a circular cylinder submitted to two-dimensional flow. A three-dimensional realistic case is then analyzed using a logarithmic boundary layer function superimposed to the present vortex profile models to simulate the wind effects on a cubic building model subject to the flow generated by a moving tornado. Different tornado paths are considered to evaluate the aerodynamic forces induced on the building surface due to the translating vortex flows simulated in this work.

2 FLOW FUNDAMENTAL EQUATIONS

In the field of Computational Wind Engineering (CWE), wind flows are generally assumed to be incompressible, turbulent and isothermal, where air is considered as a Newtonian fluid with no gravity effects (see Braun and Awruch, 2009). In this work, an arbitrary Lagrangian-Eulerian (ALE) formulation is adopted for the kinematic description of the flow field using a Cartesian coordinate system, which leads to a system of equations given by:

Navier-Stokes equations:

$$\frac{\partial v_i}{\partial t}+\left(v_j-w_j\right)\frac{\partial v_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial p}{\partial x_j}{\delta }_{ij}+\frac{\partial }{\partial x_j}\left[\frac{\mu +{\mu }_t}{\rho }\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_{\mathrm{j}}}{\partial x_i}\right)+\frac{\lambda }{\rho }\frac{\partial v_k}{\partial x_k}{\delta }_{ij}\right]\left(\mathrm{i},j,k=1,2,3\right)in{\Omega }^{\text{f}}$$ (1)

Mass conservation equation using pseudo-compressibility (Chorin, 1967):

$$\frac{\partial p}{\partial t}+\rho {\beta }^2\frac{\partial v_j}{\partial x_j}=0(\text{j}=1,2,3)in{\Omega }^{\text{f}}$$ (2)

where the flow variables are the components $v_i$ of the flow velocity vector v and the thermodynamic pressure $p$, which are described as functions of the Cartesian coordinate vector x and time $t$ within the fluid spatial domain ${\Omega }^{\text{f}}$ and during the time interval $[t_0,t_f]$. For dynamic finite element grids and turbulent flows, a mesh velocity vector with components $w_i$ and the eddy viscosity ${\mu }_t$ must be defined using special numerical approaches described later. The flow velocity components $v_i$ and the mesh velocity components $w_i$ are expressed according to the directions of the orthogonal Cartesian axes $x_i$. The physical properties that characterize the fluid are the dynamic and volumetric viscosities, $\mu$ and $\lambda$, and the specific mass $\rho$, considering that $\beta$ is the pseudo-compressibility parameter, which is sometimes associated with the speed of sound. The components of the Kronecker delta are denoted by ${\delta }_{ij}$ (${\delta }_{ij}$ = 1 if $i=j$; ${\delta }_{ij}$ = 0 if $i\ne j$).

Initial and boundary conditions are given by:

$$v_{\text{i}}\left[x_1,x_2,x_3\right]=v_{\text{i}}^{\text{0}}\left(\text{i}=\mathrm{1,2,3}\right);p\left[x_1,x_2,x_3\right]=p^{\text{0}}$$ (3)

$$v_{\text{i}}=v_{\text{i}}^{\text{*}}\text{on}{\Gamma }_v;p=p^{\text{*}}\text{on}{\Gamma }_p$$ (4)

$s_i=s_{\text{i}}^{\text{*}}=\left[-p{\delta }_{\text{ij}}+\left(\mu +{\mu }_t\right)\left(\frac{\partial v_{\text{i}}}{\partial x_{\text{j}}}+\frac{\partial v_{\text{j}}}{\partial x_{\text{i}}}\right)+\lambda \frac{\partial v_{\text{k}}}{\partial x_{\text{k}}}{\delta }_{\text{ij}}\right]n_{\text{j}}\text{on}{\Gamma }_{\sigma }\left(\mathrm{i},j,k=\mathrm{1,2,3}\right)$ (5)

where $v_{\text{i}}^{\text{0}}$ and $p^{\text{0}}$ are initial values for the velocity and pressure fields, $v_{\text{i}}^{\text{*}}$ and $p^{*}$ are the flow velocity components and pressure values prescribed on boundaries ${\Gamma }_v$ and ${\Gamma }_p$ of the spatial domain ${\Omega }^{\text{f}}$, respectively, $s_{\text{i}}^{\text{*}}$ are components of the flow traction vector prescribed on boundary ${\Gamma }_{\sigma }$ according to the Cartesian directions $x_{\text{i}}$ and $n_{\text{j}}$ are components of the unit normal vector at a point on boundary ${\Gamma }_{\sigma }$ according to the Cartesian directions $x_{\text{j}}$.

2.1 Boundary conditions for generation of vortex and tornado-like flows

According to Strasser and Selvam (2015a), the cross-section of a tornado vortex may be considered as a composition of three different regions: (1) internal laminar core, (2) transition region and (3) external turbulent region, where the flow velocity field may be decomposed into translation and tangential velocity components. Figure 1 shows some typical tangential velocity profiles utilized by different authors (see Kim and Matsui, 2017; Strasser, 2015; Strasser and Selvam, 2015a), which indicates that the tangential velocity values increase with the distance r to the vortex center where a maximum value is obtained (V_θ,max) at the critical radius $r_c$. For $r>r_c$ the tangential velocity values decrease as the distance to the vortex center is increased.

Figure 1
Analytical models for the tangential velocity profile: (a) schematic view of the vortex regions; (b) distribution of the normalized tangential velocity ($V_{\theta }/V_{\theta ,\max }$) over the relative radius ($r^{\prime }/r_c$).

The vortex and tornado-like flow fields are generated here using transient boundary conditions according to the velocity profile models employed in this work, which are defined considering two coordinate systems (see Fig. 2): (a) a fixed coordinate system XYZ, which is defined arbitrarily, and (b) a moving coordinate system X’Y’Z’, which is associated with the tornado vortex and with its origin defined at the vortex center. It is important to notice that $t^{*}$ = 0 when the tornado vortex center is located at the closest point of the tornado trajectory with respect to the origin of the fixed coordinate system XYZ.

Figure 2
Coordinate systems utilized in the definition of transient boundary conditions for vortex and tornado-like flow fields.

The position of a boundary node with respect to the vortex center is determined by the radial distance using the following expression:

$${r^{\prime }}^2={x^{\prime }}_1^2+{x^{\prime }}_2^2={\left(x_1-l_x-U_x{t}^{*}\right)}^2+{\left(x_2-l_y-U_y{t}^{*}\right)}^2\text{with}{x^{\prime }}_3=x_3=0$$ (6)

where ${x^{\prime }}_1$, ${x^{\prime }}_2$ and ${x^{\prime }}_3$ are the coordinates of a boundary node in the moving coordinate system and $x_1$, $x_2$ and $x_3$ are the corresponding coordinates in the fixed coordinate system, $U_x=U_{\infty }.\cos \varphi$ and $U_y=U_{\infty }.\sin \varphi$, where $\varphi$ is the angle between the vortex path and the longitudinal direction $x_1$ and $U_{\infty }$ is the free flow speed, which is also associated with the module of the vortex translation vector U = ($U_x$,$U_y$), $l_x$ and $l_y$ are coordinates of the closest point of the tornado trajectory with respect to the origin of the fixed coordinate sytem XYZ and $t^{*}$ = $t$ - $T_{lag}$ (see Fig. 2).

Several vortex models are available in the literature to describe the tangential velocity profile, which are classified by Strasser (2015) using the following classes: bi-regional profiles, continuous profiles and algebraic profiles. Special attention must be paid to the algebraic profiles, considering that this class of velocity profile is able to reproduce all profiles related to the remaining classes. In the present work, the following tangential velocity profiles are utilized:

• The Modified Rankine Combined Vortex Model - MRCVM (Hughes, 1952):

  $$V_{\theta }^{MRCVM}\left(r^{\prime }\right)=\omega \cdot r^{\prime }0\le r^{\prime }/r_c\le 1.0$$(7)
  $$V_{\theta }^{MRCVM}\left(r^{\prime }\right)=\omega \cdot r_c\cdot {\left(\frac{r_c}{r^{\prime }}\right)}^x{r}^{\prime }/r_c>1.0$$(8)

• The Vatistas Model (Vatistas et al., 1991):

  $$V_{\theta }^{VM}\left(r^{\prime }\right)=\omega \cdot r^{\prime }\cdot r_c{}^2\cdot {\left(\frac{2}{\left({r^{\prime }}^{2n}+r_c{}^{2n}\right)}\right)}^{1/n}$$(9)

where the Lamb-Oseen/Burgers-Rott (L-O/B-R) tangential velocity profile is obtained when $n=2$, and the MRCVM for $n=100$ and $x=1$. The Vatistas model is able to obtain more realistic vortices and to reduce the numerical error associated with the discontinuity observed in the tangential velocity profile at $r^{\prime }=r_c$ for RCVM models (Strasser and Selvam, 2015b). The maximum tangential velocity is obtained at $r_c$, where $V_{\theta ,max}\left(r_c\right)=\omega r_c$, being $\omega$ the vortex angular velocity.

In order to obtain the flow velocity components in terms of global Cartesian coordinates, the following equations are employed:

$$v_1\left(t\right)=U_x-V_{\theta }\left(r^{\prime }\right)\frac{{x^{\prime }}_2}{r^{\prime }};v_2\left(t\right)=U_y+V_{\theta }\left(r^{\prime }\right)\frac{{x^{\prime }}_1}{r^{\prime }};v_3\left(t\right)=0$$ (10)

The velocity profile models utilized in the present numerical simulations are restricted to the RCVM and L-O/B-R models, where the distribution of the flow velocity components is defined as follows:

$$\left\{\begin{array}{l}{v}_1\left(t\right)=\left[U_{\infty }\cos \varphi +\frac{\left(U_{\infty }\sin \varphi \cdot t^{*}+l_y-x_2\right)}{r\text{'}}\cdot V_{\theta }\left(r\text{'}\right)\right]Z_f(x_3)\hfill \\ v_2\left(t\right)=\left[U_{\infty }\sin \varphi +\frac{\left(x_1-l_x-U_{\infty }\cos \varphi \cdot t^{*}\right)}{r\text{'}}\cdot V_{\theta }\left(r\text{'}\right)\right]Z_f(x_3)\hfill \\ v_3\left(t\right)=0\hfill \end{array}\right.$$ (11)

where $x_1$, $x_2$ and $x_3$ are coordinates of a point of the flow field with respect to the origin of the fixed coordinate system. In order to develop a three-dimensional flow field for the vortex models, a logarithmic vertical profile is adopted, i.e.:

$$Z_f(x_3)=\frac{u^{*}}{\kappa }\ln \left(\frac{x_3+z_0}{z_0}\right)$$ (12)

where $u^{*}$ is the friction velocity, which is defined considering that $Z_f=1$ at $x_3=h$, where $h$ is the immersed object height, $\kappa =0.4$ is the von Kármán constant and $z_0=0.00375m$ is the roughness length associated with the characteristics of the terrain surface utilized in this work.

2.2 Turbulence modeling

Large Eddy Simulation (LES) is adopted here to simulate turbulent flows, where a spatial filtering process is applied to the flow equations in order to decompose the flow scales into large and small scales. In this context, large scales are directly resolved considering the elements characteristic lengths of the finite element mesh and scales below the mesh resolution are modeled using sub-grid closure models. As a result of the spatial filtering procedure, unresolved terms are obtained, which are usually defined as components of the sub-grid Reynolds stress tensor, i.e.:

$${\overline{\tau }}_{ij}^{SGS}=-\rho \left(\overline{{v^{\prime }}_i{v^{\prime }}_j}\right)=2{\mu }_t{\overline{S}}_{ij}$$ (13)

where ${v^{\prime }}_i$ are the components of the sub-grid scale velocity vector v’ and ${\overline{S}}_{ij}$ are the components of the strain rate tensor described in terms of filtered velocity components ${\overline{v}}_i$, which are given by the following expression:

$${\overline{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\overline{v}}_i}{\partial x_j}+\frac{\partial {\overline{v}}_j}{\partial x_i}\right)$$ (14)

Notice that overbars indicate a filtered (or large scale) quantity.

In the present work, the eddy viscosity ${\mu }_t$ is determined using the classical Smagorinsky’s sub-grid model (Smagorinsky, 1963), which may be expressed as:

$${\mu }_t=\rho {\left(C_S\overline{\Delta }\right)}^2\left|\overline{S}\right|$$ (15)

where $C_S$ is the Smagorinsky constant, $\left|\overline{S}\right|$ is the norm of the filtered strain rate tensor and $\overline{\Delta }$ is the filter characteristic length, which is usually associated with the element volume (${\Omega }_e^{}$) in a finite element formulation, i.e. $\overline{\Delta }={\Omega }_e^{1/3}$.

3 NUMERICAL MODEL

The flow numerical analysis is performed in this work using the explicit two-step Taylor-Galerkin scheme (see, for instance, Donea, 1984; Kawahara and Hirano, 1983) in the context of the Finite Element Method, where eight-node hexahedral elements with one-point quadrature and hourglass control are utilized in order to avoid spurious modes. In the present scheme, the flow variables are approximated over time using second-order Taylor series, which are evaluated at $t^{n+1/2}$ and $t^n$. The time-discretized equations correspond to the form obtained by the well-known Lax-Wendroff method, which was used formerly in finite difference models. The flow variables at $t^{n+1/2}$ are evaluated in the first step using the following equations:

$$\begin{array}{l}{\overline{v}}_i^{n+1/2}=v_i^n+\frac{\Delta t}{2}\left\{-\left(v_j^n-w_j^n\right)\frac{\partial v_i^n}{\partial x_j}-\frac{1}{\rho }\right.\frac{\partial p^n}{\partial x_i}+\\ \left.\frac{\partial }{\partial x_j}\left[\left(\frac{\mu +{\mu }_t}{\rho }\right)\left(\frac{\partial v_i^n}{\partial x_j}+\frac{\partial v_j^n}{\partial x_i}\right)+\frac{\lambda }{\rho }\frac{\partial v_k^n}{\partial x_k}{\delta }_{ij}\right]+\frac{\Delta t}{4}\left(v_j^n-w_j^n\right)\left(v_k^n-w_k^n\right)\frac{{\partial }^2{v}_i^n}{\partial x_j\partial x_k}\right\}\end{array}$$ (16)

$$p^{n+1/2}=p^n+\frac{\Delta t}{2}\left(-\rho {\beta }^2\frac{\partial v_j^n}{\partial x_j}\right)$$ (17)

$$v_i^{n+1/2}={\overline{v}}_i^{n+1/2}-\frac{1}{\rho }\frac{\Delta t}{4}\frac{\partial }{\partial x_i}\left(p^{n+1/2}-p^n\right)$$ (18)

The flow variables at $t^{n+1}$ are finally obtained using:

$$v_i^{n+1}=v_i^n+\Delta v_i^{n+1/2}$$ (19)

$$p^{n+1}=p^n+\Delta p^{n+1/2}$$ (20)

where:

$$\Delta v_i^{n+1/2}=\Delta t{\left\{-\left(v_j-w_j\right)\frac{\partial v_i}{\partial x_j}-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\frac{\mu +{\mu }_t}{\rho }\right)\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\right)+\frac{\lambda }{\rho }\frac{\partial v_k}{\partial x_k}{\delta }_{ij}\right]\right\}}^{n+1/2}$$ (21)

$$\Delta p^{n+1/2}=\Delta t{\left\{-\rho {\beta }^2\frac{\partial v_i}{\partial x_i}\right\}}^{n+1/2}$$ (22)

Superscripts indicate the time where the corresponding flow variable is evaluated. The time step $\Delta t=t^{n+1}-t^n$ is obtained locally using the stability condition $\Delta t=\alpha .\Delta x/(c+U_{\infty })$, where $\Delta x$ is the characteristic length of the element, $U_{\infty }$ is the reference flow speed and $\alpha$ is a safety coefficient ($\alpha <1$), while $\beta$ is the pseudo-compressibility parameter. A unique time step $\Delta t$ is utilized throughout the fluid mesh, which corresponds to the minimum value obtained among all the elements.

The weak form of the finite element equations are finally obtained applying the Bubnov-Galerkin weighted residual method and the Green-Gauss theorem on Eqs. (16) to (22), which leads to the following system of matrix equations:

$${\overline{v}}_i^{n+1/2}=v_i^n-M_D^{-1}\frac{\Delta t}{2}\left[\left(A+B\right)v_i^n+D_{ij}{v}_j^n-G_i^T{p}^n-t_i^n\right]$$ (23)

$$p^{n+1/2}=M_D^{-1}{M}^{*}{p}^n-M_D^{-1}\frac{\Delta t}{2}\left(\rho {\beta }^2{G}_j{v}_j^n\right)$$ (24)

$$v_i^{n+1/2}={\overline{v}}_i^{n+1/2}-\frac{\Delta t}{4\rho }{M}_D^{-1}{G}_i^T\left(p^{n+1/2}-p^n\right)$$ (25)

$$v_i^{n+1}=v_i^n-M_D^{-1}\Delta t\left(A{v}_i^{n+1/2}+D_{ij}{v}_j^{n+1/2}-G_i^T{p}^{n+1/2}-t_i^{n+1/2}\right)$$ (26)

$$p^{n+1}=M_D^{-1}{M}^{*}{p}^n-M_D^{-1}\Delta t\left(\rho {\beta }^2{G}_j{v}_j^{n+1/2}\right)$$ (27)

where $p$ and $v_i^{}$ are the finite element vectors containing pressure and velocity components ($i=1,2,3$) evaluated at element nodes. The element matrices $M_D^{}$ and $M^{*}$ are the lumped mass matrix and the modified mass matrix utilized by Kawahara and Hirano (1983) to stabilize the pressure field, which is given by $M^{*}=e{M}_D^{}+(1-e)M$M^*, where $M$ is the consistent mass matrix and e is a selective lumping parameter ($0\le e\le 1$). The remaining element matrices and vectors are defined as follows: $A$ and $B$ are the advection and stabilization matrices, $D_{ij}$ are the diffusion matrices, $G_i$ are the gradient matrices and $t_i^{}$ are the traction vectors referring to the boundary terms ($i,j=1,2,3$). Hexahedral elements with trilinear interpolation functions are employed here for approximation of both, the velocity and pressure fields. By using one-point integration, the finite element matrices and vectors can be evaluated analytically. However, an hourglass control numerical technique to avoid spurious modes on diffusive terms must be utilized (see Christon, 1997). Additional details on the finite element formulation employed in this work may be found in Braun and Awruch (2008) and Braun and Awruch (2009).

The mesh velocity vector $w$ is determined arbitrarily according to the mesh motion scheme adopted. In this work, the mesh motion is determined considering the components of the mesh velocity vector obtained by the following equation:

$$w_k^i=\frac{{\displaystyle \sum _{j=1}^{N_B}{w}_k^j{\left(d_{ij}\right)}^{-n}}}{{\displaystyle \sum _{j=1}^{N_B}{\left(d_{ij}\right)}^{-n}}}(i=\mathrm{1,...,}{N}_I;k=\mathrm{1,2,3})$$ (28)

where $N_B$ and $N_I$ are the number of boundary ALE nodes and internal ALE nodes, respectively, $d_{ij}$ is the Euclidian distance between an internal ALE node ($i$) and a boundary ALE node ($j$), which may be located on the fluid-structure interface or on an external boundary of the moving region, and n is a user defined parameter to control mesh flexibility. In the present work, this parameter is set to $n=3$ for all simulations where the mesh motion scheme is employed.

4 NUMERICAL APPLICATIONS

4.1 Moving circular cylinder subject to stationary vortex flow

In this example a two-dimensional circular cylinder is submitted to the flow generated by a stationary vortex, where the cylinder translates in the longitudinal direction of the computational domain with its center aligned with the center of the vortex. The main objective here is to demonstrate that the present approach is equivalent (in terms of flow-induced forces) to consider a fixed cylinder submitted to a moving vortex, which rotates and translates simultaneously, as indicated in Fig. 3. In this case, the ALE kinematic description and a special scheme for mesh motions in the FEM context are proposed to simulate the present application numerically. The present simulations are performed taking into account idealized flow conditions, whereas no three-dimensional effects are considered and a two-dimensional LES-type approach is adopted. A similar analysis was carried out recently by Guo and Cao (2019), but using the Lattice-Boltzmann and Immersed Boundary methods. Although the LES methodology is an inherently three-dimensional approach, a two-dimensional analysis can also be adopted approximately in the case of flows with homogeneous longitudinal fluctuations (Bruno and Khris, 2003).

Figure 3
Equivalent approaches for the vortex flow analysis: (a) moving vortex with stationary cylinder; (b) stationary vortex with moving cylinder.

Flows with different conditions and different numerical approaches for turbulence modeling are simulated here, where two Reynolds numbers (Re = V_∞D/ν = 10³ and 3,9x10³; see Table 1 for flow reference parameters) are investigated using LES and DNS (direct numerical simulation). The computational domains and boundary conditions utilized in the present analyses are shown in Fig. 4, considering the following flow situations: (a) uniform flow with stationary cylinder, (b) moving (translation + rotation) vortex with stationary cylinder and (c) stationary vortex with moving (translation) cylinder. For uniform flow, velocity boundary conditions are imposed at the inflow and lateral walls of the computational domain using the following values: $v_1^{}=V_{\infty }$ and $v_2^{}=0$, where $V_{\infty }$ is the free-stream flow speed. At the outflow, a uniform pressure distribution is prescribed with$p=0$. When a vortex is present, the flow velocity components are redefined according to the velocity profile models utilized here to generate the vortex flow field. At the outflow, no pressure is prescribed in this case and the no-slip condition is enforced on the surface of the immersed cylinder. Therefore, when the immersed object is stationary, all velocity components are null and when the object is translating, $v_1^{}=U_x$ and $v_2^{}=0$ are imposed, where $U_x$ is the cylinder translation velocity. In order to reproduce a two-dimensional flow field, the third component of the flow velocity vector must be prescribed considering $v_3^{}=0$ for all nodes of the finite element mesh.

Table 1
Flow parameters and fluid constants utilized in the uniform flow analysis.

Undisturbed flow speed - uniform flow (V∞)	0.001 m/s
Cylinder diameter (D)	1 m
Specific mass (ρ)	1 kg/m3
Kinematic viscosity (ν)	ν = V∞D/Re
Volumetric viscosity (λ)	0

Figure 4
Computational domains and boundary conditions utilized for the moving cylinder analysis. Measurement points - relative coordinates (with respect to the cylinder center): Point 1 (1.5D, -0.5D, 0.0); Point 2 (1.5D, 0.0, 0.0); Point 3 (1.5D, 0.5D, 0.0). Length unit: [m].

Initial conditions are also specified taking into account the different flow conditions analyzed here. A constant pressure field $p=0$ is considered initially for all simulations carried out here, while the initial conditions for the flow velocity field are considered as follows: (a) for uniform flow, $v_1^{}=V_{\infty }$, $v_2^{}=0$ and $v_3^{}=0$; (b) for flow with moving vortex: $v_1^{}=v_{1,VPM}^{}(T_{lag})$, $v_2^{}=v_{2,VPM}^{}(T_{lag})$ and $v_3^{}=0$; (c) for flow with stationary vortex: $v_1^{}=v_{1,VPM}^{}(0)$, $v_2^{}=v_{2,VPM}^{}(0)$ and $v_3^{}=0$, where $v_{1,VPM}^{}(t)$ and $v_{2,VPM}^{}(t)$ are flow velocity components corresponding to a velocity field defined by the velocity profile models adopted here for vortex flow simulation, considering that $T_{lag}$ is a time parameter determined according to the translation velocity and initial distance between the center of the immersed cylinder and the vortex center. Flow and numerical parameters utilized in the present simulations are indicated in Tables 1 and 2, respectively, and constants employed for vortex flow modeling are presented in Table 3. Notice that the fluid kinematic viscosity is given as function of Reynolds number and flow reference velocity, which is defined as the undisturbed flow speed ($V_{\infty }$) for uniform flow analysis or translation velocity ($U_x$) for vortex-cylinder interaction analysis.

Table 2
Numerical parameters for flow analysis.

Pseudo-compressibility parameter (β)	0.17 m/s
Smagorinsky’s constant - LES (CS)	0.2
Safety coefficient (α)	0.44
Time increments (Δt)	M604x128: 3x10-2 s
	M1208x256: 1.5x10-2 s
	M1812x384: 1x10-2 s
Selective lumping parameter (e)	0.9
Simulation time interval - [t0, tf]	1617 s

Table 3
Parameters for vortex flow modeling.

Vortex critical radius (rc)	3 m
Vortex origin, vertical coordinate (ly)	0 m
Cylinder (vortex) translation velocity (Ux)	0.034 m/s
Time lags (Tlag,i)	0.0; 1029.41 s
Angular velocity (ω)	5.67x10-3; 1.13x10-2 rad/s

In order to verify the influence of the mesh quality on the numerical results, three finite element meshes with different refinement levels are employed here, which are identified as follows: (a) M604x128, with 90,752 elements; (b) M1208x256, with 363,008 elements; (c) M1812x384, with 816,768 elements. The smallest element lengths for these computational grids are 1.38x10^-2 m, 6.94x10^-3 m and 4.63x10^-3 m, respectively. Guo and Cao (2019) utilized a much more dense mesh (5x10⁶ cells) to fulfill discretization requirements associated with the Immersed Boundary method. The computational domains used in the present simulations were extended laterally in order to accommodate mesh motions prescribed by the mesh motion scheme utilized here in conjunction with the ALE formulation.

Initial results are obtained considering the immersed cylinder subject to uniform flow. Predictions referring to the Strouhal number S_t and force coefficients acting on the cylinder are presented in Table 4, where numerical results obtained by other authors using 2D models are also shown. The influence of mesh quality on the numerical predictions can be observed and a comparison between LES and DNS solutions indicates slight differences when force coefficients and Strouhal number are evaluated for a Reynolds number Re = 10³. On the other hand, when a higher Reynolds number is considered (Re = 3.9x10³), the influence of turbulence modeling on the numerical results becomes clear. One can see a good agreement among the predictions presented here with respect to the Strouhal number, although some differences are observed for the time-averaged force coefficient ${\overline{C}}_x$ and vertical force coefficient $C_y^{max}$ (maximum amplitude). Experimental predictions usually indicate a drag force coefficient and Strouhal number of 1.0 and 0.21, respectively, for the range of Reynolds numbers analyzed here (see, for instance, Schlichting and Gersten, 2016), where 3D effects are significant. In order to reproduce these effects accurately, a 3D numerical modeling is required.

Table 4
Stationary cylinder subject to uniform flow: Strouhal number and aerodynamic force coefficients for Re = 10³ and 3.9x10³.

Mesh configurations	Source	${\overline{C}}_x$	$C_y^{max}$	$S_t$
604x128	Present work - Re = 103 (LES)	1.520	1.174	0.214
1208x256	Present work - Re = 103 (LES)	1.536	1.335	0.226
1812x384	Present work - Re = 103 (LES)	1.546	1.404	0.226
1812x384	Present work - Re = 103 (DNS)	1.531	1.533	0.238
1812x384	Present work - Re = 3.9x103 (LES)	1.709	1.818	0.232
1812x384	Present work - Re = 3.9x103 (DNS)	1.445	2.312	0.134
4000x1250	Guo and Cao (2019) - Re = 103	1.250	-	0.238
200x200	Wanderley and Levi (2002) - Re = 103	1.500	1.370	0.235
180x160	Hernandez and Sphaier (1999) - Re = 103	1.600	1.700	0.225
116x136	Beaudan and Moin (1994) - Re = 3900	1.74	2.165	0.263
120x145	Rajani et al. (2016) - Re = 3900	1.66	-	0.274

Flow turbulence characteristics are evaluated here using time histories of flow variables obtained at three different points in the wake region for the mesh configuration M1812x384, as indicated in Fig. 4. This analysis is performed considering uniform flow with two distinct Reynolds numbers, Re = 10³ and 3.9x10³, where LES and DNS are utilized comparatively. Results are presented in Table 5 for normalized turbulence intensity and normalized Reynolds stress components, which are calculated as follows:

$${\overline{I}}_{vi}=\frac{1}{V_{\infty }}\sqrt{\frac{1}{N}{\displaystyle \sum _{j=1}^N{\left({v^{\prime }}_i\right)}_j^2}};{\overline{{\tau }^{\prime }}}_{ij}=\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left(\frac{{v^{\prime }}_i{v^{\prime }}_j}{V_{\infty }^2}\right)}_k}$$ (29)

where ${v^{\prime }}_i$ are fluctuating flow velocity components and $N$ is the number of samples over the time interval where the flow variables are collected.

Table 5
Flow turbulence characteristics for circular cylinder submitted to uniform flow.

Point 1	Re = 103 (Uniform)		Re = 3.9x103 (Uniform)
	DNS	LES	DNS	LES
Normalized turbulence intensity - ${\overline{I}}_{v1}$	0.469	0.476	0.654	0.472
Normalized turbulence intensity - ${\overline{I}}_{v2}$	0.564	0.550	0.659	0.645
Normalized Reynolds stress - ${\overline{{\tau }^{\prime }}}_{11}$	0.220	0.227	0.427	0.222
Normalized Reynolds stress - ${\overline{{\tau }^{\prime }}}_{12}$	0.121	0.137	0.228	0.138
Normalized Reynolds stress - ${\overline{{\tau }^{\prime }}}_{22}$	0.318	0.303	0.434	0.416
Point 2	Re = 103 (Uniform)		Re = 3.9x103 (Uniform)
	DNS	LES	DNS	LES
Normalized turbulence intensity - ${\overline{I}}_{v1}$	0.154	0.154	0.588	0.155
Normalized turbulence intensity - ${\overline{I}}_{v2}$	0.990	0.949	0.881	1.037
Normalized Reynolds stress - ${\overline{{\tau }^{\prime }}}_{11}$	0.024	0.024	0.346	0.024
Normalized Reynolds stress - ${\overline{{\tau }^{\prime }}}_{12}$	0.0019	0.0008	0.0056	0.0009
Normalized Reynolds stress - ${\overline{{\tau }^{\prime }}}_{22}$	0.981	0.900	0.776	1.076
Point 3	Re = 103 (Uniform)		Re = 3.9x103 (Uniform)
	DNS	LES	DNS	LES
Normalized turbulence intensity - ${\overline{I}}_{v1}$	0.471	0.473	0.627	0.470
Normalized turbulence intensity - ${\overline{I}}_{v2}$	0.561	0.551	0.609	0.627
Normalized Reynolds stress - ${\overline{{\tau }^{\prime }}}_{11}$	0.222	0.224	0.393	0.221
Normalized Reynolds stress - ${\overline{{\tau }^{\prime }}}_{12}$	-0.126	-0.133	-0.140	-0.133
Normalized Reynolds stress - ${\overline{{\tau }^{\prime }}}_{22}$	0.315	0.303	0.371	0.394

The predictions obtained in the present investigation corroborate observations reported previously, where one can see that LES and DNS lead to similar results for uniform flow and Re = 10³. On the other hand, differences are significant for a higher Reynolds number, as expected. It is also observed that the measurement points 1 and 3 obtained equivalent values for the turbulence parameters evaluated here, except for the normalized Reynolds stress component ${\overline{{\tau }^{\prime }}}_{12}$, where different signs are identified. Nevertheless, this difference is in accordance with the alternate vortex shedding phenomenon occurring on the cylinder surface. The longitudinal turbulence parameters obtained at Point 2 are relatively small when compared with the corresponding values evaluated at Points 1 and 3, although transversal components such ${\overline{I}}_{v2}$ and ${\overline{{\tau }^{\prime }}}_{22}$ show higher values at that point. Notice that reliable predictions require here a 3D numerical modeling of the flow field. Nevertheless, the present results are useful as a first approach to the actual problem, from which important observation can be obtained, at least qualitatively.

The influence of the relative motion between vortex and cylinder is evaluated using the RCVM and L-O/B-R velocity profile models, where two conditions are investigated comparatively: (a) stationary cylinder subject to translating vortex flow ($V_{\theta ,\max }/U_x=0.5$) and (b) moving cylinder subject to stationary vortex flow, both considering a Reynolds number Re = 10³ and LES. Figure 5 shows results referring to streamline fields, which correspond to time instants when the vortex core passes the frontal, central and posterior regions around the immersed cylinder. Notice that the streamline fields obtained for the different conditions analyzed in the present investigation are quite distinct. When the immersed cylinder is stationary and the vortex is translating, one can see that the streamlines are predominantly horizontal, with curvilinear lines in the near wake of the immersed object and in the vortex core. The streamline curvature is associated with the tangential velocity $V_{\theta }$, such that the maximum curvature is found next to regions where the maximum tangential velocity $V_{\theta ,\max }$ is identified. On the other hand, when the vortex is stationary and the cylinder is moving, it is observed that the streamlines are concentric as well as the vortex core and the wake are well defined. No significant differences can be identified when the streamline fields obtained with the velocity profile models utilized here are compared each other.

Figure 5
Streamline fields obtained for different time instants and different flow conditions: (a) stationary cylinder with translating vortex; (b) moving cylinder with stationary vortex.

A scalar velocity field near the immersed cylinder is presented in Fig. 6 considering different time instants and the different flow conditions proposed for the present analysis, while the Reynolds number is set to Re = 10³ and LES is adopted. When the flow field referring to the moving vortex is observed, one can see that upstream the immersed object, the upper region of the computational field shows velocity values smaller than those observed in the lower region, which is justified considering that the approaching vortex is rotating in the counterclockwise direction. As the moving vortex enters the wake region of the immersed object, the vortex street is deviated upwards. On the other hand, when the stationary vortex flow is considered, the vortex core is well identified, with flow velocity increasing along the radial direction. As the moving cylinder is approaching the stationary vortex, the vortex street in the wake is also deviated upwards, resembling the flow characteristics observed previously.

Figure 6
Flow scalar velocity fields obtained for different time instants and different flow conditions: (a) stationary cylinder with translating vortex; (b) moving cylinder with stationary vortex.

Figure 7 shows the vorticity fields obtained here according to the flow conditions proposed (LES and Re = 10³). It can be observed that the flow fields are identical for both the flow situations analyzed in the present investigation. While the vortex core is not interacting with the circular cylinder, the flow pattern in the wake is similar to that observed in uniform flow conditions, indicating the predominance of the translational flow component. As the moving vortex approaches the immersed object, the vortex shedding behind the cylinder is modified by the translating vortex flow, which is affected by this interaction. When the vortex core coincides with the center of the immersed cylinder or when it is localized in the wake of the immersed cylinder, the vortical structures behind the cylinder are displaced vertically owing to the rotation of the moving vortex in the counterclockwise direction. One can see that the vortex flow pattern is strongly affected in the wake owing to strong interactions occurring in that region with vortices shed from the cylinder surface. Notice that eddies are also observed in the vortex core and the deviation of the vortex street upwards is now clearly identified.

Figure 7
Vorticity fields obtained for different time instants and different flow conditions: (a) stationary cylinder with translating vortex; (b) moving cylinder with stationary vortex.

Time histories of the aerodynamic force coefficients $C_x$ and $C_y$ obtained in the present simulations are shown in Fig. 8, which correspond to the condition where a stationary cylinder is subject to a translating vortex flow with $V_{\theta ,\max }/U_x=0.5$. Predictions using computational meshes with different refinement levels are compared here considering Re = 10³ and LES, where the velocity profile models RCVM and L-O/B-R are utilized to generate the vortex flow field. The moving vortex center coincides with the cylinder center when X = 0, while positions with X < 0 and X > 0 indicate that the vortex is approaching and leaving the cylinder, respectively. A first comparison is carried out taking into account two different computational meshes with different refinements levels, M1208x256 and M1812x384. Figure 8a shows that mesh M1812x384 leads to force coefficients with peak values slightly higher than those obtained with mesh M1208x256. In addition, one can see that differences between results obtained with the different meshes utilized here are not significant initially, but they increase as the vortex travels through the computational domain. A second comparison is presented in Fig. 8b, where results referring to the velocity profile models adopted here are compared using mesh M1812x384. Notice that the L-O/B-R velocity profile model leads to horizontal force values higher than those predicted using the RCVM profile. On the other hand, similar results are obtained in terms of force coefficient C_y for both models investigated here, although a significant loss of synchronization can be observed in the respective time histories as the vortex moves along its trajectory.

Figure 8
Stationary cylinder subject to moving vortex flow with $V_{\theta ,\max }/U_x=0.5$: (a) force coefficients for different mesh refinements; (b) force coefficients using the RCVM and L-O/B-R profile models.

In order to demonstrate the equivalence between the flow approaches (stationary cylinder with translating vortex; moving cylinder with stationary vortex) proposed in this work to simulate the effects of vortex flows on immersed objects, the respective time histories of force coefficients obtained here are presented in Fig. 9, taking into account a Reynolds number Re = 10³ and LES. In this case, two velocity ratio are analyzed: $V_{\theta ,\max }/U_x=0.5$ and $V_{\theta ,\max }/U_x=1.0$. The RCVM and L-O/B-R profile models are utilized to generate the vortex flow field, where the computational grid corresponds to mesh M1812x384. One can see that the force coefficients values for both the flow conditions analyzed here are almost the same when $V_{\theta ,\max }/U_x=0.5$, independent of the profile model utilized. On the other hand, some differences in the peak values and loss of synchronization can be observed between predictions obtained with the different flow approaches proposed in this work when $V_{\theta ,\max }/U_x=1.0$, where the mesh refinement plays a more important role. Furthermore, the mesh motion scheme adopted in this work may lead to element elongation in some mesh regions while the immersed object is moving, which may also lead to accuracy reductions in the flow velocity field associated with the vortex. Nevertheless, the similarity between predictions obtained with the different flow approaches utilized here is still very good. One can see that the flow-induced forces on the immersed object are more significant within the position interval -10 < X < 10 for a higher velocity ratio, where the rotational component in the vortex flow is equivalent to the translational component. On the other hand, when the flow rotational component is smaller ($V_{\theta ,\max }/U_x=0.5$), changes in the magnitude of the force coefficients are less noticeable, indicating the importance of the vortex rotation intensity in the evaluation of aerodynamic forces. In addition, it is observed that when the moving vortex (cylinder) is approaching the stationary cylinder (vortex), the force coefficient C_y is increased owing to the counterclockwise rotation of the vortex flow, while C_y decreases when the vortex (cylinder) is moving away from the cylinder (vortex).

Figure 9
Cylinder interacting with vortex flow, force coefficients for different flow conditions: (a) RCVM profile model; (b) L-O/B-R profile model.

It is important to notice that three-dimensional turbulent flow structures are present in the wake of circular cylinders for flows with the Reynolds numbers simulated in this work. In this sense, although those 3D turbulent structures are reproduced here, a 2D approach is adopted to evaluate approximately the influence of turbulence effects on the numerical results. Figure 10 shows a comparison between LES and DNS predictions, taking into account instantaneous vorticity fields for Re = 10³ and 3.9x10³, which correspond to the time instant when the vortex is passing through the circular cylinder. One can observe that the vorticity field in the wake is better defined when LES is utilized, while the wake region obtained with DNS is modified, especially for Re = 3.9x10³.

Figure 10
Instantaneous vorticity fields obtained with LES and DNS: (a) Re = 10³; (b) Re = 3.9x10³.

In order to evaluate the influence of the Reynolds number on the force coefficients acting on a circular cylinder submitted to a moving vortex, LES and DNS simulations are carried out taking into account two different flow conditions: Re = 10³ and 3.9x10³, both using a fixed cylinder and considering a moving vortex flow described with L-O/B-R model, where $V_{\theta ,\max }/U_x=0.5$ and mesh configuration M1812x384 is adopted. Figure 11 shows results referring to the force coefficients C_x and C_y given as functions of the moving vortex position. One can see that the LES and DNS predictions are similar when a lower Reynolds number is adopted, but significant differences are observed as the Reynolds number is increased. It is observed that LES leads to better approximations, as demonstrated in Fig. 10, especially when the moving vortex is traveling over the wake region (X > 0), where flow turbulence is important. Notice that DNS results are close to LES predictions for Re = 10³ and X < 0, indicating that LES is needed to obtain accurately the flow effects on the immersed structure when the moving vortex is passing through the wake region. From the numerical results obtained with the computational mesh M1812x384, the following y⁺ values are calculated: for Re = 10³, y⁺ = 1.93 (DNS) and y⁺ = 1.45 (LES); for Re = 3.9x10³, y⁺ = 2.88 (DNS) and y⁺ = 1.55 (LES).

Figure 11
Influence of the Reynolds number on aerodynamic force coefficients: (a) Re = 10³; (b) Re = 3.9x10³.

4.2 Cubic building model subject to three-dimensional tornado-like vortex

In this application, a cubic building model is subject to a flow field generated by a three-dimensional tornado-like vortex translating along the longitudinal direction of the computational domain. Four different tornado paths are considered here to evaluate the aerodynamic forces induced on the building surface due translating vortex flow, considering three paths aligned longitudinally and an oblique trajectory. The RCVM and L-O/B-R velocity profile models are adopted to generate the tornado flow field.

Figure 12 shows the computational domain utilized in the present investigation and boundary conditions for the flow fundamental equations, where the tornado trajectories proposed here are also indicated. Notice that time-dependent velocity boundary conditions (see Eqs. 6 to 12) are applied on the superior and lateral walls of the computational domain as functions of the global coordinates in order to induce the three-dimensional tornado-like flow field internally. The non-slip boundary condition is imposed on the ground and building walls. Initial conditions for the flow velocity are also defined considering the velocity profile models adopted in this work, where a velocity distribution is specified over the computational domain taking into account the initial position of the tornado vortex. Pressure initial conditions are specified considering $p=0$ for all finite element nodes of the computational mesh. Flow and numerical parameters utilized in the present simulations are given in Tables 6 and 7, respectively, and constants employed for tornado-like flow modeling are presented in Table 8. The flow field is characterized considering a Reynolds number Re = V_∞D/ν = 5.5x10⁵, where LES and the classical Smagorinsky’s sub-grid scale model are utilized for turbulence modeling.

Figure 12
Cubic building subject to three-dimensional tornado-like flow: computational domain and boundary conditions. Measurement points coordinates: Point 1 (1.5, 0.0, 0.0075); Point 2 (1.5, 0.0, 0.33); Point 3 (1.5, 0.0, 0.67); Point 4 (1.5, 0.0, 1.0). Length unit: [m].

Table 6
Cubic building subject to three-dimensional tornado-like flow: flow parameters and fluid constants.

Reference flow velocity - (V∞ = Ux + Vθ )	5.5 m/s
Building height (h)	1 m
Specific mass (ρ)	1 kg/m3
Kinematic viscosity (ν)	1x10-5 m2/s
Volumetric viscosity (λ)	0

Table 7
Cubic building subject to three-dimensional tornado-like flow: numerical parameters.

Pseudo-compressibility parameter (β)	55 m/s
Smagorinsky’s constant - LES (CS)	0.1
Safety coefficient (α)	0.4
Time increment (Δt)	5x10-5 s
Selective lumping parameter (e)	0.9
Simulation time interval - [t0, tf]	40 s

Table 8
Cubic building subject to three-dimensional tornado-like flow: parameters for tornado-like flow modeling.

Tornado critical radius (rc)	3 m
Tornado origin, vertical coordinates (ly)	-3 m (path 1)
	0 m (path 2)
	3 m (path 3)
	0 m (path 4 - oblique)
Tornado translation velocity (Ux)	1 m/s
Time lag (Tlag)	20 s
Tornado angular velocity (ω = Vθ/rc)	1.5 rad/s

A mesh quality analysis is carried out initially in order to evaluate the accuracy of the numerical formulation proposed in this work with respect to the determination of aerodynamic force coefficients induced by the tornado flow on the immersed cubic model, where two mesh resolutions are adopted: mesh M1 and mesh M2. In the present analysis, the L-O/B-R velocity profile model is adopted for tornado flow generation considering a tornado path aligned with the longitudinal axis of the reference coordinate system ($x_1$-axis, see Fig. 12). A comparison regarding the mesh configurations employed here and the mesh configuration adopted by Alrasheedi (2012) is presented in Table 9. Details of the mesh configuration utilized in the present analyses are shown in Fig. 13, where a schematic view showing the disposition of finite elements next to the immersed object is also presented for mesh M1. It is important to notice that the disposition of finite elements along the longitudinal direction $x_1$ is also adopted along the transversal and vertical directions $x_2$ and $x_3$, respectively.

Table 9
Cubic building subject to three-dimensional tornado-like flow: mesh configurations.

Source	Elements on the cube surface	Elements in the computational domain	Smallest element length (next to the cubic building)	Total nodal points
Present work - M1	35x35x35	139x139x87	0.0075 h	1,684,340
Present work - M2	40x40x40	162x162x100	0.004 h	2,560,400
Alrasheedi (2012)	40x40x20	78x78x56	0.01 h	340,704

Figure 13
Cubic building subject to three-dimensional tornado-like flow: finite element mesh and grid disposition (for mesh M1) next to the immersed object.

Figure 14 shows the numerical results obtained in the mesh quality analysis taking into account the aerodynamic force coefficients $C_x$, $C_y$ and $C_z$ ($C_x=2{\displaystyle \sum s_1}/\rho V_{\infty }^{2}A$, $C_y=2{\displaystyle \sum s_2}/\rho V_{\infty }^{2}A$ and $C_z=2{\displaystyle \sum s_3}/\rho V_{\infty }^{2}A$, where $A$ is the area of the cube face; see Eq. 5 for definition of $s_1$ and $s_2$) induced by the flow on the cubic building. One can observe that the mesh resolutions utilized here lead to very similar results while the tornado vortex is moving in regions away from the building position. On the other hand, some differences are identified within the time interval [15, 25 s], when the tornado vortex is interacting with the cubic building. In order to quantify theses differences, time-mean results are presented comparatively in Table 10 considering the time interval [15, 25 s] and the mesh configurations utilized in this work. Notice that the different mesh configurations employed in this example obtained similar time-averaged values for most of the force coefficients evaluated, indicating that mesh M1 is able to provide time-mean results with reasonable accuracy, considering the present formulation taken solely.

Figure 14
Cubic building subject to three-dimensional tornado flow: mesh quality analysis based on aerodynamic force coefficients, tornado path $l_y=0$ with the L-O/B-R velocity profile model.

Table 10
Cubic building submitted to moving vortex flow: time-mean aerodynamic coefficients.

Source	${\overline{C}}_x$	${\overline{C}}_y$	${\overline{C}}_z$
Present work - M1	0.27	0.14	0.59
Present work - M2	0.32	0.21	0.58

An additional comparison is presented in Table 11 taking into account peak values of the aerodynamic force coefficients $C_x$, $C_y$ and $C_z$ and the maximum pressure coefficient $C_{p\max }$ (normalized with respect to the reference speed $V_{\infty }=U_x+V_{\theta }$). The present results are obtained considering different mesh configurations, the L-O/B-R velocity profile model for tornado flow generation and a tornado path aligned with the longitudinal axis of the reference coordinate system, which are compared with predictions obtained by other authors utilizing numerical and experimental techniques. The pressure coefficient is calculated in this work using the following expression: $C_p=2p/\rho U_x^2$, where the tornado translation velocity is adopted as reference speed.

Table 11
Cubic building submitted to moving vortex flow: peak aerodynamic coefficients.

	Vθ/Ux	Approach	Cx	Cy	Cz	Cpmax
Alrasheedi (2012)	4.5	Numerical	0.83	0.94	1.45	-1.80
Selvam and Millett (2003)	4.5	Numerical	0.82	1.36	1.81	-2.20
Sengupta et al. (2008)	15.9	Experimental	1.82*		1.22	-
Present work - M1	4.5	Numerical	1.05	1.32	1.31	-2.08
Present work - M2	4.5	Numerical	0.89	0.82	1.38	-2.18

• ^*
  (Cx²+Cy²)^1/2

Results obtained here with the mesh configuration M2 are in good agreement with the numerical predictions presented by Alrasheedi (2012) using a finite difference model and LES. On the other hand, the numerical predictions obtained by Selvam and Millett (2003), who also employed a finite difference model and LES, are better reproduced when the mesh configuration M1 is adopted. The force coefficients calculated with mesh M1 are also very similar to experimental results presented by Sengupta et al. (2008). One can see that the present results show a reasonable agreement with predictions obtained by other authors, although a large dispersion can be observed among the reference results utilized here owing to different flow conditions adopted in the different works. It is worth mention that the tornado flow field is very sensitive to the vortex swirl ratio (ratio between angular and radial momenta) and ground conditions (smooth or rough). For additional information on this issue, see Cao et al. (2018).

Notice that the mesh resolutions utilized in the present application are similar to that adopted in the moving cylinder analysis for a smaller Reynolds number but using a finer mesh. This can be explained considering that for immersed objects with sharp edges, such as a cubic building, the flow separation is well determined, occurring along these sharp edges. As a result, the evaluation of force coefficients tends to be independent of the Reynolds number for Re > 10³, as demonstrated experimentally by Chien et al. (1951). In this case, the flow pattern around the object is quite independent of viscous influence and dynamic similarity will be attained at Reynolds numbers even well below that which should be simulated.

Figure 15 shows some results related to instantaneous pressure fields obtained from the present simulation using the mesh configuration M1 at $t=20s$, which are defined considering a vertical plane specified by $x_2=0$. The effects produced in the surroundings of the immersed building by the flow induced by the different tornado paths proposed here are investigated using the L-O/B-R velocity profile model. A distinctive suction zone is developed inside the tornado funnel with pressure magnitudes significantly higher than those identified outside. Notice that the pressure fields obtained for the tornado trajectories defined by $l_y=-3.0m$ and $l_y=3.0m$ are not equivalent, although they are equidistant to the immersed building. This behavior can be explained considering that the tornado vortex is translating and rotating in the counterclockwise direction, which leads to different interaction conditions between the translating vortex flow and the cubic building. One can see that the vortex core is usually displaced towards the region where the total velocity (translation + rotation) due to the tornado flow is minimum and the vortex core diameter for $l_y=-3.0m$ is comparatively larger than that observed for $l_y=-3.0m$.

Figure 15
Cubic building subject to three-dimensional tornado flow, instantaneous pressure fields on plane $x_2=0$ at $t=20s$ obtained with L-O/B-R profile model: (a) tornado path $l_y=0$; (b) tornado path $l_y=-3.0m$; (c) tornado path $l_y=-3.0m$.

Dimensionless pressure fields ($C_p=2p/\rho V_{\infty }^2$; see Table 6 for flow reference parameters) are presented in Fig. 16 in order to identify the position and instantaneous configuration of the tornado vortex in the computational domain at $t=20s$. Notice that the flow reference velocity $V_{\infty }$ utilized here correspond to the sum of translational and angular components of the vortex flow. The pressure iso-surfaces for $C_p=-0.35$ are shown considering different tornado trajectories and different velocity profile models, according to Alrasheedi (2012). A first comparison is performed taking into account pressure fields corresponding to the tornado path $l_y=0$ and obtained using the RCVM and L-O/B-R profile models (see Figs. 16a and 16b). Notice that the vortex structure obtained with RCVM is clearly smaller than the vortex structure predicted with L-O/B-R, indicating that the former profile model experiences higher numerical dissipation due to the sharp discontinuity observed in the RCVM velocity profile formulation. A second comparison is carried out considering the tornado paths defined with $l_y=-3.0m$ and $l_y=3.0m$ and using the L-O/B-R profile model. One can see that the pressure field is more disturbed when the tornado path defined with $l_y=3.0m$ is utilized owing to the counterclockwise rotation assumed by the tornado vortex, which leads to maximum flow velocities acting on the surface of the immersed object in this case. A final comparison is performed using the profile models RCVM and L-O/B-R and an oblique trajectory specified by 45° (see Fig. 12), where higher dissipation can be identified again when the RCVM profile model is utilized.

Figure 16
Cubic building subject to three-dimensional tornado flow, dimensionless pressure fields at $t=20s$: (a) tornado path $l_y=0$ with RCVM; (b) tornado path $l_y=0$ with L-O/B-R; (c) tornado path $l_y=-3.0m$ with L-O/B-R; (d) tornado path $l_y=3.0m$ with L-O/B-R; (e) tornado path with 45° and RCVM; (f) tornado path with 45° and L-O/B-R.

Results referring to aerodynamic forces induced by the tornado flow on the immersed object are presented in Fig. 17, where time histories of force coefficients are shown comparatively according to the different profile models and flow conditions analyzed here. The vertical yellow lines indicate the time interval when the tornado vortex is passing above the building. A first comparison considers a tornado path defined with $l_y=0$ and the predictions obtained with profile models RCVM and L-O/B-R, where the force coefficients predicted with the RCVM profile model are clearly smaller than those obtained with L-O/B-R. This result indicates again that the RCVM profile model suffers higher numerical dissipation when compared with L-O/B-R predictions. One can also observe that all force components show significant variations in their coefficient values around the physical time $t=20s$, which corresponds to the presumed impact time between the tornado vortex and the immersed body. The influence of the tornado flow on the aerodynamic forces acting on the body is notably observed within the time interval $[\mathrm{10,}30s]$.

Figure 17
Cubic building subject to three-dimensional tornado flow, aerodynamic force coefficients: (a) tornado path $l_y=0$ with RCVM; (b) tornado path $l_y=0$ with L-O/B-R; (c) tornado path $l_y=-3.0m$ with L-O/B-R; (d) tornado path $l_y=3.0m$ with L-O/B-R; (e) tornado path with 45° and RCVM; (f) tornado path with 45° and L-O/B-R.

In a second comparison, results related to the tornado paths defined by $l_y=-3.0m$ and $l_y=3.0m$ are analyzed considering predictions obtained with L-O/B-R profile model. Notice that the $C_z$ values obtained with both tornado trajectories are similar, but the remaining coefficients present different behavior over time for the different tornado paths. This is explained taking into account that the tornado vortex is translating and rotating in the counterclockwise direction, which leads to greater drag forces on the immersed body when the tornado path with $l_y=3.0m$ is considered. In addition, the maximum positive value for force coefficient $C_y$ is obtained for tornado path with $l_y=3.0m$, while the maximum negative value is obtained for tornado path with $l_y=-3.0m$. By comparing the C_z peak values obtained with $l_y=-3.0m$ and $l_y=3.0m$ and the corresponding peak value obtained with $l_y=0$, it is observed that the vertical force induced by the tornado flow on the immersed object is smaller in the former cases.

A final comparison is performed considering an oblique tornado trajectory and the profile models proposed in this work, i.e. RCVM and L-O/B-R. It is observed that all the force components are significantly reduced when compared with the respective values obtained with the tornado paths investigated previously. In general, the maximum values predicted with the L-O/B-R model are greater than those obtained with RCVM, indicating that the L-O/B-R approach can maintain a tornado vortex with larger radius without excessive numerical dissipation.

Instantaneous velocity vector fields are now presented on plane $x_2=0$ m at $t=20s$ according to the different path trajectories analyzed here and using the L-O/B-R velocity profile model. Figure 18a shows a recirculation zone in front of the immersed object, where a region with high velocity values can be identified. At the upstream vertex of the cubic model one can also see a flow separation point, which leads to high vertical velocity values in its vicinity and a recirculation zone above the building top surface. Figure 18b presents a velocity field with relatively lower values than those observed in Fig. 18c, where the tornado trajectory leads to a velocity distribution in which the combination of translation and tangential velocities obtains its maximum values owing to the vortex rotation in the counterclockwise direction. This explains the high velocity field developed around the building, especially above the building position. In this case, intense recirculation can be identified on the lateral and superior walls of the immersed body.

Figure 18
Cubic building subject to three-dimensional tornado flow, velocity vector fields on plane $x_2=0$ at $t=20s$: (a) tornado path $l_y=0$ with L-O/B-R; (b) tornado path $l_y=-3.0m$ with L-O/B-R; (c) tornado path $l_y=3.0m$ with L-O/B-R.

Flow turbulence characteristics are evaluated here using time histories of flow variables obtained at four different points in the wake region of the cubic building, as indicated in Fig. 12. Results are presented in Table 12 for turbulence intensity and Reynolds stress components, which are calculated using the expressions given by Eq. (29), but without normalization, and tornado path $l_y=0$ with L-O/B-R, velocity profile. In this case, a centered moving time-average is adopted to better represent the flow turbulence variables, considering that the influence of the moving vortex on the flow characteristics is significant within a specific time interval. The moving time-average is calculated in this work using the following expression:

Table 12
Turbulence characteristics for cubic building submitted to moving vortex flow.

	Point 1	Point 2	Point 3	Point 4
Turbulence intensity - ${\overline{I}}_{v1}$	0.287	0.707	0.753	0.744
Turbulence intensity - ${\overline{I}}_{v2}$	0.232	1.113	1.160	1.038
Turbulence intensity - ${\overline{I}}_{v3}$	0.010	0.675	0.860	0.792
Reynolds stress component - ${\overline{{\tau }^{\prime }}}_{11}$	0.0825	0.499	0.567	0.553
Reynolds stress component - ${\overline{{\tau }^{\prime }}}_{12}$	-0.0131	-0.170	-0.321	0.00429
Reynolds stress component - ${\overline{{\tau }^{\prime }}}_{13}$	-0.000538	-0.0170	-0.106	-0.138
Reynolds stress component - ${\overline{{\tau }^{\prime }}}_{22}$	0.0540	1.238	1.345	1.077
Reynolds stress component - ${\overline{{\tau }^{\prime }}}_{23}$	0.000512	0.0409	0.0989	-0.0328
Reynolds stress component - ${\overline{{\tau }^{\prime }}}_{33}$	0.000107	0.456	0.739	0.628

$${\overline{q}}_m^n(t_i)=\frac{{\displaystyle \sum _{j=i-\frac{n-1}{2}}^{j=i+\frac{n-1}{2}}q(t_j)}}{n}$$ (30)

where ${\overline{q}}_m^n(t_i)$ is the moving time-average associated with a flow variable $q$ and evaluated at a time instant $t_i$, which is defined considering $n$ data samples. Notice that the number of data samples $n$ is an odd number and represents a discrete time interval where the time-average is calculated, such that the time instant $t_i$ is always located at the center of the discrete time interval. The discrete time interval $n$ is obtained here taking into account the nearest integer calculated from $n={(f_v.\Delta t)}^{-1}$, where $f_v$ is the vortex shedding frequency and $\Delta t$ is the time increment (see Table 7). In the present analysis, the time interval $[\mathrm{15,}30s]$ is adopted for evaluation of flow turbulence characteristics, as demonstrated in Fig. 19. Time-histories of the flow variables evaluated in the wake of the cubic model are presented in Fig. 20.

Figure 19
Centered moving average applied to the flow velocity component $v_1$ evaluated at point 4.

Figure 20
Time-histories of flow variables evaluated at different points in the wake of the cubic building.

5 CONCLUSIONS

In the present work an efficient finite element formulation was proposed to simulate vortex and tornado-like flows by using velocity profile models and linear hexahedral elements with one-point quadrature techniques. In addition, an ALE approach with automatic mesh motion scheme was utilized to describe interactions between vortex flows and immersed objects. Flow fields induced by tornado-like vortices were generated here considering time-dependent boundary conditions imposed on the computational domain according to the RCVM and L-O/B-R profile models, where LES was adopted for turbulent flows. Some applications with idealized flow conditions were initially simulated in order to verify the numerical formulation presented here and a final example was analyzed considering a cubic building subject to three-dimensional flow conditions.

In a first application, a two-dimensional circular cylinder was submitted to the flow generated by a stationary vortex, where the cylinder translates in the longitudinal direction of the computational domain. It was demonstrated that this approach is equivalent to consider a fixed cylinder submitted to a moving vortex, which rotates and translates simultaneously. The ALE formulation with the mesh motion scheme proposed here was able to simulate the flow characteristics expected in this case and the equivalence between the approaches was proved considering comparative predictions of force coefficients and vorticity fields. However, some differences were observed as the tangential velocity was increased indicating that mesh refinement plays a more important role when the tangential velocity is higher and the vortex core is small. In addition, mesh distortions due to the mesh motion scheme may lead to accuracy reductions in some mesh regions of the flow field.

In the second application, a cubic building subject to three-dimensional tornado-like vortex flow was analyzed, where four different tornado paths were considered. The aerodynamic forces induced by the flow on the building surface and flow conditions around the immersed object were obtained using the RCVM and L-O/B-R velocity profile models. It was observed that the vortex structures obtained with RCVM were generally smaller than the vortex structures predicted with L-O/B-R. The vortex diameter was continuously reduced during the time interval where the tornado vortex traveled through the computational domain and a significant change in the vortex structure was identified after it interacts with the immersed body. When the tornado vortex was passing the building zone, very high pressure suctions were developed on the building walls and a significant increase in the aerodynamic forces acting on the immersed body was observed according to the vortex rotation direction and rotation intensity. From the numerical investigations performed here, one can conclude that the flow rotational component has a significant influence on the magnitude of flow-induced forces for bodies immersed in a tornado-like flow, where the velocity ratio $V_{\theta ,\max }/U_x$ plays a major role. In a future work, the ALE approach proposed in this paper should be applied to three-dimensional tornado flow conditions in order to validate the present approach considering relative motions between immersed structures and tornado vortex.

Acknowledgements

The authors would like to thank the National Council for Scientific and Technological Development (CNPq, Brazil) and Brazilian Federal Agency for Support and Evaluation of Graduate Education (CAPES, Brazil) for the financial support. The present research was developed using computational resources provided by the National Center for Supercomputing (CESUP/UFRGS, Brazil), High Performance Computing Center (NACAD/UFRJ, Brazil) and National Center for High Performance Computing (CENAPAD/UNICAMP, Brazil).

• https://doi.org/10.1590/1679-78257419

References

• Alrasheedi, N.H.M. (2012). Computer Modeling of the Influence of Structure Plan Areas on Tornado Forces. University of Arkansas.
• Beaudan, P., Moin, P. (1994). Numerical Experiments on the Flow past Circular Cylinder at Subcritical Reynolds Number. Report N° TF-62, Department of Mech. Engr., Stanford University.
• Braun, A.L., Awruch, A.M. (2008). Finite element simulation of the wind action over bridge sectional models: application to the Guamá River Bridge. Finite Elements in Analysis and Design 44: 105-122.
• Braun, A.L., Awruch, A.M. (2009). A partitioned model for fluid-structure interaction problem using hexahedral finite elements with one-point quadrature. International Journal of Numerical Methods in Engineering 79: 505-549.
• Bruno, L., Khris, S. (2003). The validity of 2D numerical simulations of vertical structures around a bridge deck. Mathematical and Computer Modelling 37: 795-828.
• Burgers, J.M. (1948). A mathematical model illustrating the theory of turbulence. Advances in Applied Mechanics 1: 171-199.
• Cao, S., Wang, M., Cao, J. (2018). Numerical study of wind pressure on low-rise buildings induced by tornado-like flows. Journal of Wind Engineering and Industrial Aerodynamics 183: 214-222.
• Chien, N., Feng, Y., Wang, H-J., Siao, T-T. (1951). Wind-Tunnel Studies of Pressure Distribution on Elementary Building Forms. Iowa Institute of Hydraulic Research, State University of Iowa, Iowa City.
• Chorin, A.J. (1967). A Numerical Method for Solving Incompressible Viscous Flow Problems. Journal of Computational Physics 2(1): 12-26.
• Christon, M.A. (1997). A domain-decomposition message-passing approach to transient viscous incompressible flow using explicit time integration. Computer Methods in Applied Mechanics and Engineering 148(3-4): 329-352.
• Donea, J. (1984). A Taylor-Galerkin method for convective transport problems. International Journal for Numerical Methods in Engineering 20(1): 101-119.
• Guo, X., Cao, J. (2019). An IB-LBM investigation into the aerodynamic coefficients in relation to the rotation intensity of a tornado-like wind. Computers & Mathematics with Applications 78(4): 1206-1226.
• Hernandez, J.E.R., Sphaier, S.H. (1999). A Projection method combined with a finite volume method for unsteady Navier-Stokes equations compared with benchmark solutions. Hybrid Methods in Engineering 1(4): 1-26.
• Hughes, L.A. (1952). On the low-level structure of tropical storms. Journal of Meteorology 9(6): 422-428.
• Kareem, A., Kijewski, T. (1996). 7th US National Conference on Wind Engineering: A Summary of Papers. Journal of Wind Engineering and Industrial Aerodynamics 62(2-3): 81-129.
• Kawahara, M., Hirano, H. (1983). A finite element method for high Reynolds number viscous fluid flow using a two-step explicit scheme. Internationa Journal for Numerical Methods in Fluids 3: 137-163.
• Kim, Y.C., Matsui, M. (2017). Analytical and empirical models of tornado vortices: A comparative study. Journal of 1080 Wind Engineering and Industrial Aerodynamics 171: 230-247.
• Lewellen, W.S. (1976). Theoretical models of the tornado vortex. Proceedings, Symposium on Tornadoes: Assessment of Knowledge and Implications for Man, Texas Tech University, Lubbock, Texas.
• Rajani, B.N., Kandasamy, A., Majumdar, S. (2016). LES os flow past a circular cylinder at Re = 3900. Journal of Applied Fluid Mechanics 9(3): 1421-1435.
• Rankine, W.J.M. (1891). A Manual of Applied Mechanics, 13^th ed., London: Charles Griffin and Company.
• Rott, N. (1958). On the viscous core of a line vortex. Zeitschrift für angewandte Mathematik und Physik 9 (5-6): 543-553.
• Schlichting, H., Gersten, K. (2016). Boundary-Layer Theory, 9^th Edition, Springer-Verlag, Berlin.
• Selvam, R.P. (1985). Application of the boundary element method for tornado forces on buildings. Texas Tech University.
• Selvam, R.P. (1993). Computer modeling of tornado forces on buildings. Proceedings: The 7 th U.S. National Conference on Wind Engineering, Edited by: Gary C. Hart, Los Angeles, June 27-30, 605-613.
• Selvam, R.P., Millett, P.C. (2003). Computer modeling of tornado forces on a cubic building using large eddy simulation. Journal of Arkansas Academy of Science 57: 140-146.
• Selvam, R.P., Millett, P. (2005). Large eddy simulation of the tornado-structure interaction to determine structural loadings. Wind and Structures 8 (1): 49-60.
• Sengupta, A., Haan, F.L., P.P. Sarkar, Balaramudu, V. (2008). Transient loads on buildings in microburst and tornado winds. Journal of Wind Engineering and Industrial Aerodynamics 96(10-11): 2173-2187.
• Smagorinsky, J. (1963). General circulation experiments with the primitive equations. Monthly Weather Review 91(3): 99-164.
• Strasser, M.N. (2015). The aerodynamic and dynamic loading of a slender structure by an impacting tornado-like vortex: The influence of relative vortex-to-structure size on structural loading. University of Arkansas.
• Strasser, M.N., Selvam, R.P. (2015a). Selection of a realistic viscous vortex tangential velocity profile for computer simulation of vortex-structure interaction. Journal of the Arkansas Academy of Science 69: 88-97.
• Strasser, M.N., Selvam, R.P. (2015b). The variation in the maximum loading of a circular cylinder impacted by a 2D vortex with time of impact. Journal of Fluids and Structures 58: 66-78.
• Sullivan, R.D. (1959). A two-cell vortex solution of the Navier-Stokes equations. Journal of the Aerospace Sciences 26 (11): 767-768.
• Vatistas, G.H., Kozel, V., Mih, W.C. (1991). A simpler model for concentrated vortices. Experiments in Fluids 11(1): 73-76.
• Wanderley, J.B.V., Levi, C.A. (2002). Validation of a finite difference method for the simulation of vortex-induced vibrations on a circular cylinder. Ocean Engineering 29: 445-460.
• Wilson, T. (1977). Tornado structure interaction: A numerical simulation. Report: UCRL-52207, Distribution Category: UC-11, 80, Lawrence Livermore Laboratory, University of California, Livermore.
• Wood, V.T., White, L.W. (2011). A new parametric model of vortex tangential-wind profiles: development, testing, and verification. Journal of the Atmospheric Sciences 68(5): 990-1006.

Edited by

Editor: Pablo Andrés Muñoz Rojas

Publication Dates

• Publication in this collection
  29 May 2023
• Date of issue
  2023

History

• Received
  19 Dec 2022
• Reviewed
  10 Mar 2023
• Accepted
  24 Apr 2023