Application of Computer Graphics Flow Visualization Methods in Vortex Rope Investigations

Abstract

Computer graphics visualization techniques for application on data from Computational Fluid Dynamics (CFD) simulations of the vortex rope, a phenomenon present in hydraulic turbines operating in off-design conditions, were devised. This included not only objects for visualization (what to visualize) but also methods of the visualization itself (how to do it). By means of advanced methods based particularly on volume rendering of Eulerian fields in combination with Lagrangian objects, various phenomena were captured, such as the motion of the vortex rope or the backflow zone. The data came from simulations using a scale-resolving hybrid turbulence model, the Stress-Blended Eddy Simulation. In such detailed simulations and other applications involving complex three-dimensional structures, proper visualization methods are needed to leverage the content captured in the resultant data.

1. Introduction

Computer graphics flow visualization belongs to the field of data visualization. The main goal is to convert the data into visual representations, mostly 2D pictures viewed on a screen or printed on paper. Generally, the data are 3D, multivariate, and time-dependent. The process is, in such complex cases, typically irreversible, meaning that it is not possible to obtain all of the data back from the resultant pictures, and a portion of the information is inevitably lost. This gives rise to the central questions of data visualization: what to visualize and how in order to capture all the relevant features properly. Here, we give a brief overview of currently used computer graphics flow visualization techniques. For a more comprehensive description, we refer to Reference [1].

The process begins with the generation of the data. They may come from a numerical simulation (Computational Fluid Dynamics (CFD)) or an experiment. The second step comprises data modification. One can define and compute new variables from the existing ones, restrict the domain of interest by clipping, or interpolate the data to a different grid or data points. New objects can also be generated. These include streamlines and stream surfaces or tubes, vortex lines and vortex tubes, particles and their trajectories, streaklines, and material lines or surfaces. The third step is the generation of objects that will be rendered in the next step. Three-dimensional scalar fields can be displayed on a series of 2D slices. One can also extract isosurfaces. Streamlines and other line objects can be visualized as tubes. Various glyphs can be used for visualization of particles, vectors, or even tensors [2]. The data values on these objects can be coded by their color, transparency, material properties, textures, or dimensions. For an example with application on a vortex rope case, we refer to Reference [3], where quantities related to vorticity transport were visualized on pathlines using variable radius, segmentation, and different colors. The final step is rendering the scene and displaying the result. Scalar fields can be rendered directly by volume rendering. It is based on a fictitious light particle moving through the data field. Physically-based phenomena, such as emission, scattering, and absorption, can be considered [4]. Their parameters are based on the data values, i.e., we map the data values to a parameter, e.g., the emission strength, using a defined transfer function. In the end, we obtain a 2D field to be displayed on the screen. The development of this technique is driven particularly by medicine with applications to computed tomography (CT) or magnetic resonance imaging (MRI) data [5]. The art of volume rendering resides in the selection of the emission, scattering, and absorption parameters, where it is advisable to make unimportant parts of the data domain completely transparent and find such settings for the rest that provide the best notion of how the scalar values are distributed in the 3D space. For surface rendering, there exist a variety of methods ranging from simple targeted at speed to photorealistic developed primarily for computer-generated imagery (CGI) with applications including movies and video games.

Swirling flows are typical for many engineering applications ranging from aircraft, process engineering, turbomachinery to combustion. These flows are rather susceptible to instabilities, one of them with important technical implications is vortex breakdown, an abrupt change in the flow structure of the swirling flow. Generally accepted definition of vortex breakdown comes from Leibovich [6], who defined it as a disturbance characterized by formation of an internal stagnation point on the vortex axis followed by reversed flow in a region of limited extent. In addition, the definition of Lugt [7] is often cited, who stated that vortex breakdown is a closed region or spiral within swirling flow with streamline divergence near the axis. It was also Lugt, among others, who noted that the nature of the spiral vortex breakdown changes with inlet swirl and that pressure increases along the vortex axis, which can be further risen by streamline divergence in a diffuser. The pressure increase is accompanied by an axial velocity drop and the formation of a stagnant region. Since its first documented observation in 1957 [8], an enormous effort was invested in study of this phenomenon and lots of papers were published (see Reference [9,10] or [11] for extensive reviews). It is not surprising that most of them contain some form of visual representation of the spiral vortex breakdown. Visualization often plays a significant role in qualitative understanding and improves insight into complex flow behavior.

The vortex rope, a case of spiral vortex breakdown, is an unwanted phenomenon encountered in the draft tubes of water turbines [12]. It is a coherent structure, i.e., it persists its characteristics for longer times than turbulent eddy timescales. Due to its coherence, it may cause significant vibrations and noise. To overcome these effects and provide a wider operating range, efficient flow control methods aimed at its disintegration need to be employed [13,14,15,16]. Knowledge of its behavior is a prerequisite, leading to studies of this phenomenon by means of CFD [17,18,19] or stability analysis [20]. Such calculations produce data that need to be visualized.

