Characterising Single and Two-Phase Homogeneous Isotropic Turbulence with Stagnation Points

: It has been shown that, for dense, sub-Kolmogorov particles advected in a turbulent ﬂow, carrier phase properties can be reconstructed from the particles’ velocity ﬁeld. For that, the instantaneous particles’ velocity ﬁeld can be used to detect the stagnation points of the carrier phase. The Rice theorem can therefore be used, implying that the Taylor length is proportional to the mean distance between such stagnation points. As this model has been only tested for one-dimensional time signals, this work discusses if it can be applied to two-phase, three-dimensional ﬂows. We use direct numerical simulations with turbulent Reynolds numbers Re λ between 40 and 520 and study particle-laden ﬂows with a Stokes number of St = 0.5. We conﬁrm that for the carrier phase, the Taylor length is proportional to the mean distance between stagnation points with a proportionality coefﬁcient that depends weakly on Re λ . Then, we propose an interpolation scheme to reconstruct the stagnation points of the particles’ velocity ﬁeld. The results indicate that the Rice theorem cannot be applied in practice to two-phase three-dimensional turbulent ﬂows, as the clustering of stagnation points forms very dense structures that require a very large number of particles to accurately sample the ﬂow stagnation points.


Introduction
Turbulent flows laden with inertial particles are widely encountered in nature, playing a preeminent role in particles dispersion in the atmosphere, rain formation and marine snow sedimentation, among others [1,2]. They are also relevant for several industrial flows, such as fuel or coal combustion, fluidized beds reactors and separation techniques. One of the main challenges to characterizing these flows is the need to simultaneously resolve the particle positions and velocities and the flow velocity field at their scale [3,4]. All these configurations involve highly turbulent three-dimensional flows, which can be highly inhomogeneous and unsteady, where possible finite-size effects from particles may also be present. In this work we focus on a simplified case: homogeneous isotropic turbulent flows (HIT) laden with point-like inertial particles.
The stagnation points of velocity fields in turbulent flows present several relevant characteristics that can be used to gain further understanding of these systems. For instance, the zero-crossings of fluctuating one-dimensional velocity signals have been extensively studied, as they can be used to quantify the Taylor microscale λ of homogeneous isotropic turbulence via the Rice theorem [5][6][7]. As a consequence, these structures have been intensively studied over the last years, in works that cover the energy cascade of turbulence [8,9] and atmospheric flows [10,11], among others. They present several advantages, such as the zero-crossings of a velocity signal being robust when the flow is unsteady and/or the calibration of probes is not guaranteed. Furthermore, it has recently been shown that they can also be used to quantify the integral length scale L [12,13].
While most works focus on zero crossings, others have considered the case of stagnation points (STPS), defined as the set of velocity nulls satisfying v(x n ) = 0, where v is the fluid velocity field [14,15]. In particular, Goto and Vassilicos [16] generalized the Rice theorem, finding a relation between the number density of STPS and λ, with B a constant that may vary with Re λ due to the dependency of the non-Gaussianity of velocity derivatives with this parameter. While the study from Goto and Vassilicos [16] confirmed the validity of this theorem, it did not explore a sufficiently large range of Reynolds numbers based on the Taylor scale Re λ to report on the dependency B(Re λ ).
Although these studies concern single-phase turbulent flows, it has recently been proposed that the Rice theorem can be applied to particle-laden turbulent flows. The work by Mora et al. [17] developed an experimental method to estimate the carrier-flow turbulent kinetic energy dissipation rate ε in the presence of inertial sub-Kolmogorov particles at moderate Re λ . Its foundations rely on the unladen flow dissipation calculation using the Rice theorem, and the density of zero crossings n s . Moreover, the results from such a model apply, in principle, also to three-dimensional particle velocities depending on the simplified equation of motion, with v p the particle velocity and u(x p , t) the carrier's flow velocity evaluated at the particle's location x p , and τ p the particle viscous response (defined in the next section). This simplified model relies on two conditions: the diameter of the particles must be smaller than the Kolmogorov lengthscale of turbulence η, and their density must be much larger than the carrier's flow density. The Fourier transform of Equation (2) yields, v p =û iωτ p + 1 .
As a consequence, the particle field velocity is a low-pass filtered version of the carrier phase one, with a cut-off frequency of f c = τ −1 p /2π, or f c τ η = (2πSt) −1 . The cut-off frequency therefore depends on the Stokes number of inclusions, defined as St = τ p /τ η , with τ η = (ν/ε) 1/2 the flow Kolmogorov time scale (ν is the fluid kinematic viscosity). We can then deduce from Equation (3) that if the cut-off frequency f c is large enough to resolve the dissipation scales, n s should be recovered. Thus, it is possible to deduce the value of λ from the particles' velocities. As stated above, while this model has been developed for 1D signals and zero crossings, Equation (3) is already defined for three-dimensional velocities, and the zero-crossings number density can also be redefined as the stagnation points number density [16].
We can therefore conclude, in principle, that the generalized Rice theorem and the model from Mora et al. [17] can be combined to deduce the carrier phase value of λ using inertial particles. This rationale can also give access to other small-scale quantities, such as ε and η, among others. Beyond its fundamental interest, this could also be used to quantify the carrier phase properties in experiments of two-phase turbulent flows. Indeed, to resolve the carrier phase simultaneously with the inclusions velocities in such conditions is beyond the possibilities of current experimental techniques. Finally, to quantify these properties of the carrier phase would also help to detect the presence of two-way coupling between the inertial particles and the carrier flow.
The first two conditions for the applicability of this model are: (i) that Equation (3) holds, and (ii) that the cut-off frequency verifies the relation f c τ η = (2πSt) −1 > 10 −2 . This condition is proposed based on the fact that Vassilicos and collaborators [8,9] found that, when low-pass filtering zero crossings n s with cut-off frequencies at least one order of magnitude larger than the Kolmogorov length-scale, such low-pass filtered velocity records were still able to resolve the value of λ. The last condition, (iii) is to have enough particles to sample all the stagnation points present in the flow. This latter condition has the added constraint that stagnation points are known to form dense clusters [14], thus making the sampling of these points more difficult.
This work aims at studying the applicability of the aforementioned model for stagnation points. We use direct numerical simulations (DNSs) with random forcing, that avoid experimental errors that may contaminate the counting of stagnation points [7]. The results could be extended to experimental fields of particles advected by turbulent flows. We will first focus on verifying the applicability of the generalized Rice theorem (Equation (1)) to our DNSs, which span a wide range of Re λ that goes between 40 and 520. We will particularly focus on the dependency of B with Re λ , not discussed in previous works, and in verifying the applicability of Equation (1) to instantaneous velocity fields. Then, on a second part, we will use the method of Mora et al. [17] to study DNSs of two-phase flows with Re λ = 240, using up to 10 7 tracers (i.e., St = 0) that follow the streamlines of the flow, and inertial particles (St = 0.5) that evolve according to Equation (2). Our results indicate that the model may not apply unless an extremely large number of particles is injected, as the clusters of stagnation points require a very large spatial resolution (or particle densities) to be resolved.

