Temporal vortex tracking in isotropic turbulence

An algorithm to allow the temporal vortex tracking of the intense vorticity structures is introduced and applied to direct numerical simulations (DNS) of homogeneous isotropic turbulence (HIT). The algorithm allows the identification of the individual life events of the worms, including events of splits into several structures and of merging between different structures. Results obtained from a relatively small DNS, with Taylor based Reynolds number equal to Reλ = 54, for a simulation lasting more than 6 integral time scales, shows that the complement of the CDF of the life time of the IVS is well approximated by an exponential distribution function. The mean life time of the IVS which is approximately 3.6τη , where τη is the Kolmogorov time, and the birth rate of these structures is close to 4.5/τη (number of new structures per Kolmogorov time). A population model is built to predict the number of existing structures at any time, which gives a number of existing IVS equal to ≈ 16, which agrees well with the results from the DNS.


Introduction
It is well know that intense vorticity structures (IVS) or 'worms' abound in virtually all turbulent flows.These structure are important for the study of the viscous dissipation mechanism and to internal intermittency [1], and their study is also motivated by their role in the transport, mixing and diffusing at the small scales of motion [2,3].
The IVS are typically defined as regions of particularly strong concentrated vorticity and are often mentioned as having a life time which is big compared with the characteristic time scale of the flow [4].Several methods have been employed to identify these structures as discussed in [5,2].These include the use of iso-surfaces of vorticity, which has been mostly used in homogeneous and isotropic turbulence [6,7,8,9], and quantities derived from the velocity gradient invariants [2,10,11].However, the criteria to define the IVS that we will use here is the one proposed by [12,13] who defined these structures as the regions of particularly strong vorticity, defined by flow points with the highest enstrophy representing 1% of the total, [12,13].
Although the IVS have been primarily studied in isotropic turbulence e.g.[6,7,8,9,14,15,16,12,5,13,17,18] other flow types have also been studied such as mixing layers [11], channel flows [10] and jets [19,3].It turns out that many of the characteristics of the IVS observed in isotropic turbulence, such as the size of the vortex core radius, are also observed in other flow types.
For instance, it has been shown that the vortex core radius of the worms is of the order of the Kolmogorov micro-scale ⟨R ivs ⟩/η ≈ 4 − 5, while the mean tangential velocity is close to the root-mean-square velocity ⟨U ivs ⟩/u ′ ≈ 1, both in isotropic turbulence [13,17] and also in mixing layers [11], channel flows [10].Recently, Ghira et al. [17] have shown that the length of the IVS also scales with the Kolmogorov time, with ⟨L ivs ⟩/η ≈ 60, in isotropic turbulence.
Despite the large number of these studies, a detailed temporal analysis of the life events and life times of these structures is still lacking.In the present work we develop an algorithm to allow the temporal vortex tracking of the intense vorticity structures and we apply this to a single simulation of isotropic turbulence with relatively small Reynolds number as illustration.Larger simulations using the same procedure will be reported elsewhere.The algorithm allows the identification of the individual life events of the worms, including events of 'splits' into several structures and of 'merging' between different structures.This allows us to rigorously determine the life times of the IVS in isotropic turbulence.
This paper is organised as follows.The next section (Section II) describes the DNS used in the present work.Section III describes the temporal tracking algorithm to compute the IVS lives, while Section IV analyses the results obtained in the present work.The work ends in section V with an overview of the main results and conclusions.

Direct numerical simulations of homogeneous isotropic turbulence
For the purposes of this work we employ a single simulation of forced Homogeneous Isotropic Turbulence (HIT) obtained by solving the incompressible Navier-Stokes equations in a triple periodic domain using a pseudo-spectral solver.We use the Alvelius' forcing scheme [20] and let the turbulence evolve until statistically steady conditions are reached.Statistical stationarity is obtained by imposing a power input through the imposed forcing that balances the viscous dissipation rate, P = ϵ, where P = f i u i is the power of the input forcing and ϵ = 2νs ij s ij , is the viscous dissipation rate, while s ij = (∂u i /∂x j + ∂u j /∂x i )/2 is the rate-of-strain tensor.
Table 1 summarises the physical and computational parameters of this simulation in the statistically stationary conditions used in the subsequent analysis.The simulation consists on a small HIT case using N 3 = 128 3 collocation points, for a Taylor based Reynolds number of Re λ = 54 and a resolution of k max η = 2.0.The analysis is done using a total of N f = 1400 instantaneous fields and lasts for close to 7 integral time scales and more then 50 Kolmogorov time scales.
Figure 1 shows the evolution of the kinetic energy and dissipation rate normalised by their initial values at statistically steady conditions, where the snapshots analyzed by the Wormtracker were taken.The figure shows the typical variability present in similar (statistically stationary) simulations of HIT.As expected the viscous dissipation exhibits proportionally greater variations than the kinetic energy, however the instantaneous value of ϵ never exceeds 10% its mean value during the simulated times.

