Structure and dynamics of small-scale turbulence in vaporizing two-phase flows

Improving our fundamental understanding of multiphase turbulent flows will be beneficial for analyses of a wide range of industrial and geophysical processes. Herein, we investigate the topology of the local flow in vaporizing forced homogeneous isotropic turbulent two-phase flows. The invariants of the velocity-gradient, rate-of-strain, rate-of-rotation tensors, and scalar gradient were computed and conditioned for different distances from the liquid–gas surface. A Schur decomposition of the velocity gradient tensor into a normal and non-normal parts was undertaken to supplement the classical double decomposition into rotation and strain tensors. Using direct numerical simulations results, we show that the joint probability density functions of the second and third invariants have classical shapes in all carrier-gas regions but gradually change as they approach the carrier-liquid interface. Near the carrier-liquid interface, the distributions of the invariants are remarkably similar to those found in the viscous sublayer of turbulent wall-bounded flows. Furthermore, the alignment of both vorticity and scalar gradient with the strain-rate field changes spatially such that its universal behaviour occurs far from the liquid–gas interface. We found also that the non-normal effects of the velocity gradient tensor play a crucial role in explaining the preferred alignment.


Scientific Reports
| (2021) 11:15242 | https://doi.org/10.1038/s41598-021-94334-x www.nature.com/scientificreports/ of turbulent primary atomization and mixing. Using the same numerical parameters and physical properties, Bouali et al. 9 studied the effect of vaporization on scalar dissipation rate (SDR) statistics using a passive scalar to represent evaporation and mixing processes in two-phase turbulent flows. Hasslberger et al. 10 investigated the topological features of turbulent fluid flows in an atomizing liquid jet by the probability distribution of topology volume fractions and size distributions. Recently, the methodology of studying local topology in the context of multiphase flows has been applied to analyse flow structures in turbulence generated by rising bubbles 11 and in decaying droplet-laden isotropic turbulence 12 .
The aim of the present work is to build upon these studies by characterizing flow structures in the vicinity of evaporation fronts in a forced HIT, using classical Hermitian/skew-Hermitian decompositions of the velocity gradient tensor and the recently proposed Schur decomposition approach introduced by Keylock 13,14 , which splits the velocity gradient tensor into normal (characterized by the eigenvalues) and a non-normal (characterizing the tensor asymmetries) contributions. The main purpose of this study is to i) improve the fundamental understanding of the mechanisms of evaporating droplets and turbulence interaction, and the one of the most useful tools is the investigation of the local topology through the use of Schur and Cauchy decomposition and ii) provide insight regarding the limitations of the traditional Cauchy decomposition. Insights obtained in this work are expected to be beneficial for the development of subgrid-scale models and data-driven modelling methodologies.

Velocity gradient tensor decomposition
Cauchy-Stokes decomposition. The general pattern of the three-dimensional streamlines in the immediate vicinity of an observer moving at any point within the flow field can be determined by considering at the nature of the eigenvalues, ( 1 , 2 , 3 ) , of the velocity gradient tensor A ij = ∂u i /∂x j evaluated at that point, from which the symmetric strain-rate S A and skew-symmetric rotation-rate A tensors are obtained as S A ij = A ij + A ji /2 and � A ij = A ij − A ji /2 , respectively. Ordering i as follows, 1 > 2 > 3 , and knowing that for incompressible flows 1 + 2 + 3 = 0 , the first eigenvector ( e 1 ) is associated with the extensive eigenvalue ( 1 > 0 ), the third eigenvector ( e 3 ) with the compressive eigenvalue ( 3 < 0 ) and the second, intermediate, eigenvector ( e 2 ) can be either compressive or extensive. The three invariants P A , Q A , and R A of A ij in incompressible flow are expressed as: Similarly, the second and third invariants of S A are defined by its characteristic equation as follows: and The only invariant of A is: where ω A i = ε ijk A kj is the vorticity field. The invariants defined above are usually investigated in joint PDF s combining two invariants. The most common combinations consist on the maps of (Q A , R A ) , (Q S A , R S A ) and (Q � A , −Q S A ) . Figure 1 summarizes sketches of each of these maps, indicating the corresponding physical (1) (2) (4) Figure 1. Sketch of the invariant maps of (Q A , R A ) , (Q S A , R S A ) and (Q � A , −Q S A ) showing the physical and topological features related to each zone. In (a), SFS stands for stable-focus-stretching topology, UFS for unstable-focus-stretching topology, SNSS for stable node-saddle-saddle topology, and UNSS for unstable-nodesaddle-saddle topology.
Scientific Reports | (2021) 11:15242 | https://doi.org/10.1038/s41598-021-94334-x www.nature.com/scientificreports/ meaning of each particular location. The invariants R A and Q A are of great importance from a topological point of view, since the sign of the discriminant function A = 27R A,2 + 4Q A,3 , which characterizes a different streamline flow pattern, is dependant on their magnitude. In the (Q A , R A ) space, we distinguish four non-degenerate topology types allowing us to infer the relationships between local flow topology and physical mechanisms,e.g., enstrophy production and dissipation rate production (see Fig. 1a). The (Q S A , R S A ) map is a useful tool for analysing the geometry of local strain in the fluid. In particular, since Q S A = −ε 0 /4ν , large negative values of Q S A are associated with regions of intense kinetic energy dissipation (see Fig. 1b). Finally, the (Q � A , −Q S A ) map is particularly useful to analyse topology associated with the dissipation of kinetic energy (see Fig. 1c). For a detailed description of each zone in the three maps, the reader is referred to Blackburn et al. 15 .
Schur decomposition. The complex Schur decomposition applied to the velocity gradient tensor A ∈ C n×n , is given by where * denotes the Hermitian of a complex matrix, U ∈ C 3×3 is a unitary matrix (i.e., UU * = I ), and T ∈ C 3×3 is upper-triangular. Thus, the characteristic polynomials of A and T are the same and the roots of the latter are simply the diagonal elements of T . This decomposition is valid for any complex or real matrix. The matrix T in (5) is given by: where L is a diagonal matrix whose elements correspond to the eigenvalues of A , ( 1 , 2 , 3 ) , and is an upper triangular matrix that represents the non-normal part of A (characterizing the matrix asymmetries). The velocity gradient tensor A can be decomposed into its normal, B , and non-normal, C , contributions 13 as follows: The linear algebra tool used to accomplish this is the complex Schur transform (5). To ensure the uniqueness of the decomposition of A , the ordering of the eigenvalues of A is fixed in an increasing order, i.e., 11 (resp. 33 ) represents the most negative (resp. the most positive) eigenvalue of A . Applying transformation (5) to A obtain the following expressions: Thus, we note that B acts locally and is related to dynamics driven by the eigenvalues of A , whereas C is associated with vortical structures and non-local effects in the flow field. In order words, this decomposition explicitly gives i) information related to the eigenvalues of A , which are associated with local dynamics, and ii) the nonsymmetric structure in A , which is induced by non-local effects. This information is not included in the conventional Cauchy-Stokes decomposition of A . Note that non-normality could be measured by either �AA * − A * A� or , in which �φ� = tr φφ * ≡ 3 i,j=1 |φ ij | 2 is the Frobenius norm of a generic quantity φ . The following concepts and formula, which associate the Schur decomposition components with the invariants of A , are detailed in Keylock 13 . As a second step, we consider the strain-rate and rotation-rate tensors of B and C . The resulting expressions are: where In this case, the Frobenius norm of the strain and rotation tensors are as follows: Since the eigenvalues of A and B are identical, the invariants of tensor A are identical to those of B . Therefore, we obtain: and where E p φ = 1 4 ω i S � ij ω j is the enstrophy production rate for the generic quantity φ . Taking into consideration the above identities, we note that the dynamics of A are dictated only by the normal part, B . The second and  www.nature.com/scientificreports/ third invariants contain non-normal contributions when �C� � = 0 . Two new terms arise in both the strain and enstrophy productions for the third invariant, i.e., The strain production is the result of the sum of the normal strain production, the interaction production, and the non-normal production, while the enstrophy production is the sum of the normal enstrophy production, the interaction production, and the non-normal production.