Visualizations in published CFD studies of the vortex rope typically rely on simple techniques. The vortex rope itself is visualized usually as an isosurface of a suitable quantity, most often the static pressure. Sometimes, a vortex identification method is used, e.g., the Q-criterion by Hunt [21], the ${\lambda }_2$-criterion by Jeong and Hussain [22], or the swirling strength by Zhou et al. [23], to name the most popular. For a review of their application in hydraulic machinery, see Reference [24]. Rajan and Cimbala [25] compared the static pressure and the mentioned vortex identification methods on a case with the vortex rope corresponding to the FLow INvestigation in a Francis Draft Tube (FLINDT) project [26]. They used a suitably selected isosurface and found out that visualizations based on the static pressure capture only the upper part of the vortex rope, while the vortex criteria capture the whole structure and its disintegration into smaller vortices downstream.

In this paper, we extend the previous works and use more elaborated vortex identification methods. These include the ${\mathit{M}}_{\mathrm{Z}}$-criterion by Haller [27] and the residual vorticity by Kolář [28]. Another aim of this paper is to devise computer graphics flow visualization methods suitable for studies of the vortex rope and possibly other similar phenomena. We also propose a way of visualizing the whole 3D field of the finite-time Lyapunov exponent, a parameter connected with so-called Lagrangian coherent structures [29]. The methods are demonstrated on data from CFD simulations of flow downstream of a swirl generator that cover a wide range of operating points given by the swirl number.

2. Vortex Rope Visualizations

2.1. Test Case

The vortex rope phenomenon is usually studied in model turbines or swirl generators. The former mentioned provide conditions similar to those in actual water turbines, but they are usually much more expensive compared to simple swirl generators. We have designed such a generator without any moving parts, capable of reaching a wide range of operating conditions in terms of Reynolds and swirl numbers. This is achieved by two separate inlets, axial and tangential, leading into a volute with a central hub, and the operating point is given by the component flow rates. The geometry is sketched in Figure 1. An advantage is that the operating point can be changed smoothly during operation, but there is also a disadvantage of vortices generated in the mixing layer between axial and tangential flows that can affect the studied vortex rope, especially in cases with lower swirl numbers.

2.2. CFD Simulation

Geometry of the generator and the diffuser was filled with 7 million hexahedral cells. While several non-confomal interfaces are between mesh blocks within the generator, the diffuser mesh is completely conformal with $y^{+}<5$. This parameter is given by

$$y^{+}=\frac{\rho v_{\tau }y}{\mu },v_{\tau }=\sqrt{\frac{{\tau }_{\mathrm{w}}}{\rho }},{\tau }_{\mathrm{w}}={\left.\mu \frac{\partial v_{\mathrm{t}}}{\partial y}\right|}_{y=0},$$

(1)

where y is the wall normal distance, $\rho$ the density, $\mu$ the dynamic viscosity, and $v_{\mathrm{t}}$ the velocity component parallel to the wall. A slice through this mesh is presented in Figure 2, which illustrates the maximum effort spent to achieve orthogonality and equilateral edges of computational cells. Velocity boundary conditions were prescribed on the axial and tangential inlets and a constant static pressure on the outlet from the computational domain. Modeling of the turbulent flow was carried out with Stress-Blended Eddy Simulation (SBES) turbulence model [30]. SBES is a hybrid URANS-LES approach, based on an improved asymptotic shielding of RANS boundary layer against LES modification. Implementation of SBES within ANSYS Fluent 18.1 is employed, where the RANS part (mostly near-wall regions) is solved using SST k-$\omega$ model and the LES part (core of the flow containing a vortex rope) with wall-adapting local eddy-viscosity model (WALE LES). Its formulation of eddy viscosity is based on the following concept:

$${\nu }_{\mathrm{t}}^{\mathrm{SBES}}=f_{\mathrm{s}}{\nu }_{\mathrm{t}}^{\mathrm{RANS}}+(1-f_{\mathrm{s}}){\nu }_{\mathrm{t}}^{\mathrm{LES}},$$

(2)

where $f_{\mathrm{s}}$ is the shielding function, and ${\nu }_{\mathrm{t}}$ is the eddy viscosity. As can be seen, the shielding function controls the blending of RANS and LES closures. It was designed with respect to a list of requirements, some of them being preventing the boundary layers from being solved by LES (otherwise, much finer meshes are needed) and quick switching to LES after their detachment. For more information, we refer to Reference [30]. We used the following solver settings: PISO algorithm (Pressure-Implicit with Splitting of Operators) for pressure–velocity coupling, least squares cell-based scheme for gradients discretization, PRESTO! (PREssure STaggering Option) for pressure discretization, bounded central differencing for momentum discretization, second order upwind scheme for turbulent kinetic energy and turbulent dissipation discretization, and, finally, bounded second order implicit scheme for temporal discretization.

As the simulations cover different operating points, we introduce the swirl number, a parameter quantifying the intensity of swirl in circular cross sections. It is an integral parameter with different definitions used in the literature. The following definition [31] of $\mathrm{Sr}$ is used throughout the paper:

$$\mathrm{Sr}=\frac{1}{R}\frac{{\int }_{S}r\left|v_{\phi }{v}_{\mathrm{z}}\right|\mathrm{d}S}{{\int }_S{v_{\mathrm{z}}}^2\mathrm{d}S},$$

(3)

where S is a selected circular cross section, R its wall radius, $v_{\phi }$ the tangential velocity, and $v_{\mathrm{z}}$ the axial velocity. A plane right under the hub tip was selected for the evaluation of $\mathrm{Sr}$.

Visualization of swirling flow within the diffuser was carried out for 2 snapshots, which correspond to instantaneous flow fields characterized by $\mathrm{Sr}\approx 0.45$ and $\mathrm{Sr}\approx 0.85$. The total flow rate (and hence also the Reynolds number) was held constant (5 L/s).

Computational data presented in this paper are part of extensive research aimed on investigation of vortex rope behavior and its control, where visualization is one of partial goals. The reader is referred to Reference [32,33] for details about experiments and CFD simulations.

2.3. Software and Methodology

Visualizations presented in this paper were created in ParaView 5.8.0 and Blender 2.82a (both are open-source tools). ParaView is a software for data analysis and visualization built on VTK (Visualization Toolkit) [34] and was employed for the pre-rendering operations (generating 3D surface objects, such as tubes for streamlines, isosurfaces, etc., and generating series of images for volume rendering). Finite-time Lyapunov exponents, residual vorticity, and ${\mathit{M}}_{\mathrm{Z}}$-criterion were computed using our own Python programmable filters. Blender is a 3D creation suite offering tools for the whole pipeline from modeling to rendering, animation and video editing. Its Cycles renderer offers more control over lighting and material effects and well combines multiple volume renderings. We used it for rendering the objects generated in ParaView. The 3D surface objects were exported from ParaView in X3D format. Importing them to Blender and setting up for rendering is straightforward, with many guides publicly available. On the other hand, volume rendering of data fields in Blender is more complicated. A guide through this procedure is, therefore, given in Appendix A.

2.4. Vortex Identification Methods

The central task of vortex rope visualizations is how to identify the vortex rope itself. The simplest way is to exploit the fact that static pressure decreases towards the vortex rope core. This approach is especially suitable for quick explorations, as the static pressure is readily available from CFD.

Another way is to use vortex identification methods. Although the intuitive notion of a vortex as a portion of fluid rotating around a common axis (vortex core) is clear, there is no widely accepted rigorous definition capable of clearly distinguishing the vortices among other possible flow structures. Kolář [28] presented a review of the most popular vortex identification methods and summarized the requirements on vortex definitions: objectivity (independence of the reference frame motion), validity for compressible and variable-density flows, and ability to determine the local intensity of the swirling motion and the swirl orientation, among others. The problem resides in the inability of the methods to comply with all the requirements and in the arguability of some of them.

Vortices are naturally connected with the vorticity, a quantity defined as the curl of the flow velocity

$$\mathsf{\Omega }=\nabla \times \mathit{v}.$$

(4)

It is a measure of the local spinning motion, precisely twice its angular velocity. This motion, however, can also be caused by purely shearing flows as is the case of the boundary layers. The Q-criterion overcomes this problem by defining vortices as the regions where the antisymmetric part $\mathit{R}$ of the velocity gradient tensor prevails over its symmetric part $\mathit{S}$ in the sense of the Frobenius norm, i.e., $\parallel \mathit{A}\parallel =\sqrt{{\sum }_{i,j}{A_{ij}}^2}$. We note that the nonzero elements of the antisymmetric part correspond to the components of the vorticity vector. The Q-criterion is defined as follows [21]:

$$Q=\frac{1}{2}\left({\parallel \mathit{R}\parallel }^2-{\parallel \mathit{S}\parallel }^2\right).$$

(5)

This is also the second invariant of the velocity gradient tensor. This tensor and its properties gave rise to other vortex criteria.

It can be shown that if the streamlines form closed or spiraling shapes, the characteristic equation of the velocity gradient tensor has one real and two complex solutions, leading to a pair of complex conjugate eigenvalues:

$${\mathsf{\Lambda }}_1\in \mathbb{R},{\mathsf{\Lambda }}_2=a+b\mathrm{i},{\mathsf{\Lambda }}_3=a-b\mathrm{i}.$$

(6)

The swirling strength is then defined as the magnitude of their imaginary part, i.e., b. The idea is that the vortices are the regions of b > 0, with b serving also as a measure of their strength [23].

Another popular vortex criterion, the ${\lambda }_2$-criterion, is based on different considerations. It is known that swirling motion itself leads to a local pressure minimum at the vortex axis. We can take the gradient of the Navier–Stokes equation assuming incompressible fluid and no external forces:

$$\nabla \left(\frac{\partial \mathit{v}}{\partial t}+\left(\mathit{v}·\nabla \right)\mathit{v}-\nu \Delta \mathit{v}\right)=-\frac{1}{\rho }\nabla \left(\nabla p\right).$$

(7)