Wormtracker
As described in the introduction the intense vorticity structures (IVS) are defined by flow points (from the DNS grid) with the highest vorticity magnitude covering 1% of the total flow domain.
The algorithm used to detect and quantify these structures ('wormtracker') is based in the one from Jiménez et al. [12,13]), with some improvements already described in Ghira et al. [17].The only aspects of the wormtracker that are relevant to the present work are related to the detection, and automatic assignment/numbering of these structures.Other aspects of the wormtracker, such as the computation of the associated kinematic and geometrical aspects (e.g the IVS' tangential velocity and vortex radius), are detailed in Ghira et al. [17].
In order to detect these structures an histogram of vorticity magnitude is computed for each instantaneous field of the simulation (N f fields in total) in order to obtain the reference vorticity magnitude threshold ω ivs corresponding to the 1% highest vorticity magnitude.Any point for which the local vorticity magnitude ω = (ω i ω i ) 1/2 (ω i is the vorticity vector) is greater than ω > ω ivs is assigned to an individual structure.The workmtracker uses the direction of the local vorticity vector of each individual structure to build the full axis of the structure and all the IVS are numbered, starting (number one) with the structure with the strongest vorticity, to the structure with smallest vorticity.As in [12,13,17] structures with less than 20 axis points (from the ones in the DNS grid that have been assigned to an IVS axis) are discarded from the statistical sample, however, unlike in previous works, the wormtracker algorithm runs together with the Navier-Stokes solver so that there is no need to store and analyse each one of the individual instantaneous fields in subsequent post-processing, i.e. the wormtracking is done during the actual simulation.This proved to be essential to allow the time tracking of each one 5th Madrid Turbulence Workshop Journal of Physics: Conference Series 2753 (2024) 012011 of the individual structures during very long times, which is essential to make a proper temporal analysis of the IVS.

Connectivity
In order to connect structures of two contiguous snapshots, an algorithm was developed based on spacial overlapping principles.
Given a snapshot at time t s , all identified structures are assigned an index/number, with the number 1 being assigned to the structure having the most intense vorticity point, followed by 2 for the second most intense structure, and so on, until the last, less intense vortex axis has been numbered.Next, using a scalar variable, we assign the same index/number to this scalar in all mesh points within a distance of one ∆x (mesh spacing) from each point along the filament in every direction.The result is a scalar variable identifying in space the position of a filament and its representative cloud.For instance, the scalar field will display the value of 1 to all the points of the cloud surrounding the most intense axis (the axis containing the most intense vorticity point).This choice is based on the assumption that a filament is not expected to move away, appreciably, from its cloud within one timestep ∆t due to the Courrant Number restriction.
For the next snapshot, at time t s + ∆t, a new set of structures is identified and a new indexing/axis numbering takes place, again starting from the axis of the most intense vortex structure.However, by overlapping the new structures with respect to the scalar field of clouds from the previous snapshot, we are able to connect and refer back in time this new indexing.In other words, we can assign a global index/numbering that identifies unambiguously each structure from its 'birth' to its 'death' and relate it to the local indexing as function of time using all the existing/instantaneous snapshots from a given simulation.
This procedure allows us to have a spacial recognition about the positions of the structures in subsequent timesteps.Nevertheless, we still need to build some criteria to make comparisons.We start by defining the overlap operator O ij , which represents (in %) for a given cloud i how much of a worm j belongs to it (where only the points belonging to some cloud are taken into account).For instance, consider a given worm j which has n points along its filament, and suppose that from these n axis points, a number of x points belong to cloud a, y axis points to cloud b and z axis points do not overlap with any other existing cloud (n = x + y + z).In this case the the matrix elements of the overlap operators O aj and O bj are defined by O aj = x/(n−z) and O bj = y/(n − z), respectively.We then define Complete Superposition when O ij > 90% and Incomplete Superposition if 10% ≤ O ij ≤ 90%, where by superposition we mean that worm j is overlapping cloud i.
Using these definitions and metrics we can now look for the most simple events in the temporal evolution of a worm, namely, i) the continuation of the life of the worm (Same Worm), ii) Pure Splitting, iii) Pure Merging, iv) Death, and v) Birth.Same Worm events occurs whenever there is a Complete Superposition of a cloud by one structure.A Pure Splitting means that a cloud has more than one structure with Complete Superpositions.A pure merging happens when more than one cloud merge into one structure with Incomplete Superpositions on them and, at the same time, all the clouds involved must satisfy the same merging criteria.A Death is detected if one cloud does not have any structure overlapping.Finally, a Birth occurs because a structure is not overlapping any existing cloud.These criteria are summarised in Table 2.
As mentioned before, the above connectivity algorithm allows one to build a global worm index using every one of the local indexes obtained for each time step/instant of the simulation.With this procedure all statistics taken at each instant of time can be traced back to each specific individual structure, as a time evolving entity.The most basic entity arising from the connectivity algorithm is the Lifetime Map, which is shown in Figure 2 for the present simulation, whose parameters have been presented in Table 1. This > 0 Pure Merging (all clouds must satisfy) = 0 = 0 Death Note that the present analysis starts from the statistically stationary state, which explains the existence of about ≈ 20 axis already at ∆t = 0.The first detected splitting event is shown in Figure 3, in exactly the two (contiguous) times when that event was recorded.There we can see (Figure 3 exists with the corresponding cloud.In the very next instant (Figure 3(b)) we see two different structures instead.The connection is made by comparing those two structures with the cloud presented in Figure 3(a).By the end of the second timestep, the scalar variable with the cloud information is reconstructed to allow the continuation of the analysis in the next time step.

Results
In this section we will use the temporal vortex tracking algorithm described to estimate the lifetime and the population of the IVS.Analysis of the other types of events defined above, as well as the computation of these quantities at higher Reynolds numbers will be reported elsewhere.

The life time of the intense vorticity structures
From the data gathered, we can compute a simple estimation for the Cumulative Distribution Function (CDF) of the worms' lifetime, P (T /τ η < t/τ η ), where T is the lifetime random variable, t is the measured lifetime and τ η the Kolmogorov's timescale.In this process we discard the 5th Madrid Turbulence Workshop Journal of Physics: Conference Series 2753 (2024) 012011 Figure 4. Complementary probability of the CDF of the worms' lifetimes, where data for lifetimes smaller than a given value t 0 have been removed.T is the random variable describing the lifetime t and τ η is the Kolmogorov's timescale.In the CDF for t 0 /τ η = 0 all the lifetimes have been used.
contribution from structures that are 'alive' both at t = 0 and t = t f (final time), because we do not know their past and future history, and their inclusion would affect the estimate.This is justified because the simulation used here was run for a very long time, with more than 6 integral time scales (see table 1).Since the CDF is computed for each known (measured) value of the lifetime data the use of thresholds in the detection method described above may lead to some numerical artifacts in the results.We expect to see effects from these artifacts especially for very small timescales.Therefore, we investigate their effect by recomputing the CDF eliminating structures living less than a given time t 0 .Figure 4 shows 5 CDFs obtained by removing from the statistical sample lifetimes smaller than times between 0 ≤ t 0 /τ η ≤ 5.In the CDF for t 0 /τ η = 0 no removal was done i.e. all the lifetimes have been used (note that in this figure the vertical axis represents not the CDF but the complementary probability of the CDF).
The figure shows that using different cutoffs t 0 /τ η affects the CDFs shape, however all CDFs converge to virtually the same function for t 0 = 0.5τ η .We therefore proceed with our analysis by eliminating all the structures displaying lifetimes smaller than t 0 /τ η = 0.5, as it is likely that events lasting less than half a Kolmogorov time are included in artificial numerical 'noise'.This is done with figure 5 which displays the CDF for t 0 ≤ 1τ η , which we label as the 'true' probability function for a structure to live more than a certain time t.Noting that this function approximates well to a line with a constant slope, we conclude that the CDF is well approximated by an exponential function.This function, whose parameters are indicated in the enclosed label, is also added in the figure.From this we can now estimate the mean life time of the IVS by computing the reciprocal of 0.32 (from the fitted curve) and summing it to the cutoff chosen value.We obtain τ /τ η = 3.61 for the average life time of these structures.This value is much Figure 5. CDF of the worms' lifetimes.T is the random variable describing the lifetime t and τ η is the Kolmogorov's timescale.t 0 is the removed small time intervale that has been eliminated from the data.smaller than the integral time scale usually claimed as being the typical life time of the IVS [21].Our explanation, consubstantiated in several animations, is that even though many of these long life times structures do exit (and maybe these are the structures that we retain in our mind from movies and visualisations) they represent relatively rare events.

The population of the intense vorticity structures
Using the life time data available from the simulation it is possible to build a population model, to estimate the number of existing 'alive' IVS, as well as the number of Births and Deaths of these structures for a given simulated time.We postulate that the population model follows a well know law used in particle physics.In short, we assume that the number of individuals since a given time N (t) is governed by the following differential equation, where τ is the mean lifetime of the individuals (worms) and N B the accumulated number of births, (thus dN B /dt represents their birth rate).Assuming a constant value of this rate dN B /dt = C we write, where C is the birth rate constant.The solution of this equation is, where N 0 = N (0) is the initial number of individuals/worms, at the initial time when the time tracking starts.With the filtered data used in our analysis, where we have discarded the structures at t = 0 (see above), we have N 0 = 0 and therefore we can estimate the birth rate constant C by fitting our experimental curve of accumulated number of births to the model N B = Ct.Moreover, for sufficiently long simulating times, the model predicts a saturated number of individuals N (∞) = Cτ existing at a given (long) simulated time.Table 3 summarises all the results, including the mean lifetime estimated in the as previous section.We see that for this relatively small simulation we have about ≈ 16 worms alive at any given instant, and the normalised birth rate constant is Cτ η ≈ 4.5.These values are well supported from the DNS data.Finally, we again use the lifetime map data to compute a series of statistics such as the number of births and deaths of worms, and the number merges and splitting events.Table 4 gives those numbers for the final simulated time.The number of existing structures at t = 0 and t = t f is also shown.Due to statistically steady conditions of our simulation, we expect the number of IVS to remain approximately constant and thus comparable to what our model predicts at saturated conditions.This is indeed the case as shown in the table.We start/end the analysis with approximately the same number of structures, 24 and 19, respectively.During the simulations we had approximately 900 births and deaths, and the number of splits (48) is slightly larger than the number of merger (16).
Table 4. Summary of some statistics for the data: Number of colocation points (N 3 ); Taylor micro-scale based Reynolds number (Re λ ); Initial number of individuals (N (0)); Final number of individuals (N (t f )); Final accumulated number of births (N B (t f )); Final accumulated number of deaths (N D (t f )); Final accumulated number of merges (N M (t f )); Final accumulated number of splittings (N S (t f )).

Conclusions
Intense vorticity structures (IVS) or 'worms' are known to populate virtually all turbulent flows and a detailed temporal analysis of the life events and life times of these structures is still lacking.
In the present work an algorithm to allow the temporal vortex tracking of the intense vorticity structures is introduced and applied to direct numerical simulations (DNS) of homogeneous isotropic turbulence (HIT).
The algorithm allows the identification of the individual life events of the worms, including events of splits into several structures and of merging between different structures.Results obtained from a relatively small DNS, with Taylor based Reynolds number equal to Re λ = 54, for a simulation lasting more than 6 integral time scales, shows that the complement of the CDF of the life time of the IVS is well approximated by an exponential function.
This allows one to compute the mean life time of the IVS which is approximately 3.6τ η , where τ η is the Kolmogorov time, while the birth rate of these structures is close to 4.5/τ η (number of new structures per Kolmogorov time).A population model is built to predict the number of existing structures at any time instant giving a number of existing IVS equal to 16, that agrees well with the results from the DNS.
A more involved discussion of these events, including more advanced analysis methods and much higher Reynolds numbers will be published elsewhere.

Figure 1 .
Figure 1.Evolution of the kinetic energy (K) and dissipation rate (ϵ) normalized by their initial values K 0 and ϵ 0 , respectively.During this time window the Wormtracker analyzes the instantaneous fields at each timestep of temporal integration.

Figure 2 .
Figure 2. Zoom of the lifetime map showing the history of several worms.Every horizontal line represents a worm that is 'alive'.The vertical axis shows the global index/number assigned to each structure while the horizontal axis represents the elapsed time in timestep units.The red lines show the connection among branches arising from splitting events.

Figure 3 .
Figure 3. Visualisation of the first splitting, showing (a) the structure undergoing that splitting and (b) the resulting branches for one IVS located at one of the corners of the computational domain (the total size of the computational domain is much larger).Note that in (a) the cloud is also shown for the respective filament.That cloud is compared with the positions occupied by the two structures appearing in (b), thus identifying this chain of events.In (b) we also see the two clouds surrounding each new filament.The process of cloud reconstruction takes place at the end of each snapshot analysis for further time tracking.

Table 2 .
Classification criteria based on the number of structures contributing to a specific range of values of the overlap operator O ij .This table should be interpreted with respect to each one of the existing clouds i.Note that a cloud with given merging criteria will only be considered as merged if all the involved clouds satisfy the same criteria.
figure shows a zoom of the lifetime map showing the history of several worms.Every

Table 3 .
Estimates from DNS data and from the population model: Number of collocation points (N 3 ); Taylor micro-scale based Reynolds number (Re λ ); Mean lifetime normalized by the Kolmogorov's timescale (τ /τ η ); Birth rate normalized by the Kolmogorov's timescale (Cτ η ); Number of individuals predicted by the population model at saturated conditions (N (∞) = Cτ ).