Methods
Governing equations. In the following, the material properties of both phases, i.e., density and viscosity, are assumed to be constant. Subscripts ℓ and g denote the liquid and gas phases, respectively. In a computational cell, the local fraction of liquid volume is denoted α . The governing equations of the two-component incompressible flow model consisting of the three-dimensional Navier-Stokes equations for the mixture velocity u and pressure p are as follows: and the advection equation for volume fraction α is: where the mixture density ρ = (1 − α)ρ g + αρ ℓ and dynamic viscosity µ = (1 − α)µ g + αµ ℓ are computed from α and constant material parameters ρ ℓ , ρ g , µ ℓ , and µ g . S denotes the strain-rate tensor ( S ≡ [∇u + (∇u) ⊤ ]/2 ), f is the forcing term that prevents the turbulence from decaying, and f σ is the force per unit of volume due to surface tension calculated as f σ = σ κnδ(s) , where σ and κ are the surface tension coefficient and the radius of curvature of the interface between the two fluids (i.e., the droplet and the surrounding fluid), n is the unit vector normal to the interface and directed towards the fluid, with respect to which the interface is concave (droplet), and δ is the Dirac δ-function required to impose f σ only at the interface between the two fluids, and s is a normal coordinate centred at the interface, such that s = 0 at the interface. To take into account the influence of passive vaporization on the concentration fields in the gas phase, an additional equation describing the evolution of a representative scalar ̟ of a passive evaporation process with a low temperature level at the interface is set to be equal to the normalized vapour mass fraction, i.e., ̟ = ϒ v /ϒ vs , where ϒ v denotes the vapour mass fraction and ϒ vs is the saturation value of that fraction based on the Clausius-Clapeyron equation. The scalar ̟ is set to one at the interface following the methodology described in Duret et al. 16 ; it then evolves in the gas phase by means of convection and diffusion as follows: where diffusion coefficient D is set to D = µ g /ρ g , corresponding to a unity Schmidt number. At the initialization of evaporation in the simulation, ̟ is set to zero in the gas and unity at the interface and liquid. The scalar equation is solved for the entire domain. In the liquid, ̟ is fixed at unity during the simulation.
Numerical methods. The above governing equations (14) are solved based on the kernel of the two-component incompressible flows solver Aphrós. For a detailed description of the numerical algorithm used, readers are referred to Karnakov et al. 17 . The equations are discretized on a uniform Cartesian mesh using a finite volume method based on Chorin's projection for pressure coupling 18 and the Bell-Colella-Glaz scheme 19 . The Poisson equation is iteratively solved using the Hypre library 20 . To transport the interface, the solver uses the volumeof-fluid method VOF/PLIC with piecewise linear reconstruction 21 . In the VOF/PLIC method, the normals are estimated using the mixed Youngs-centred scheme, which is a hybridization of Youngs' scheme and the height functions, combining the advantages of both methodologies. To compute the distance from the interface, the marching cubes algorithm 22 is used to estimate the signed distance function , which measures the shortest distance to the liquid-gas interface, such that = 0 at the iterface, � > 0 indicates liquid region and � < 0 indicates gas region.
To sustain the turbulent kinetic energy of the velocity field at a desired level during the numerical simulation, the control-based linear forcing approach of Bassenne et al. 23 , which was initially developed for the gas phase, is here extended here for the liquid-gas mixtures. This approach is used to achieve statistically steady target values of kinetic energy and energy dissipation rate after several large-eddy turn over times. The method of Bassenne et al. 23 consists of introducing low-pass filtering to the flow as follows: where K = �u i u i �/2 and ε = � µ ρ ∂u i ∂x j u i ∂x j � are turbulent kinetic energy and the energy dissipation rate, respectively. K ∞ and ε ∞ denote their corresponding steady mean values, τ ∞ is the integral time based on K ∞ and ε ∞ as Numerical set-up. We initialized the flow with sufficiently enough liquid structures (more precisely, eight droplets were confined into a periodic domain whose sum is equal to the prescribed liquid volume fraction as in Duret et al. 16 ). Subsequently. we applied the forcing procedure to reach a given level of turbulence intensity. Upon obtaining a statistically-converged field, the passive evaporation process was activated. This strategy ensures that the flow almost reaches equilibrium between kinetic energy and surface tension/viscous forces without generating either i) violent small-scale atomization that is difficult to treat numerically or ii) unphysical shear stress across the surface. The shortcomings of this approach are that i) it is highly time consuming, and ii) it is constrained by the moderate values of Reynolds numbers lower than those achievable for single phase flows. The setup under investigation is characterised by several key dimensionless parameters. Herein, all these parameters were considered when the gas reached a 40% saturated state. The physical parameters used are similar to those in previous work 8,24 . The main difference between previous approaches and that used in the present study is the value of the forcing kinetic energy (doubled in this work to K ∞ = 7.2 m 2 /s 2 ). Moreover, the forcing energy dissipation rate is set to ε ∞ = 0.52 × 10 6 m 2 /s 3 . The computational domain has a size L 3 , with a uniform Cartesian grid featuring 512 3 points. This resolution is also twice as large as that used in previous works, since the focus of the present analysis is to analyse small-scale features. The initial density and dynamic viscosity ratios are set to 30, leading to the same kinematic viscosity in both the liquid and gas phases. Surface density σ is set to 0.0135 N.m −1 and the box is filled with φ ℓ = 10% . For both the liquid and gas phases, the Weber numbers can be defined in several ways. In this study, similarly to Dodd & Ferrante 25 , the numbers are defined based on velocity fluctuations, thus We = ρK ∞ L/σ . The Reynolds number is defined as Re = ρ √ K ∞ L/µ and the liquid Ohnesorge number as Oh ℓ = √ We ℓ /Re ℓ . Another relevant non-dimensional parameter, which classically characterizes single-phase homogeneous isotropic turbulence, is the Taylor microscale Reynolds number Re , which is defined as follows Re = 2K ∞ 5ρ g /3µ g ε 0 , where = 10µ g K ∞ /ρ g ε ∞ is the Taylor microscale. This simple configuration is representative of liquid-gas flows that may be encountered in application such injection and atomization 8 . The Kolmogorov length scale is defined as η g = ν 3 g /ε ∞ 1/4 . The resolution parameter �x/η g is 1.78, where x is the uniform grid spacing in each coordinate direction. This estimation ensures that small scale features are well resolved in the present simulation 26 . The corresponding set of fluid parameters is summarized in Table 1.
Assessment of mesh convergence. Unlike the single phase flow, there is no theoretical framework to define the smallest length scale in turbulent liquid-gas flows. Consequently, to ensure that the turbulence is deemed developed and realistic, only a mesh convergence study can be performed. In this study, the assessment of grid convergence is evaluated on the basis of four metrics which are based, respectively, on the: (i) total surface area per unit volume , (ii) temporal evolution of the normalized vapour mass fraction, (iii) velocity powerspectra and (iv) norm of the strain rate tensor. The first metric, is related to the mean liquid-gas surface density that determines the total amount of surface area divided by the simulation volume, denoted as 24 (see Table 2). We noted that this metric quickly converged at higher resolutions, by 1% less when the mesh resolution of 256 3 was modified to 512 3 . The second metric can describe the mixing and evaporation processes, which are estimated from the temporal evolution of the normalized vapour mass fraction ̟ . The temporal evolution of the mean scalar value ̟ and maximum value of the scalar RMS of ̟ up to the time required for the gas to reach 40% of its saturation state is reported on Fig. 2. The high-resolution case exhibits a similar evolution of ̟ and max (̟ RMS ) as in the lower-resolutions cases. The cases 256 3 and 512 3 exhibit small discrepancies, indicating, therefore, that the small-scale effects are resolved on the higher resolution mesh. A grid convergence based on the velocity spectra were performed, as shown in Fig. 3. A large range of the spectra remains unchanged between the three grid resolutions. However, the power spectra of 256 3 resolution and 512 3 resolution overlap in almost all scale ranges, which indicates the convergence of velocity power spectra under this present grid refinement. The last criteria is the norm of the strain rate tensor S A ij S A ij , which can be considered as a good candidate to assess the small-scale features of the velocity field. This metrics converges to a plateau when the resolution reaches 256 3 (see Table 3). Therefore, a resolution of 512 3 could be considered as being fully appropriate for the present study.   www.nature.com/scientificreports/ The invariants of the velocity-gradient tensor and the scalar ̟ in the carrier fluid were analysed as the droplet interface approached. The second-order central difference is used to estimate all the gradients at every point in the flow, except near the liquid-gas interface, where the central difference stencil would lead to mixing between the liquid and gas velocities. In that scenario, the second-order one-sided scheme approximation was used. To discriminate the possible influence of the liquid-gas interface on both the velocity and scalar fields topologies, three layers were considered based on their distance from the interface as follows: −2.8 < �/η g < 0 (near the droplets surface), −7 < �/η g < −2.8 (intermediate region between the carrier-fluid and the droplet) and �/η g < −7 (far from the droplets). These regions were labelled as + c (c for close), + i (i for intermediate) and + f (f for far), respectively. We used the Kolmogorov length scale as the viscous scale reference quantity for normalizing distances from the interface. Dodd & Jofre 12 used a viscous length scale based on the velocity and length scales near the droplet surface in their decaying isotropic turbulence laden with droplets of δ g = ν g ρ g /τ g , where τ g = µ g 4ε 0 /15ν g . The relationship between η g and δ g is δ g = (15/4) 1/4 η g ≈ 1.39η g . Therefore, the near surface to the droplet is twice the viscous length scale. The occupancy of near the droplets surface grid points (expressed as a percentage of realizations) is almost 6% , while the intermediate and far regions occupy 10% and 74%, respectively.
Typical instantaneous snapshot of the liquid-gas turbulence flowfield is shown in Fig. 4 for the time at which the gas reaches 40% of the saturation state. Despite having a relatively small Weber number, the carrier-liquid interface undergoes complex turbulent breakup and appears highly wrinkled and tortuous, comprising elongated ligaments, drops, bag-like structures and bubbles, although the liquid was laminar at the initialization of the simulation.