Numerical Simulations
Our study was conducted using DNSs at five different values of Re λ using the GHOST code [18,19]. These simulations follow standard practices regarding their temporal integration, dealiasing procedures, and have an adequate spatial resolution of the smallest scales, i.e., κη 1 (where κ = N/3 is the maximum resolved wavenumber in Fourier space and N the linear spatial resolution [20]). The Kolmogorov lengthscale η is defined as η = (ν 3 /ε) 1/4 . Fully dealiased pseudospectral methods with second-order Runge-Kutta methods for the time stepping are used. The 3D simulation domain for all datasets has dimensions of 2π × 2π × 2π. All relevant parameters can be found in Table 1. Table 1. Relevant parameters from the DNS used in this study. N is the number of points in the DNS in one axis, such that N 3 is the total number of grid points in the simulation domain. L/(2π) is the integral lengthscale in units of the domain linear size 2π. η is the Kolmogorov dissipation scale. Re λ is the Reynolds number based on the Taylor microscale λ. "# snapshots" is the number of snapshots of the vector fields used for the analysis, and # STPS snps the averaged number of STPS (i.e., stagnation points) over the total number of snapshots. Numerical simulations solve the incompressible Navier-Stokes equations for the velocity u with a random solenoidal forcing f, where p = p/ρ (with p the pressure and ρ a uniform mass density), which is obtained from the incompressibility condition ∇ ∇ ∇ · u = 0. In Equation (4), Du/Dt = a is the Lagrangian acceleration of the fluid elements. We define the r.m.s. velocity as u = |u i | 2 1/2 (where u i is a Cartesian component of the velocity and Einstein notation is used), the Taylor microscale is λ = (15νu 2 /ε) 1/2 , and the integral scale is L = π/(2u 2 ) E(k)/k dk (where E(k) is the isotropic energy spectrum).
The solenoidal forcing f is given by a superposition of Fourier modes with random phases in the shell with wavenumber k = 1. A new random forcing was generated every 0.5 large-scale turnover time, and the forcing was linearly evolved from its previous state to the next state along this period of time. This results in a continuous and slowly evolving random forcing with a correlation time of 0.5 turnover times, which at the largest resolution considered has an integral scale L/(2π) ≈ 0.309, and which will be useful for simulations with inertial particles, as discussed below. The simulations also use the largest Reynolds number attainable at their spatial resolution, with κη ≈ 1 (see Table 1).
We use five numerical datasets, labeled in the following as "DNS-N", where N is the linear resolution of each dataset. The Taylor-based Reynolds number, Re λ = u λ/ν, spans more than one decade. We have Re λ ∈ [40, 520] for spatial resolutions of 64 3 , 128 3 , 256 3 , 512 3 , and 1024 3 grid points (see Figure 1). We stored enough snapshots of the vector fields to have adequate global statistics. For all datasets, we applied the method proposed by Haynes and collaborators [21][22][23] to compute the stagnation points. This method goes through each cell of the DNS domain and uses the velocity values of the eight cell's corners. If there is a change of sign in all the three velocity components, a local trilinear interpolation function is created with the corner's velocity values. Then the Newton method is used to find any velocity nulls within the cell. If there is no change of sign in any of the three velocity components then no velocity nulls should be contained in the cell for a well resolved DNS. Both elliptic and hyperbolic stagnation points were considered. More details about this method can be found in references [14,21]. For DNS-512 we also have data of tracers and inertial point particles without gravity. Particles are integrated following Equation (2), which can be written as: These equations are integrated with a high-order Runge-Kutta method to evolve in time and a high-order three-dimensional spatial spline interpolation to estimate the fluid velocity u(x p ) at the particle position (see [24,25] for details). Simulations with particles are conducted as follows: first, a DNS of the Eulerian flow is conducted without particles, until a turbulent steady state is reached. Then, particles are injected with a uniform random distribution in space, and with the same initial velocity as the velocity of the fluid element at the particle position. Particles are integrated for several turnover times (in the case of tracers) or for several particle relaxation times (in the case of inertial particles) before data starts to be collected for the analysis.
In order to apply the method described in the introduction, the particles' velocities v p were interpolated using the 'griddata' function in the interpolate module of the SciPy library. Particles' velocities were thus interpolated on the DNS Eulerian grid points with a linear interpolation. For points outside the particles' position range, on the boundaries of the DNS domain, the nearest method was employed, obtaining a synthetic velocity field from the particles for each snapshot. We therefore proceeded to apply the method from Haynes [21] and collaborators here too, to detect the stagnation points.
As discussed above, we used two types of particles: inertial particles with St = 0.5 and tracers with St = 0. While tracers are expected to sample the flow uniformly, the former have been reported to cluster [14]. For each type of particle, two DNSs were run, one with 10 6 particles and another with 10 7 . This will allow us in the next sections to analyse the influence of the number of particles in the convergence of the Rice theorem, and to verify if it actually applies for our datasets (we remind the reader that particles are injected only for the N 3 = 512 3 run, i.e., with ≈ 1.3 × 10 8 Eulerian grid points). In the following, our datasets with particles will be labelled as: As each run was performed independently, and given that we used random forcing, the temporal evolution of the single-phase flow is not expected to be identical for all datasets.