If the Hessian matrix $\nabla \left(\nabla p\right)$ has two positive eigenvalues, then there exists a plane in which p reaches its local minimum in the point of interest. As the next step, the velocity gradient is decomposed into the antisymmetric part $\mathit{R}$ and the symmetric part $\mathit{S}$. Some of the terms on the left-hand side vanish, as they correspond to the vorticity transport equation. We then arrive at the following form:

$$\frac{\mathrm{D}\mathit{S}}{\mathrm{D}t}+{\mathit{S}}^2+{\mathit{R}}^2-\nu \Delta \mathit{S}=-\frac{1}{\rho }\nabla \left(\nabla p\right).$$

(8)

It is well-known that the local pressure minimum does not necessarily guarantee the existence of a vortex, and vice versa [28]. For this reason, the contribution of the material derivative of $\mathit{S}$ and the viscous term to the Hessian matrix is discarded, and the ${\lambda }_2$-criterion requires two eigenvalues of ${\mathit{S}}^2+{\mathit{R}}^2$ to be negative. The eigenvalues are usually sorted in ascending order, and a vortex is defined as the region where ${\lambda }_2$ is negative [22]. The absolute value of this eigenvalue can serve as a measure of the local swirl intensity.

A more elaborate vortex criterion was proposed by Kolář [28]. It tackles the already mentioned drawback of the vorticity, i.e., admitting nonzero values not only for purely rotating flows but also for purely shearing flows. The proposed method is based on a triple decomposition of the velocity gradient into so-called effective shear, residual strain and residual rotation parts. The last part leads to the residual vorticity, i.e., a portion of the original vorticity related only to shearless motion. A bottleneck of this method is the effective shear extraction procedure. First, the decomposition of the velocity gradient tensor is defined as follows:

$$\nabla \mathit{v}=\nabla {\mathit{v}}_{\mathrm{shear}}+\nabla {\mathit{v}}_{\mathrm{residual}},$$

(9)

$$\nabla {\mathit{v}}_{\mathrm{residual}}=\mathsf{{\rm Y}},{\mathsf{{\rm Y}}}_{ij}=\left(\mathrm{sgn}\frac{\partial v_{\mathrm{i}}}{\partial x_j}\right)min\left(\left|\frac{\partial v_i}{\partial x_j}\right|,\left|\frac{\partial v_j}{\partial x_i}\right|\right).$$

(10)

The shear tensor is by definition purely asymmetric, and this procedure is proposed so that this property is satisfied. The results, however, depend on the choice of the frame of reference. We seek the frame where the Frobenius norm of the residual tensor is minimized. The shear part in this so-called best reference frame (BRF) is denoted as the effective shear. Kolář also derived a relation between $\parallel \nabla \mathit{v}\parallel$ and $\parallel \nabla {\mathit{v}}_{\mathrm{residual}}\parallel$:

$$\parallel \nabla {\mathit{v}}_{\mathrm{residual}}{\parallel }^2+4\left(|S_{12}{R}_{12}|+|S_{23}{R}_{23}|+|S_{31}{R}_{31}|\right)={\parallel \nabla \mathit{v}\parallel }^2.$$

(11)

Considering the fact that $\parallel \nabla \mathit{v}\parallel$ remains unchanged in all frames, the BRF is defined as the frame where $|S_{12}{R}_{12}|+|S_{23}{R}_{23}|+|S_{31}{R}_{31}|$ is maximized. The computation of the residual vorticity therefore involves time-consuming optimization. Once the BRF is found, the residual tensor is computed, and then the residual vorticity is extracted from its antisymmetric part and returned back to the original frame.

An advantage of vorticity-based criteria is that the vorticity determines also the swirl orientation. The axial vorticity can be used to determine the sense of rotation of the vortex rope (see the following section). The azimuthal vorticity and its gradient can be used to detect the onset of vortex breakdown [35].

We have visualized the vortex rope by volume rendering of these quantities in Figure 3. The colorbar range was determined with the target of reaching maximal similarity of the representations. The opacity transfer function represented by the white line is of the same shape for all cases to ensure comparability of the scaling. Two snapshots with different swirl numbers are depicted to assess also the dependence of each variable on the operating point.

The static pressure turns out to be highly sensitive to the swirl number (decreasing rapidly with increasing swirl number). For the case with the higher swirl number, it gives qualitatively different results compared to the other quantities. It highlights a whole spiraling structure rather than individual vortices. Together, the results show that the flow in the diffuser in this case exhibits a huge spiraling structure (highlighted by the static pressure) which consists of twin spiraling vortices under the hub and many smaller vortices downstream (highlighted by the other quantities). The values of the Q-criterion and ${\lambda }_2$-criterion scale rapidly, resulting in less-pronounced weaker vortices compared to the swirling strength or vorticity fields. The triple decomposition attributes around 2/3 of the total vorticity in the vortex rope to shearing motion, leading to an interesting result where even smaller vortices are clearly distinguishable. Overall, all the vortex identification methods highlighted structures of practically the same shape but with a different scaling of the values quantifying their strength. Definition of the right threshold remains a general problem for most of the vortex criteria.