Results and discussion
Invariants of the velocity gradient, rate-of-strain, and rate-of-rotation tensors considering distance from the interface. In this section, we study the dynamics, topology, and geometry of the gaseous fluid using the invariants of the velocity gradient, rate-of-strain, and rate-of-rotation tensors. All joint PDFs presented in this work cover 95% of events observed in the flow at each distance from the liquid-gas surface. Firstly, we topologically classified the local flow regions using the standardized difference between S A and A as proposed by 27 and κ S A , A = 1 correspond to pure rotational flow, pure shear flow and pure elongational flow, respectively.
, the region is dominated by eddies, and if −1/3 < κ S A , A < 1/3 , the region is dominated by shear. As shown in Fig. 5a, the flow becomes increasingly shear-dominated near the carrier-liquid interface. This is quantified by plotting the histograms of the flow topology parameter κ S A , A for different regions, as shown in Fig. 5b. All calculated curves exhibit a dominant peak in correspondence at κ S A , A = 0 , suggesting a dominantly shear flow. Close to the carrierliquid interface, the right tail drops to zero for κ S A , A = 1 more quickly than for the other two regions, on the other hand, when far from the interface, the PDFs of κ S A , A indicate the increased presence of purely elongational flow.
(Q A , R A ) diagrams. The joint PDFs of the (R A , Q A ) map exhibit a characteristic teardrop shape at + f with narrow cusp in R A = 0 and A = 0 (see Fig. 6c), in which Q A and R A are particularity well-correlated  www.nature.com/scientificreports/ in the following regions: i) the upper left quadrant SFS ( Q A > 0 and R A < 0 ) indicating a predominance of biaxial stretching of the fluid elements and ii) the lower right quadrant UNSS ( Q A < 0 and R A > 0 ) associated with a predominance of enstrophy production by vortex stretching. The observed teardrop shape has previously been observed in single phase HIT, boundary layers, mixing layers and channel flows 28 . The modification of turbulence outside the liquid boundary layer is undetected when looking only at fine scale-scale motions. The isocontours lines corresponding to a constant probability density tend to broaden and increase monotonically as the carrier-liquid interface is approached (see Fig. 6a and b), indicating an increase in magnitude of velocity gradient due to the liquid-gas interface inertia causing deviations in gaseous velocity and, thus, increased gaseous velocity gradient. Note also that i) an increase in the correlation between Q A and R A in the UFC and SNSS was observed when approaching the liquid-gas interface and ii) there is an almost equal probability of enstrophy-dominated regions (UFC and SFS), leading to symmetric (R A , Q A ) map, similar to the single phase high-expansion regions noted in forced compressible HIT 28 and shocked turbulence 29 . Recently, Hasslberger et al. 10 proposed a model predicting that the total percentage of volume fraction is expected to be roughly 80% for focal topologies (40% for UFC and 40% for SFS) and 20% for nodal topologies (10% for SNSS and 10% for UNSS). The configuration investigated in the present study deviates slightly from the idealized model but can be considered approximatively valid for comparison purposes. We also noticed that, in terms of the probability of occurrence in each station, the region + i exhibits a transitional shape with a more pronounced tail than at + and larger velocity gradients than at + f . Similar symmetric behaviour between R A and Q A was reported by Hasslberger et al. 11 for bubble induced turbulence. Moreover, the analytical solution around a bubble obeying Stokes flow yields qualitatively similar behaviour.
(Q S A , R S A ) diagrams. The joint PDFs of (Q S A , R S A ) are depicted in Fig. 7 in order to analyse the geometry of straining of the fluid element in every station. Equation (2) indicates that locations with Q S A ≤ 0 are associated with high level of dissipation. Moreover, equation (3) implies that R S A has the same sign as 2 . Thus, R S A > 0 (resp. R A < 0 ) implies the production (resp. destruction) of S A,2 , and is thus associated with sheetlike structure (resp. tubelike structures). The most striking feature of the joint PDFs of (Q S A , R S A ) is the decreasing preference for a rate-of-strain topology of the type saddle-saddle-unstable-node (two positive eigenvalues and one negative, i.e., expansion of fluid elements) when approaching the liquid-gas interface (see Fig. 7b, c). The  www.nature.com/scientificreports/ probability of finding R S A > 0 becomes slightly greater than finding R S A < 0 in the vicinity of the carrierliquid interface and, thus, the joint PDF of (R S A , Q S A ) becomes slightly tilted towards R S A > 0 and tends to be symmetric along R S A = 0 as seen in Fig. 7a; this implies the presence of more contraction in the region. It is also notable that the lines of constant probability increases when approaching the carrier-liquid interface, resulting in a higher flow dissipation rate close to the interface. Mapping of the joint PDFs of Q S A and R S A attains, at least qualitatively, the characteristic shape, for example, in single-phase mixing layers 30 , channel flow 15 , and homogeneous isotropic turbulence 31 . This diagram highlights that the information gleaned from the (Q A , R A ) diagram is insufficient to characterize turbulent flows since the similar 'teardrop structures' masks important turbulence information. Fig. 8 indicate that the bulk of the data tend to be concentrated near the origin, whereas high gradient motions are distributed differently, i.e., as functions of distance from the carrier-liquid interface. In particular, in the vicinity of the interface, most data fall on the 45 • line through the origin, representing a region in the flow with high dissipation accompanied by high enstrophy density (see Fig. 8a). This is consistent with Fig. 7a, which shows a local vortex sheet, i.e., the rate of strain is dominated by the velocity gradient within the sheet. In the closest region from the carrier-liquid interface, the joint PDF of (Q � A , −Q S A ) shows a marked tendency to be aligned with the vertical line defined by Q A = 0 . This attests to a predominance of dissipation (strain production) over enstrophy, i.e., irrotational dissipation. Far from the liquid-gas interface, the peak value of Q A reaches larger values than the peak value of Q S A within the same probability contour (this was also observed in single phase HIT 31 ), whereas the opposite occurs close to the interface (a similar behavior was observed in the log-law and wake regions of channel flow simulations by Blackburn et al. 15 ).
To more accurately analyse the contribution each topology of (Q � A , −Q S A ) maps on the gaseous flowfield, Kevlahan et al. 32 introduced three different flow regions based on enstrophy and dissipation: (i) the convergence region, Q � A < −Q S A /2 , where irrotational strain is high compared with the rotational strain, (ii) the shear region, −Q S A /2 ≤ Q � A ≤ −2Q S A , where irrotational and rotational strain are approximately equal, and (iii) eddy region, Q � A > −2Q S A , where rotational straining is high compared with irrotational straining. Far from the carrier-liquid interface, 20% of the flow volume is in the eddy zone, 45% in the shear zone and 35% in the convergence zone. These values are quite similar to those in the joint PDFs obtained for single-phase incompressible HIT by Ooi et al. 31 (the percentage of eddy, shear and convergence zones in that study were 19.2%, 44.6% and 36.2%, respectively). In the vicinity of the interface, both eddy and convergence regions have equal probabilities, while most data fall within the shear region.
Normal and non-normal effects on the velocity gradient dynamics. We first estimated the relative importance of the normal and non-normal contributions for the decomposition given by A = B + C . For �C� ≈ 0 , it is legitimate to model the behaviour of A using its eigenvalues. For the flow topology parameter κ S A , A , we plot in Fig. 9 the standardized difference κ B,C = (�B� − �C�)/(�B� + �C�) to estimate the contributions of normal and non-normal effects to the dynamic of the velocity gradient. Normal effects are considered to be important when values of κ B,C are close to 1, whereas they become negligible as κ B,C approaches −1 . In the gaseous flow field far from the carrier-liquid interface, the overall mode of the distribution of κ B,C was slightly negative, with a median close to zero, indicating that the non-normal effects in the dynamics of two-phase homogeneous isotropic turbulent flow are significant to A relative to the eigenvalues. Similar tendencies have previously been noted for single phase HIT 13 . However, there is a tendency for more strongly negative values of κ B,C close to the carrier-liquid interface, suggesting that tensor non-normality plays a greater role in these dynamics when moving into the boundary layer close to the carrier-liquid interface.
When conditioning the overall PDF of κ B,C by four regions of joint PDF of (R A , Q A ) , as shown in Fig. 10, we observed that i) the contribution of C becomes higher where R A < 0 regardless of the distance from the interface, ii) the departure from the zero mode was reduced far from the interface, and iii) the left skewness   www.nature.com/scientificreports/ �S B � > �S C � + � C � , or S B > 2 S C . Therefore, adopting an eigenvalue-based description of the flow is essentially effective in the vicinity of the Vieillefosse tail in the lower hand quadrant of the (Q A , R A ) plane, which acts as an attractor for the dynamics of the restricted Euler (eigenvalue-based) dynamics of A (which amounts to neglecting the anisotropic portion of the pressure Hessian tensor) 13,33 . To evaluate the effect of non-normality on the production of enstrophy and dissipation presented in (12), we again used the standardized differences Fig. 11). Note that � B � = 0 throughout the SNSS and UNSS regions ( � A < 0 ), and � C � = �S B � . The dynamics of the rate-of-rotation and rate-of-strain tensors are dictated only by the non-normal part, C , particularly when close to the carrier-liquid interface, as inferred by the strong negative skewness of both κ S B ,S C and κ B ,S C . This suggests that the properties of vortical structures reported in the literature, such as topology and alignment, are poorly described by the eigenvalues and depend on non-normal contributions.
When the results of κ S B ,S C are partitioned by the four regions of (R A , Q A ) space (see Fig. 12), we found that when R A < 0 , the distribution of S C dominates over that of R A > 0 (consistent with the modes of κ S B ,S C ). As anticipated from the analysis of Fig. 11b, we observed a difference between the SNSS and UNSS regions in terms of their modal values of κ S B ,S C , which were 0.07 and 0.19, respectively. Close to the carrier-liquid interface, the non-normal effects are predominant throughout the entirety of (R A , Q A ) space, with the least prominent effects found in regions characterized by � A < 0 . Overall, when enstrophy production dominates, there is a greater probability for larger shear effects than the eigenvalue-associated effects related to S A and A .
Factors affecting enstrophy production and strain self-amplification. In this section, we analyse the factors that influence important velocity gradient terms, in particular, enstrophy production and strain self amplification. For incompressible flows, the alignment between vorticity ω A and the eigenvectors of the strainrate tensor S A is presumably the primary mechanism controlling the transfer of turbulent kinetic energy from large scales to small ones through vortex stretching 34 . Analyses of the ω A − S A coupling are therefore required to understand the nature of fluctuations in strain-rate in turbulent flow. The evolution of vorticity and strain fields are strongly connected in a non-linear manner. In particular, the dynamics of both fields are primarily driven by self-amplification demonstrated by the transport equations for enstrophy and strain-rate norm squared as follows 35 : ij is the viscous stress tensor and T cap is an extra potential enstrophy production term, which occurs due to the the capillary force itself 11 , and Note that equations (18) and (19) are valid only for single phase flows without interfaces. The production of enstrophy and strain are both dominated by the terms ω A i ω A j S A ij and S A ij S A jk S A ki , which are expressed using the eigenvectors and eigenvalues of S A as follows:  responding to the alignment between the vorticity vector for ω A and the i th eigenvector for the strain rate tensor of A , such that ˚ A i and e A i are the normalized eigenvalues and eigenvectors of S A , respectively. The eigenvectors are scaled to unit length, whereas the eigenvalues are normalized by the magnitude of strain rate ˚ A i = A i /�S A � . The eigenvalues are purely real (because S A is symmetric) and are in decreasing order: Furthermore, the relation ∇ · u = A i holds, in the case of incompressible flows this leads to A 1 > 0 and A 3 < 0 . Incompressibility also dictates that the maximum value obtained by any one normalised strain rate is skews towards positive values. The tendency observed at the interface signifies the prevalence of vortex stretching over vortex compression, which is considered to be a universal characteristic of small-scale features in turbulent flows 2 . A similar shape was noted in the inner layer of compressible turbulent boundary layers 36 .
To understand the coupling between enstrophy production and strain self-amplification, we used the joint Fig. 14). Far from the interface, the contour lines suggests that the carrier gas have a propensity to be located in the top-left quadrant (i.e., , demonstrating that strong production of enstrophy inhibits strong production of strain and vice versa. Moreover, for Near the interface, the shape of the joint PDF shows more strongly compressed and squeezed contour lines indicating that enstrophy production coincides with dissipation production, while enstrophy destruction coincides with dissipation destruction. This behaviour is related to the fact that strain rate and vorticity are simultaneously amplified along the vortex sheets.

Alignment of vorticity vector and strain-rate eigenvectors. The alignment of ω A with the principal strain rate axes directly influences the sign of ω
indicate the production of enstrophy due to the local alignment of vorticity with extensive strain rate eigenvectors, whereas negative values of ω A i ω A j S A ij correspond to the attenuation of enstrophy by vortex compression. Numerous studies, using both experiments and DNS, have revealed that the vorticity vector preferentially aligns with the eigenvector Figure 13.  www.nature.com/scientificreports/ corresponding to the intermediate eigenvalue of the strain tensor, this is considered counter-intuitive because the vorticity should preferentially aligns with largest eigenvector. This can be understood by the conservation of angular momentum 37 . It has also been highlighted that the intermediate eigenvalue is predominantly positive 5 The geometrical analysis far from the carrier-liquid interface reported in Fig. 15 clearly shows that vorticity ω A i) preferentially aligns with the intermediate strain direction e A 2 , ii) exhibits an arbitrary alignment to the extensive direction e A 1 , and iii) tends to be normal to the most compressive principal direction e A 3 . The alignment of vorticity with the principal strain direction agrees with the findings of both single phase HIT 38 , and weak shear layers 5 . However, this alignment behaviour changes in premixed flames, especially for small Lewis numbers and high Damköhler numbers 39 . Close to the carrier-liquid interface, we found that i) the parallel alignment of ω A with the intermediate strain direction e A 2 becomes more likely, ii) a higher frequency of orthogonality is observed between ω A and the compressive strain direction e A 3 , and iii) no special orientation of ω A and the extensive strain direction does not hold, as evidenced by a strong tendency towards orthogonal alignment between ω A and e A 1 with the same magnitude as that obtained between ω A and e A 3 . A similar change of alignment was observed when approaching the wall in a compressible turbulent boundary layer flow along a flat plate in the absence of either a streamwise mean pressure gradient 40 , or downstream of normal shock wave 29 . This behavior is also obtained under certain conditions in premixed turbulent combustion within the reaction zone 39 . The additive decomposition of A into its normal and non-normal parts provides a deeper exploration of previous analysis of turbulence alignment structure, since we can define two additional vorticity vectors ω B (for B ) and ω B (for C ), and six strain rate eigenvectors e B and e C . A subset of these permutations are shown in Fig. 16. Previous results have indicated that vorticity in homogeneous turbulent flows is preferentially aligned to the eigenvector corresponding to the intermediate eigenvalue of the strain tensor S A . Indeed, it has been assumed 41 that ω A stretches in the direction of any eigenvector corresponding to a the strain tensor with a positive eigenvalue, and will eventually become aligned with the eigenvector of the most positive A 1 . The latter dominant alignment is valid when only the normal contribution of A is taken into account, as shown in Fig. 16a. However, a stronger alignment is found between A and A 2 when considering only the non-normal contribution of A , with a noticeable anti-alignment between A 1 and A 3 (see Fig. 16c). This is consistent with the finding of Keylock 13 for single phase HIT. Multiscale analysis by Leung 34 , Doan et al. 42 and Ahmed et al. 43 also indicate that We noted also that i) the alignment between ω A and A 2 is less strong in the vicinity of the liquid-gas interface, whereas ii) the orthogonality between ω A and A 3 and A 1 is much stronger in this region. The results shown in Fig. 16d and b, i.e., the cosine angle between ω B and B i and between ω C and C i , are similar to those shown in Fig. 16a and c.