Validation of the Generalized Rice Theorem in Single-Phase Turbulent Flows
We first verify the validity of the generalized Rice theorem for our DNSs. To this end, in Figure 2 we plot the prefactor B in the Rice theorem (defined in Equation (1)) as a function of Re λ . Figure 2a is deduced for each instantaneous Eulerian velocity field from the DNSs at N = 512, with Re λ also computed instantaneously for each snapshot (n s f denotes the density of stagnation points in the Eulerian fluid velocity). On the other hand, Figure 2b shows the value of B for all resolutions N in the DNSs, and averaged over all snapshots corresponding to that run as detailed in Table 1 (Re λ is also computed from the averaged characteristic quantities). Our results are in good agreement with the generalized Rice theorem, as we find B is of order unity and has a small dependency with Re λ . Furthermore, our results are consistent with the study from Goto and Vassilicos [16] on a similar flow. Note that while this study focused on the range Re λ ∈ [60, 150], in Figure 2b we extend the validity of the theorem to Re λ ∈ [40, 520]. Our DNSs also show that B is a slowly decreasing function with Re λ , and therefore presents the opposite trend as the one found for zero-crossings of the fluctuating velocity [8] (but consistent with the constant value found in [9], where the same flow was studied for a large range of Re λ ). This result is surprising, as it suggests that the contribution of small-scale intermittency effects decreases when Re λ is increased.
The good collapse of all fields shown in Figure 2a,b suggests that the generalized Rice theorem is valid not only when averaging in time Eulerian fields, but also for instantaneous realisations. To confirm this feature, Figure 3a shows the temporal evolution of n −1/3 s f and λ for DNS-512-1 and DNS-512-2 for the Eulerian fluid velocity field. It can be observed that both curves present almost identical trends, with values of B that remain almost constant (see Figure 3b). We therefore conclude that the generalized Rice theorem applies to our datasets, for a wide range of Re λ as well as for instantaneous velocity fields.