The vortex identification methods presented so far all have a common drawback: they are only Galilean invariant, i.e., their results do not depend on translation of the reference frame but do depend on its rotation. There is usually not a single natural choice of a reference frame, which leads to the need of an objective (independent of the frame motion) definition of a vortex [27].

Haller [36] showed that the deviation of the vorticity from its spatial mean is objective. In our case, the spatial mean in the diffuser turns out to be negligible and so can be omitted. The same principle can be applied to the residual vorticity. Haller [27] also proposed another objective vortex definition which is based on Lagrangian considerations. He partitioned the spatial domain into hyperbolic, parabolic, and elliptic zones based on the behavior of particle trajectories. The complete theory is beyond the scope of this work, only the computation procedure will be presented. We note that this criterion is applicable on incompressible fluid only. The strain acceleration tensor is given by

$$\mathit{M}=\frac{\mathrm{D}\mathit{S}}{\mathrm{D}t}+{\left(\nabla \mathit{v}\right)}^{\mathrm{T}}\mathit{S}+\mathit{S}\nabla \mathit{v}.$$

(12)

The elliptic zone is the region where $\mathit{M}$ is indefinite with respect to vectors belonging to certain defined subspace. To determine this, we need to compute the eigenvalues s of the strain rate tensor $\mathit{S}$. Two of them have the same sign, and the last one has the opposite sign. They are ordered as follows:

$$\mathrm{sgn}\left(s_1\right)=\mathrm{sgn}\left(s_2\right)\ne \mathrm{sgn}\left(s_3\right),|s_1|\ge |s_2|.$$

(13)

With the substitutions

$$a=-\frac{s_1}{s_3},\widehat{\mathit{M}}={\mathit{E}}^{\mathrm{T}}\mathit{M}\mathit{E},$$

(14)

where $\mathit{E}$ is the matrix of eigenvectors of the strain rate tensor, function $m\left(\alpha \right)$ is defined as follows:

$$\begin{array}{cc}\hfill m\left(\alpha \right)& =(1-a){\widehat{M}}_{11}{cos}^2\alpha +a{\widehat{M}}_{22}{sin}^2\alpha +a(1-a){\widehat{M}}_{33}\hfill \\ \hfill & +\sqrt{a(1-a)}\left(2\sqrt{1-a}{\widehat{M}}_{13}cos\alpha +2\sqrt{a}{\widehat{M}}_{23}sin\alpha +{\widehat{M}}_{12}sin2\alpha \right).\hfill \end{array}$$

(15)

The elliptic zone is the region where $\min \left(m\left(\alpha \right)\right)<0$. Testing this condition requires a suitable numerical method. A vortex is then defined as a set of trajectories remaining in the elliptic zone for the whole integration time. Unfortunately, as Haller pointed out, this criterion requires many numerical operations including second order spatial differentiation and computation of particle trajectories, leading to numerical errors that consequently incorrectly disqualify a significant portion of particles from being contained in a vortex. For this reason, a relaxed definition is used: a trajectory is contained in a vortex, if the ratio of the time spent in the elliptic zone to the integration time is above a selected threshold. This criterion is referred to as the ${\mathit{M}}_{\mathrm{Z}}$-criterion.

2.5. Hybrid Methods for Investigations of the Vortex Rope

2.5.1. The ${\mathit{M}}_{\mathrm{Z}}$-Criterion Combined with the Residual Vorticity

The ${\mathit{M}}_{\mathrm{Z}}$-criterion by definition gives no information about the intensity of the swirling motion. For this reason, we suggest combining with the residual vorticity. From the computational domain, we select only those data points (particles) that satisfy the ${\mathit{M}}_{\mathrm{Z}}$-criterion. These points are then colored according to the residual vorticity and rendered by volume rendering (Figure 4). At first, one needs to choose appropriate integration time. For each case, we present a longitudinal slice of an area under the hub tip with the vortex rope, colored by ${\mathit{M}}_{\mathrm{Z}}$, for three different integration times T. Generally, as T grows, more particles get into the hyperbolic zone, leading to the decay of their ${\mathit{M}}_{\mathrm{Z}}$ value. For the case with the lower swirl number, we observe that particles starting from the vortex rope core remain in the elliptic zone for long times, while particles starting from its surroundings spend considerable portion of time in the hyperbolic zone. There is a sharp gradient of ${\mathit{M}}_{\mathrm{Z}}$ values that clearly highlights the vortex core edge. For the case with the higher swirl number, the situation is more complicated, as there are many smaller vortices concentrated in the spiraling structure. Together with shorter time scales, this increases the uncertainty of the results. The volume renderings with the threshold of ${\mathit{M}}_{\mathrm{Z}}>0.99$ contain the vortex rope in the first case, but miss many structures with a high residual vorticity, especially the twin spiraling vortices under the hub tip, in the second case. A lower threshold of ${\mathit{M}}_{\mathrm{Z}}>0.9$ is fulfilled by substantially more particles in the second case, leading to much more structures visible in the volume rendering. In the first case, the change is less significant.

2.5.2. Combination of the Static Pressure and the Residual Vorticity