Contribution of A
i to enstrophy production and strain self-amplification.. We investigated the statistics of the strain-rate eigenvalues in equations (20) and (21) in order to shed more light on the mechanisms of enstrophy production and strain self-amplification. Table 4 shows the values of which represents the positive source term in strain product equation (19). Far from the carrierliquid interface, there are two extensive strain-rate eigenvalues, which are comparable to both experimental and numerical HIT. This tendency indicates that the positive sign of 3 . As a consequence, the strain self-amplification process is associated with the stretching and compressing velocity gradients into sheets. When approaching the interface, the values of The PDFs of the strain-rate eigenvalues as shown in Fig. 17a. In the vicinity of the carrier-liquid interface, the PDFs of ˚ A 1 and ˚ A 3 are almost mirror images of the other, while the PDF of the normalized intermediate eigenvalue is symmetric around zero, consistent with increasing two-dimensionality of the flow close to the interface. The PDF of the intermediate strain rate is clearly skewed towards positive values far from the interface, indicating that its behaviour is predominantly driven by vortex stretching and, thus, the amplification of enstrophy. This observation is a characteristic feature of small scale turbulence 38 . To satisfy the incompressibility constraint of the flow, the PDF of A 3 exhibits a negatively skewed shape, which breaks its mirror image symmetry with A 1 . The statistics of A 2 were further investigated by plotting the dimensionless quantity β A = √ 6˚ A 2 in Fig. 17b. Following this normalization proposed by Ashurst et al. 38 , β A should lie within the range [−1, 1] owing to the incompressibility constraint, which is evidenced in our present study. The PDF of β A demonstrates, as expected, that i) the values of β A are mostly positive far from the interface, and ii) its right skewness decreases close the interface, becoming more symmetric with respect to zero.
Disaggregating the behaviour of β A in the (Q A , R A ) map shows that, with increasing distance from the interface, β A is strongly positive along the densely populated right discriminant line (UNSS), exhibiting lower  . 18c). Moreover, the dimensionless values β A are likely to be more expansive in the strain-dominated streamlines (UNSS and SNSS ) than in the rotation-dominated streamlines (SFS and UFC ). However, close to the interface (see Fig. 18a and b), the contributions of the nodal streamline region become sufficient to symmetrize the global PDF of β A as shown in Fig. 18c. We also identified a noticeable increase in the occurrence of β A in UFC topology close to the interface.
To obtain a more complete understanding of the distribution of β A in the second and third invariants of A , we introduced the normalized second and third invariants of the normalized velocity gradient tensor b ij = A ij /�A� as suggested by Das & Girimaji 46 as q A = Q A /�A� 2 , r A = R A /�A� 3 . The main difference between the previously presented (Q A , R A ) space of Chong et al. 4 (featuring normalization by the mean value of the magnitude of vorticity) and the (q A , r A ) space is that the latter provides a mathematically bounded phase plane allowing to deduce, among other features, a streamline shape as observed in (Q A , R A ) map. Indeed, it can be shown 46 that −1/2 ≤ q A ≤ 1/2 and − √ 3/9 ≤ r A ≤ √ 3/9 . The lower panel of Fig. 18 shows the conditional normalized mean intermediate strain rate eigenvalue β A in the (q A , r A ) map. We note a close correspondence between the results shown in the upper and lower panels. Moreover, β A seems to show negligible dependence on the second normalized invariant q A , exhibiting instead a monotonic increase with respect to the third invariant q A far from the interface (see Fig. 18f). As we approach the interface (see Fig. 18d and e), the probability of being in the first quadrant becomes higher, which explains the symmetrization of the distribution of β A observed in Fig. 18a. Statistics of the passive scalar gradient. The scalar dissipation rate (SDR) ≡ N ̟ := D∇̟ · ∇̟ plays a crucial role in determining the rate of micro-mixing and, thus, in the modelling of turbulent reacting flows. The scalar dissipation rate is also used for premixed combustion modelling (see Chakraborty et al. 47 and references therein). The focus of this section is the turbulence-scalar interaction (TSI) := −2ρD∇̟ · S A · ∇̟ 48 . This quantity, characterizing the coupling between the velocity and scalar fields in the SDR transport equation, is a leading order source term in the SDR evolution 29 . When expressed in the eigenframe e A i of strain rate S A , the TSI term becomes −2ρD∇̟ · S A · ∇̟ = −2ρN ̟ A i cos 2 e A i , n ̟ , where n ̟ = ∇̟/�∇̟ � denotes the normalized scalar gradient. We note that the sign of the TSI will be dictated by the cosines | cos e A i , n ̟ | ≡ ζ e A i ,n ̟ and associated eigenvalues A i . Following our discussion of the contribution of A i , we here discuss the contribution of ζ e A i ,n ̟ . Geometrical alignment in a single-phase HIT was first reported in Ashurst et al. 38 . It was found that the fluctuating scalar-gradient vector preferentially aligned with the compressive strain-rate eigenvector e A 3 , which is perpendicular to the intermediate strain-rate eigenvector e A 2 . In Fig. 19, we show, for the carrier-liquid, changes in alignments between the scalar gradient, n ̟ , and both principal strain directions, by means of space-averaged PDFs of the cosine of the angles between those directions ζ e A i ,n ̟ at different distances from the interface. The preferred orientation of n ̟ is reflected along each of the principal strain eigenvectors. Indeed, far from the interface, as inferred from Fig. 19c, the largest contributions to the spatially mean square value of n ̟ are along e A 3 , while alignments with the other two strain direction are very similar to the results of Ashurst et al. 38 for single phase HIT. As the interface is approached, the magnitudes of the alignment change gradually. The orientation of the scalar gradient, lies mostly at about 45 • to both the e A 1 and e A   49 near the wall of a turbulent channel flow.
To assess the importance of normal and non-normal effects on the observed alignment in Fig. 19, we disaggregated this alignment between the scalar gradient n ̟ and the eigenvectors e B i and e C i corresponding to strain-rate tensors of B and C , respectively, as portrayed in Fig. 20. The most striking observation is that the alignments of ζ e A i ,n ̟ are primarily driven by the non-normal contribution. For the three regions with respect to the liquidgas interface, the normal effects of A modify the magnitude of the alignments, whereas the general alignment ζ e A i ,n ̟ emerges as a direct consequence of the alignment between the scalar gradient and the eigenvectors of the strain rate of non-normal tensor C.
The PDFs of the alignment between n ̟ and e A i in the four topologies classified in the map (Q A , R A ) are shown in Fig. 21. Preferential alignment of n ̟ with e A 3 was noted for the four topologies independent of their distance from the carrier-liquid interface; these where ordered left to right from most to least aligned as UNSS > SNSS > SFS > UFC . This suggests that non-focal strain-dominated topologies enhanced the alignment ζ e A 3 ,n ̟ to a greater extent than the focal topologies. Moreover, the alignment of n ̟ to both e A 1 and e A 3 at about 45 • is much predominant in the focal topologies SFS and UFC . Close to the interface, we note that i) the alignment ζ e A 1 ,n ̟ exhibits a decorrelated trend in non-focal topologies, and ii) n ̟ is statistically perpendicular to e A 2 in the topologies. The different distributions observed in Fig. 21 were confirmed by plotting the conditional mean ζ e A i ,n ̟ in the (q A , r A ) diagram close and far from the liquid-gas interface (see Fig. 22). Specifically, the strong Figure 19. PDFs of alignments between n ̟ , and eigenvectors e A i , corresponding to the eigenvalues A i of the strain-rate tensor, conditioned on different distances from the interface. www.nature.com/scientificreports/ alignment ζ e A 3 ,n ̟ in the UNSS is visible far from the interface, and its probability becomes higher in the focal topology close to the interface. Furthermore, the transverse alignment of n ̟ with e A 2 is mostly concentrated in the negative and positive vicinity of r A ≈ 0.