Validation of the Rice Theorem in Two-Phase Turbulent Flows
We now proceed to study the applicability of the generalized Rice theorem to particleladen flows. It can be easily seen that for both St = 0 and St = 0.5 the condition f c τ η = (2πSt) −1 > 10 −2 holds. Additionally, the applicability of Equation (2) is trivially valid in our case, as particles evolve according to it in the simulations. Figure 4 shows the reconstructed (i.e., interpolated) x velocity field component using 10 7 particles with St = 0 or St = 0.5, compared to the actual Eulerian flow velocity component. Fields are very similar, although not identical. The squared point-wise differences between the three velocity fields are presented in Figure 5. Discrepancies between the flow field and inertial particles are expected (Figures 4c and 5b), as Equation (3) implies that the STPS may be preserved but not the velocity values elsewhere. Furthermore, this is also expected when comparing the flow velocity and the tracers (Figure 4a,b respectively, or see Figure 5a), as the latter are also expected to preserve the null points but, depending on the number of particles present in the flow, could result in a coarse-grained reconstruction of the former. Moreover, as expected, the tracers' reconstructed field shows fewer differences with the flow field than the inertial particles (see Figure 5a,b).
We can therefore study the validity of Equation (1) in our DNSs. A first test is to compare the number of STPS detected for the carrier phase and for the interpolated (tracer or inertial) particle velocity fields (Figure 6a). Surprisingly, we see that for any of the particle sets, a smaller number of STPS are detected, and even when injecting 10 7 tracers we have 30% fewer STPS in the interpolated field (Figure 6a). This surprising result points towards the inapplicability of the model from Mora et al. [17] to our datasets. This is confirmed when comparing n −1/3 sp for the tracers' interpolated field and λ in Figure 6b, as we find values of B to always be larger than those found for the Eulerian flow velocity in Figure 3b. Nevertheless, Figure 6b suggests that, while the generalized Rice theorem does not apply exactly, some trends are still recovered. In other words, a calibrated value of B(St, N part , . . . ) (where N part is the total number of particles) could be used if data are available from, e.g., numerical simulations.  Black and red circles correspond to STPS detected from the tracers and inertial particles fields, respectively. As the amount of STPS in a 2D slice can be small, in all panels we show all STPS in a slice with 6 grid points in x (i.e., the STPS detected in x ∈ π/2 ± 6π/512). As the only condition that may be violated is to have enough particles to sample all stagnation points in the flow, we will now analyse this hypothesis. For all N = 512 DNSs, we find that the carrier phase has values of n −1/3 s f of around 0.3. This implies that we have around 30 stagnation points per unit of volume (with a total volume of (2π) 3 in our DNSs). Conversely, we have 4 × 10 3 and 4 × 10 4 (inertial or tracer) particles per unit of volume when injecting, respectively, 10 6 and 10 7 inclusions.
These densities imply that, in principle, all stagnation points should be resolved by the interpolated fields. However, as we discussed in the introduction, this consideration does not take into account the clustering of stagnation points. Indeed, Figure 7a,b shows the presence of strong clusters of STPS in the flow. Using similar DNSs, a previous work [14] showed that stagnation points form very dense clusters, and that clustering increases significantly with Re λ (as shown in Figure 7c). Such clustering (against, e.g., an homogeneous spatial distribution of points) implies that to resolve all stagnation points would require injecting many more particles, as dense regions of stagnation points also require larger densities of particles to be resolved with our interpolation scheme. This can also be seen in Figure 4, where the interpolated fields for both tracers and inertial particles do not recover all stagnation points in dense regions. We remark that tracers are instead distributed homogeneously in space, while inertial particles are also known to form clusters in turbulent flows, but the positions of such clusters are not directly related to the positions of those formed by the stagnation points [14]. Besides this effect, another important effect that can explain why the particles see fewer zeros than the number of STPS in the Eulerian field is associated with the stability of stagnation points in 3D flows [26]. Even in the simpler 2D case, most instantaneous stagnation points can be classified into elliptic or hyperbolic types, depending on their local stability. Tracers and inertial particles can be expected to spend longer times around elliptic points, whereas hyperbolic points should quickly push nearby tracers and inertial particles along their unstable manifolds. In 3D, the possible topologies of velocity field nulls are more complex, but any 3D stagnation point with an unstable manifold should have the same effect. As a result, STPS are sampled differently depending on their topology, and some STPS will be less sampled than others. This can be another important reason behind the lower number of detected zeros from the 3D particles' fields.