Figure 3 shows that, for the higher swirl number, the static pressure highlights the whole region of spiraling motion, whereas other criteria highlight small vortices. This led us to the idea of combining them using volume rendering, similarly as isosurfaces of a given quantity are often colored by another quantity. Two possibilities are depicted in Figure 5. If the transparency is set according to the static pressure and colors according to the residual vorticity, we see the spiraling region mostly in black, indicating low values of the residual vorticity. Of the vortices with a high residual vorticity, only those in the front are visible. The advantage resides in the fact that the vortices are better distinguishable in this representation than in the volume rendering of the static pressure field in Figure 3. If the transparency is set according to the residual vorticity and colors according to the static pressure, we get information about the pressure in smaller vortices. Generally, combining two quantities in one representation this way can elucidate their mutual relation.

2.5.3. Combination of the Axial Vorticity and Vortex Identification Methods

The axial vorticity is the parameter used to distinguish vortices rotating clockwise from those rotating counterclockwise with respect to the diffuser axis. The vortex criterion of interest must be adjusted to return positive values only. This adjusted criterion is then multiplied by the sign of the axial vorticity. Here, we define a modified residual vorticity as follows:

$${\widehat{\mathsf{\Omega }}}_{\mathrm{RES}}=\parallel {\mathsf{\Omega }}_{\mathrm{RES}}\parallel \mathrm{sgn}\left({\mathsf{\Omega }}_{{\mathrm{RES}}_{\mathrm{z}}}\right).$$

(16)

Visualizations are presented in Figure 5. A different technique was employed: the structures were rendered as light-emitting, no external lighting was applied. It better highlights the structures thanks to the contrast between them and the dark diffuser wall in the background. Negative axial vorticity prevails, and the vortex rope is rotating counterclockwise. Clockwise rotating vortices are present, as well. They play an important role when the swirl number is low. There is a bundle of vortices rotating in both directions under the hub tip. In the situation in Figure 5 (bottom right), a clockwise rotating vortex is breaking a counterclockwise rotating filament that previously gained some intensity and became the most significant vortex in the bundle.

2.6. Finite-Time Lyapunov Exponent

The finite-time Lyapunov exponent (FTLE) is an objective (independent of reference frame) Lagrangian quantity connected with Lagrangian coherent structures (see Reference [29] for a comprehensive review). Its applications comprise, e.g., studies of general turbulence [37], flow over an airfoil [38], motion of animals [39,40], and oceanic eddies [41]. The FTLE is a measure of separation of neighboring particles in time, i.e., the local maximizing curves (ridges) of the FTLE field may indicate either locally maximal stretching or locally maximal shear [37]. When integrating forward in time, high values of FTLE signalize repelling structures. Similarly, backward integration leads to attracting structures. Details on the procedure can be found in Reference [42]. Consider a function that for each initial position ${\mathit{x}}_0$ of a particle, initial time $t_0$ and final time t returns the position $\mathit{x}$ of the particle in the final time t. This so-called flow map can be written as

$${\phi }_{t_0}^t\left({\mathit{x}}_0\right)=\mathit{x}(t_0,{\mathit{x}}_0,t).$$

(17)

Now, take two initial positions: ${\mathit{x}}_0$ and ${\mathit{x}}_0+\mathrm{d}\mathit{x}$. Flow map ${\phi }_{t_0}^t({\mathit{x}}_0+\mathrm{d}\mathit{x})$ can be expressed using the Taylor series. Assuming infinitesimally small $\mathrm{d}\mathit{x}$, we obtain

$${\phi }_{t_0}^t({\mathit{x}}_0+\mathrm{d}\mathit{x})-{\phi }_{t_0}^t\left({\mathit{x}}_0\right)=\Delta \mathit{x}=\nabla {\phi }_{t_0}^t\left({\mathit{x}}_0\right)\mathrm{d}\mathit{x}.$$

(18)

This equation represents the distance $\Delta \mathit{x}$ in time t between particles seeded in points ${\mathit{x}}_0$ and ${\mathit{x}}_0+\mathrm{d}\mathit{x}$ in time $t_0$. Such $\mathrm{d}\mathit{x}$ is searched that maximizes the Euclidean norm of $\Delta \mathit{x}$:

$$∥\Delta \mathit{x}∥=\sqrt{\mathrm{d}\mathit{x}·\left[{\nabla {\phi }_{t_0}^t\left({\mathit{x}}_0\right)}^{\mathrm{T}}\nabla {\phi }_{t_0}^t\left({\mathit{x}}_0\right)\mathrm{d}\mathit{x}\right]}=\sqrt{\mathrm{d}\mathit{x}·\left(\mathit{D}\mathrm{d}\mathit{x}\right)},$$

(19)