Conclusions
Direct numerical simulation of a forced turbulent two-phase evaporating flow is performed in order to characterize the invariants of the velocity gradient tensor and local flow topology and thus study the dynamics of turbulent scalar interactions. The topological and dissipating behaviour of the flow were investigated for three regions: in the vicinity of liquid pockets, far from these regions, and in an intermediate region as a means of quantifying the transition between the aforementioned regions. Far from the liquid-gas interface, the gaseous flow is almost unaffected by the presence of the liquid phase. Indeed, the joint PDFs of velocity-gradient, rateof-strain, and rate-of-rotation tensors far from the liquid-gas interface closely resemble canonical and universal single-phase homogeneous isotropic turbulence. However, close to the interface with the liquid-pockets, the flow topology exhibits a boundary-layer-like flow with a predominance of vortex sheets. In the context of LES, these findings suggest viscous regions behaving as turbulent wall-bounded flows, and thus, from the point view of turbulence modelling, a strategy that adapts subgrid modelling to capture droplets interface is plausible. The enstrophy production rate, strain self-amplification and scalar-turbulence interaction were also investigated for topologies classified by the velocity gradient tensor. We found that the mechanisms active in the vicinity of the liquid pockets were different from those acting in the far region; again, we noted a commonality between the present findings and turbulent boundary-layer-like flow. Our results highlight the necessity of considering nonnormal effects on the velocity gradient tensor as a means of explaining alignment properties. As demonstrated herein, the application of eigenvalue-focused approaches to turbulence dynamics is not sufficient to understand www.nature.com/scientificreports/ these mechanisms. It is crucial to consider non-normal effects in conjunction with eigenvalues in order to fully understand the complex dynamics of the velocity gradient tensor. On another note, the statistics presented in the current analysis, are obtained by conditioning on the distance from the two-phase interface. Similar approaches were adopted in several previous analyses 12,50 . However, the present investigation provides more details on the dynamics of the turbulence-scalar interaction in the context of a two-phase flow, through the introduction of the non-normal effect statistics in the aforementioned fashion, which is, to the best of the authors' knowledge, have not been considered before. Direct implications of these statistics are discussed in the terms of enstrophy and strain rate generation/dissipation mechanisms which are relevant in the framework of high-fidelity LES model development for multi-phase flows. In this regard, the present analysis points to the importance of considering non-normal effects in future sub-grid scale models to improve their predictive capabilities. Despite the fact that the current analysis is conducted on a canonical configuration, the conclusions made on the basis of small-scale statistics are likely to be valid irrespectively of the flow configuration nature. Nevertheless, conditional statistics on the distance, with respect to the interface, should be considered in a qualitative sense. Indeed, the transition from the interface effect to a standard singlephase behaviour takes place over a distance that may depend on the mean flow nature.
These present results could be helpful in the construction of reliable SGS models. Typically, the characteristic SGS energy dissipation is an important criterion to assess a model performance. In this respect, the deviatoric part of a new SGS model including the non-local contribution (in the sense of the additional Schur decomposition) could be formulated as where the over bar signifies a spatially filtered quantity, and Favre filtered quantity is denoted by a tilde, C B and C C denote the normal coefficient and the non-normal coefficient, respectively, and is the filter width. The model given by (22) can be downgraded into the standard Smagorinsky model 51 in the regions where the non-local effects are insignificant. Furthermore, it has the potential to offer a better prediction of the SGS energy dissipation in comparison to the standard approach especially in non-normal dominant regions. An example of such regions, is at the vicinity of the liquid-gas interface where non-local effects dominate. Figure 22. Conditional average ζ e A i ,n ̟ |q A , r A at the vicinity of the interface (a-c) and far from the interface (d-f).