Conclusions
Throughout this work we used DNS data to study the applicability of the generalized Rice theorem to single and two-phase flows. Our results can be summarised as follows:

•
We verified the validity of the generalized Rice theorem for our dataset, which covers the range Re λ ∈ [40, 520]. Furthermore, we showed that the prefactor B in Equation (1), that quantifies the non-Gaussianity of velocity derivatives, has a weak dependence with Re λ , and that it tends to decrease when this parameter is increased.
• Furthermore, we showed that the generalized Rice theorem applies for time-averaged three-dimensional velocity fields, but also for instantaneous realizations. • We proposed an interpolation scheme to reconstruct the stagnation points using the particles' velocity field. Our results indicate that the Rice theorem cannot be applied in practice to two-phase three-dimensional turbulent flows, as the clustering of stagnation points forms very dense structures that require a very large number of particles to accurately sample the flow stagnation points. Even with 10 7 tracers or inertial particles, we did not manage to apply the Rice theorem satisfactorily. • We find that this lack of resolution of stagnation points is consistent with the strong clustering of STPS, as it implies the presence of very dense regions of these points, which require the injection of a very high number density of particles to be resolved. Another possible explanation for the lower number of STPS detected with the particles' velocity field is the local stability of 3D STPS with unstable manifolds. • While the number of the carrier phase STPS is always larger than the one obtained when using the interpolation scheme proposed here, we do find that they evolve over time following similar trends. This feature requires further study to be validated.
In conclusion, our study suggests that the generalized Rice theorem and the rationale from Mora et al. [17] cannot be used in a practical way to reconstruct the carrier phase from particles' measurements in turbulent two-phase flows. Its application would require a number of particles that would make such study extremely hard using modern experimental techniques.