where $\mathit{D}$ is denoted as the Cauchy-Green tensor. Assume the standard eigenvalue problem $\mathit{De}=\lambda \mathit{e}$, where $\mathit{e}$ is an eigenvector and $\lambda$ an eigenvalue. Then, $∥\Delta \mathit{x}∥$ is maximized, if the direction of $\mathrm{d}\mathit{x}$ is the direction of the eigenvector ${\mathit{e}}_{\max }$ of $\mathit{D}$ corresponding to the largest eigenvalue ${\lambda }_{\max }=max\left(\lambda \right)$. One can then rewrite this equation using $\mathrm{d}{\mathit{x}}_{\max }={\mathit{e}}_{\max }∥\mathrm{d}{\mathit{x}}_{\max }∥$ and $\mathit{D}\mathrm{d}{\mathit{x}}_{\max }={\lambda }_{\max }{\mathit{e}}_{\max }∥\mathrm{d}{\mathit{x}}_{\max }∥$ as follows:

$$∥\Delta {\mathit{x}}_{\max }∥=\sqrt{{\lambda }_{\max }}∥\mathrm{d}{\mathit{x}}_{\max }∥=∥\mathrm{d}{\mathit{x}}_{\max }∥exp\sigma |t-t_0|,$$

(20)

where the finite-time Lyapunov exponent is given by

$$\sigma =\frac{1}{|t-t_0|}ln\sqrt{{\lambda }_{\max }}.$$

(21)

The flow map can be considered in both positive and negative time, justifying the absolute values of the time difference in the previous equations.

To its advantage, the FTLE was found to be relatively insensitive to imperfect velocity data [43]. Its principal disadvantage, however, is the cost associated with the computation of particle trajectories for each seed point in the domain in order to compute the flow map. The distribution of the seed points needs to be fine enough to reach the desired level of details, leading to a high number of trajectory computations for each time instant. To this point, efficient methods for computing a sequence of FTLE fields in time were proposed [44]. The Eulerian representation of the FTLE is then obtained by spatial interpolation.

FTLE fields are usually rich in fine structures, making the visualization challenging in 3D. We again rely on volume rendering. The results also depend on the chosen integration time. As the time increases, the structures become sharper and clearer [37]. In Figure 6, we present visualizations of the FTLE fields in a longitudinal slice and in the whole diffuser. In the longitudinal slice, we observe significant FTLE ridges surrounding the vortex core, and a local minimum of FTLE in its center. This indicates solid-body-like rotation. In the pictures of the whole diffuser, an isosurface of the Q-criterion is added to complement the volume renderings of FTLE fields. It fits well into the visible FTLE structures and enables us to see their relation to the vortex rope core. The opacity transfer function makes the regions of lower FTLE values completely transparent and sharply increases at the end of the range to better highlight the most significant structures. The combination of the forward and backward FTLE fields was made within a single object using Mix Shader tool in Blender software.

2.7. Visualization of the Features of the Vortex Rope

The features of flows with a developed vortex rope to be captured by visualizations are:

• the vortex rope itself,
• backflow zone,
• swirling flow at the outer wall, and
• unsteadiness and regularity of the flow.

Here, methods are devised to tackle these challenges. Figure 7 shows visualizations for different operating points addressing the first three features. The vortex rope is identified by the modified residual vorticity. Together with the backflow zone, it was visualized by volume rendering, which makes it possible to distinguish the values of these variables in the selected visibility range, unlike when using only isosurfaces. The scene is complemented with portions of streamlines (grayscale) seeded throughout the whole domain, giving a notion of the velocity field. On two short lines along the wall at the inlet of the diffuser, we seeded a series of streamlines (green) that inform us about the swirl intensity. Overall, this scene shows in a single series of pictures, how different aspects of the diffuser flow change with the operating point. The fact that they are combined in a single picture reveals their mutual relations, e.g., how the backflow zone is related to the vortex rope and its motion, or how the vortex rope shape is related to the shape of the spiraling streamlines in the outer zone of the diffuser.

As for the last point, we present two scenes in Figure 8. The first one depicts a case with a high swirl number, where a dominant frequency in the discrete Fourier transform (DFT) spectrum was detected. The spectrum is included in the figure. It is in fact the average of frequency spectra computed in 140 points distributed throughout the diffuser. The dominant frequency is of the precessing motion of the spiraling structure and allows us to extract eight snapshots in the same phase and render them simultaneously. We used volume rendering of the static pressure and included also center lines to better highlight the structure shape variation. We note that the definitions aimed to identify the vortex core lines, such as the Levy criterion [45], Sujudi and Haimes algorithm [46], and their modifications [2], would in this case produce many lines belonging to the smaller vortices contained in the whole spiraling structure. Nonetheless, in cases with a single pronounced spiraling vortex rope, they can be utilized, as well. It is analogous to the behavior of vortex criteria in Figure 3. This is why we, in this case, resorted to a method based on the static pressure field, as this variable was shown to clearly highlight the whole spiraling structure. We shall generally have in mind, as was already mentioned, that the local pressure minimum does not necessarily guarantee the existence of a vortex. The algorithm that we proposed for computing the center lines works as follows:

• Extract the cells with the pressure value below certain threshold. The threshold should be selected so that the vortex rope is present, while any other structures are excluded. This limits this approach to states with a distinct vortex rope only. In other words, we wish to extract the vortex rope from the whole domain.
• Perform a plane slice passing through the top part of the vortex rope. The plane normal and the vortex rope axis should optimally be collinear.
• Compute a suitable average (e.g., pressure-weighted) of the cell centers. This will be the first point of the resultant center line.
• Compute the center of the second plane slice. It lies in the direction of the first plane normal at a given distance from the first point of the center line. The second plane normal is the same as the first plane normal. Perform the slice.
• Compute the average of the cell centers. This will be the second point of the resultant center line.
• From now on, the algorithm runs automatically. The normal of the next plane slice is given by the direction vector joining the last two points of the center line. Its center lies at a given distance from the last point in the direction of the plane normal. The next point of the center line is again defined as an average of the cell centers.

This algorithm works as long as the slice passes through the vortex rope only once. If the vortex rope is so twisted that this is not the case, the correct cluster of cells needs to be determined to obtain proper results. Otherwise, the resultant point would lie halfway between the actual core points. A similar situation actually happens in the upper part, where two intertwined vortices (sometimes even three or four were observed) are present, and, in the slices, there are two significant pressure minima that correspond to their cores. The proposed algorithm produces a single line passing halfway between them, which is in this case the desired outcome, as we are interested in the shape of the whole spiraling structure rather than individual vortices. Enhancements treating possible undesired cases, however, are out of the scope of this work.

The results give a notion of the vortex rope regularity. In this case, significant differences in its pitch are observed. Portions of streamlines initialized from the same point inform about the regularity of the velocity field. We can observe significant irregularity that is lesser at the walls and greater at the axis.

The second picture shows the motion of the vortex rope. It is typically visualized by a series of pictures, but, in some cases, it can be beneficial to combine more snapshots in one picture. Here, we have a single pronounced vortex rope visualized by the Q-criterion. Vortex core lines (black in Figure 8) based on the definition by Levy [45] are also included. We computed them using a ParaView plugin by VCG Heidelberg [47]. As the time increases, the vortex rope fades out, which helps with orientation in time. It is clearly visible that some parts of the vortex rope are moving with a lower axial velocity than others. Such inconsistencies may eventually lead to the formation of a vortex ring [32,48].

We note that the presented technique of combining volume rendering of a vortex criterion with vortex core lines can be beneficial also in visualizations of a single time instant. As Kida and Miura [49] pointed out, vortex cores may overlap, making it hard or even impossible to distinguish between individual vortical structures, especially when using an isosurface representation. Combining volume rendering or transparent isosurfaces with the vortex core lines is a possible way of tackling this issue.

3. Conclusions

Suitable methods of computer graphics flow visualization for application on the vortex rope and possibly other phenomena were pursued in this paper. Common methods of vortex identification were tested to capture the vortex rope, as well as more elaborated ones: the ${\mathit{M}}_{\mathrm{Z}}$-criterion and the residual vorticity based on a triple decomposition of the velocity gradient tensor. The latter mentioned led to interesting results in terms of scaling with the operating point (i.e., with varying swirl number) and distinguishability of individual structures.

Hybrid methods were proposed to combine the advantages of the aforementioned criteria. A way of visualizing the whole 3D field of the finite-time Lyapunov exponent (FTLE), a method connected with Lagrangian coherent structures, was also proposed. The fact that it contains many fine structures makes its visualization challenging. Application of the knowledge from FTLE studies on the vortex rope is yet to be done and might help in an effort to understand its emergence and behavior. New ways to visualize the main features of flow with a vortex rope were also devised. The aim was to combine most of them into a single picture. A single series of such pictures for different operating points brings an overall view on the flow behavior and the mutual relation of the visualized features.

An important enabler of the techniques presented in this paper is volume rendering. This method of 3D data visualization is used abundantly in medicine, while, in applied fluid mechanics and specifically in hydraulic machinery, its applications are rather seldom. Just like isosurfaces, as the prevailing method used to visualize various structures in the flow field, it also has the ability of highlighting the structures by making unimportant parts of the domain transparent. Unlike opaque isosurfaces, however, the volume rendering by variable level of transparency retains the information regarding the distribution of the values of the quantity of interest in the 3D domain. Thanks to this, it can also be combined with vortex core lines, making it easier to distinguish between overlapping vortex cores.

All proposed methods become increasingly useful with growing trend toward hybrid URANS-LES and LES simulations, where abundance of coherent vortices arises and flow becomes visually unclear or even chaotic. It is precisely these methods that enable to extract the key features that play an important role in flow dynamics, such us the disintegration and decay of the vortex rope. Moreover, suitability of general-purpose visualization packages in the field of fluid mechanics, particularly in application to vortex dynamics, was illustrated on several examples. Smart application of these tools makes possible capturing more phenomena in one figure to highlight mutual interaction and interplay of several flow features. Thus, visual representation becomes more effective and helpful.

Computational data presented in this paper are part of extensive research aimed on investigation of vortex rope behavior and its control. The presented techniques of vortex rope identification and visualization will be used to thoroughly investigate its behavior in the whole operating range of our swirl generator. Proper understanding of fluid flow and its response to various flow control methods is generally essential in tackling modern challenges of reaching higher efficiency, wider operating range, lower noise, and other goals.