Extensive strain along gradient trajectories in the turbulent kinetic energy field

Based on direct numerical simulations of forced turbulence, shear turbulence, decaying turbulence, a turbulent channel flow as well as a Kolmogorov flow with Taylor-based Reynolds numbers Reλ between 69 and 295, the normalized probability density function of the length distribution of dissipation elements, the conditional mean scalar difference ⟨Δk∣l⟩ at the extreme points as well as the scaling of the two-point velocity difference along gradient trajectories ⟨Δun⟩ are studied. Using the field of the instantaneous turbulent kinetic energy k as a scalar, we find good agreement between the model equation for as proposed by Wang and Peters (2008 J. Fluid Mech. 608 113–38) and the results obtained in the different direct numerical simulation cases. This confirms the independence of the model solution from both the Reynolds number and the type of turbulent flow, so that it can be considered universally valid. In addition, we show a 2/3 scaling for the mean conditional scalar difference. In the second part of the paper, we examine the scaling of the conditional two-point velocity difference along gradient trajectories. In particular, we compare the linear s/τ scaling, where τ denotes an integral time scale and s the separation arclength along a gradient trajectory in the inertial range as derived by Wang (2009 Phys. Rev. E 79 046325) with the s·a∞ scaling, where a∞ denotes the asymptotic value of the conditional mean strain rate of large dissipation elements.


Introduction
One of the many approaches in turbulence research is to study geometrical structures in terms of critical points in the flow field. Gibson [1] analyzed the behaviour of zero gradient points and minimal gradient surfaces in turbulent scalar fields. He argued that these points are of importance to the problem of turbulent mixing. Wray and Hunt [2] subdivided the flow field into four types of space-filling regions, characterizing them by the second invariant of the velocity derivative tensor Q and the pressure p. More examples of studies of geometrical structures in turbulent fields can be found in [3]- [5]. Dávila and Vassilicos [6] used the spatial distribution of stagnation points in homogenous isotropic turbulence to show that it has a fractal structure of dimension D s = 2 and found that the Richardson constant is an increasing function of the number density of the stagnation points. Goto and Vassilicos [7] relate the energy dissipation rate coefficient to the stagnation point structure of homogeneous isotropic turbulence in order to prove its non-universality.
Based on the extreme points of turbulent scalar fields, i.e. points of vanishing scalar gradient, Wang and Peters [8,9] developed the theory of dissipation elements, which arise as natural geometries in turbulent scalar fields when these are analyzed by means of gradient trajectories. Starting from every grid point, trajectories along the ascending and descending gradient directions can be calculated, which inevitably end at extreme points. All points that share the same two ending points define a finite volume, which is called a dissipation element. These elements are parameterized by two values, namely the linear length l between the extreme points and the scalar difference k at the extreme points, where k is the scalar to be used in the present paper.
Based on this theory, space-filling elements are identified, which allow the reconstruction of statistical properties of the field as a whole in terms of conditional statistics within the elements. Examples of such an analysis can be found in [9,10], where in addition the model constants of the ε-equation in the widely applied k-ε turbulence model are determined via dissipation element analysis in the ε-field. However, taking instead the mixture fraction Z as the underlying scalar field allows the physical interpretation of dissipation elements in the context of the flamelet approach in non-premixed combustion; see [11] for details.
From the definition it follows that the temporal evolution of dissipation elements in turbulent fields is inherently connected with the evolution of their ending points, which are separated by a mean linear distance l m that is of the order of the Taylor microscale λ [8]. While strain and diffusion lead to continuous distortion of an element as a whole, the creation or annihilation of extreme points leads to their abrupt formation or disappearance. These different 3 effects have successfully been considered in a modelled evolution equation for the probability density function (pdf) for the length distribution P(l) of dissipation elements, which in its normalized formP(l) is assumed to be independent of both the Reynolds number and the type of turbulent flow. In this equation, two parameters appear, cf [9], which can be identified as a splitting and on attachment frequency of the elements and thus indirectly correspond to the lifetime of extreme points; see [12,13] for further details. A third important parameter employed in [9] for the description of the pdf of the length distribution is the normalized rate of strainã. The strain rate is defined as the difference of the velocity at the ending points, projected in the direction of the connecting line between the two extreme points. It will be shown below that there is a linear scaling of the mean absolute value of the velocity difference with the curvilinear distance along gradient trajectories for large elements.
The basis for this scaling is a paper by Wang [14], who studied the two-point correlation of scalar gradients along gradient trajectories. Starting from the governing equation for the passive scalar, a decorrelation assumption for the product of two-point scalar gradients and the velocity difference in the triple correlation term is used to derive a positive linear scaling of the velocity difference with curvilinear separation distance along the trajectory s and the inverse of the integral time scale τ . Furthermore, Wang shows that this result is in good agreement with data obtained from direct numerical simulations (DNSs) of homogeneous shear turbulence. It is argued that due to a conditioning of the statistics on gradient trajectories, regions of large extensive strain, which smooth the scalar field, are preferentially extracted, thereby allowing a gradient trajectory to extend over large distances. For relatively short distances in the viscous range however, with s < λ, no analytic expression for the scaling of the velocity difference has been derived so far. Based on an argument involving the curvature term, which arises in the derivation of the governing equation along trajectories, a negative value for the velocity difference is qualitatively concluded in [14].
The intention of the present paper is twofold. In a first step, we will examine the validity of the model for the pdf of the linear element length as derived in [9] as well as the Kolmogorov scaling of the conditional scalar difference for different types of turbulent flows. In a second step, the universality of the two-point velocity difference u n along gradient trajectories is analysed. To this end, we examine the validity of Wang's scaling for long trajectories.
In contrast to previous work, all gradient trajectory and dissipation element analyses presented in this paper are based on the field of the instantaneous turbulent kinetic energy k. This appears reasonable as it allows a direct interaction between the underlying velocity field and the examined scalar field, as well as an evaluation of the energy budget along gradient trajectories. In addition, the field of the turbulent kinetic energy can easily be evaluated from investigations of the three-dimensional (3D) velocity field and consequently simplifies an experimental verification of the theoretical results; see for instance [15].
As numerical test cases, DNSs of various turbulent flows have been performed. Details of the different flow setups are given in section 2. In section 3, we present the results of the dissipation element analysis. The two-point velocity difference along gradient trajectories is discussed in section 4, and the conclusion drawn is presented in section 5.

Direct numerical simulations
DNSs of six different types of turbulent flows, namely homogeneous shear turbulence (cases 1 and 2), homogeneous isotropic forced turbulence (case 3), homogeneous isotropic decaying 4 turbulence (case 4), a channel flow (case 5) and a Kolmogorov flow (case 6), were performed. Cases 1-4 and case 6 were run on the JUGENE Supercomputer of the Research Center Jülich using up to 16 384 CPUs for the calculations and case 5 was run on the HHLR Darmstadt. In all cases the Kolmogorov length η = (ν 3 /ε) 1/4 was resolved, which is particularly important for dissipation element analysis. Taylor-based Reynolds numbers Re λ = u λ/ν between 69 and 295 have been investigated, with the Taylor microscale λ = √ 10νk/ε and the turbulence intensity u = √ 2k/3. For the DNS of cases 1-4 and 6, the incompressible Navier-Stokes equations were solved in a cubic box of size 2π with periodic boundary conditions employing pseudo-spectral methods. The DNS code was parallelized to run on supercomputers using a highly efficient MPI parallelized 2D decompositioning and a 3D-FFT library [16] as the kernel of the code. This enables high-accuracy fast-Fourier transformations used for spatial discretization. Aliasing errors were removed by isotropic truncation applying the 2/3 rule.
In detail, homogenous shear turbulence has been performed for cases 1 and 2 with a mean velocity gradient of S = 1.5 on a grid with 2048 3 and 1024 3 grid points. The temporal advancement is performed by a third-order Runge-Kutta method. The convective term of the Navier-Stokes equations is formulated in the skew-symmetric form in order to reduce aliasing errors and to improve numerical stability and accuracy, cf [17]. Periodic boundary conditions, which are required for the application of spectral methods, cannot be satisfied when a mean gradient is present. To overcome this problem, a coordinate transformation of all dependent variables to a moving frame attached to the mean flow is performed [18]. Since the computational frame gets distorted with advancing time a remeshing procedure is applied to keep the distortion in an appropriate range. The total amount of cpu time needed for one integral time step was 3795 626 cpu hours for case 1 and 216 177 cpu hours for case 2.
DNSs of a homogeneous isotropic forced turbulence have been performed for case 3 on 1024 3 grid points. A pseudo-spectral method developed and implemented by Ruetsch and Maxey [19] is applied. The nonlinear terms are solved using an explicit Adams-Bashforth method, while the linear terms are solved by an implicit Crank-Nicholson method, both of which are of second order. The forcing is implemented using a method developed by Eswaran and Pope [20]. The total amount of cpu time needed for one integral time step was 42 513 cpu hours for case 3.
The algorithm of case 3, but without forcing, was used for the DNS of a homogeneous isotropic decaying turbulence (case 4) on 1024 3 grid points. It employed 327 840 cpu hours for the whole calculation. The initial velocity field is random and isotropic and is generated such that it satisfies a prescribed energy spectrum. The initial energy spectrum is taken from [21] and has the form where The constant κ p is the wave number at which E(κ) has its maximum and is set to κ p = 10. We use σ = 4 in equation (1). After an initial development where the dissipation increases, the decay of the kinetic energy and the dissipation follows a power-law decay k/k 0 = (t/t 0 ) −n and ε/ε 0 = (t/t 0 ) −n−1 , respectively. The decay exponent is found to be n = 1.4 and thus lies well in the range of values reported in the literature; cf [22]- [25]. The Reynolds number at which the flow field is evaluated is Re λ = 71. For case 6 a turbulent Kolmogorov flow was calculated on 1024 3 grid points taking 36 113 cpu hours for one integral time. In contrast to cases 1 and 2, this flow exhibits regions of strong shear as well as regions with weak shear. In order to obtain a statistically steady solution, the value of the first mode of the velocity component u in the x-direction is kept at a constant value in Fourier space, which leads to a sinusoidal mean velocity profile; see figure 1. As this imposed profile introduces a single characteristic length scale to the turbulent flow, it naturally bounds the size of the largest energy containing eddies and thus, in contrast to homogeneous shear turbulence, allows for a statistically steady flow.
In case 5, DNSs of a turbulent channel flow were performed; see table 2. Figure 2 illustrates the geometry and the coordinates of the channel flow. The numerical code uses spectral methods to solve the 3D time-dependent incompressible Navier-Stokes equations and was developed at KTH, Stockholm; see [26,27] for details. The code adopts a spectral method with Fourier decomposition in the streamwise and spanwise directions and discretization with Chebyshev polynomial series in the wall normal direction. Apart from the wall-normal integration the numerical method is similar to the one used in [28]. Time integration is performed using a third-order Runge-Kutta scheme for the advective and forcing terms and second-order Crank-Nicolson for the viscous terms.
We examine a pressure-driven incompressible turbulent flow field in a channel with constant viscosity, with no rotation and no external forces passing through the channel bounded by two parallel walls. The governing equations for such a flow satisfy the Navier-Stokes equations in non-dimensional form with no-slip boundary conditions at the channel walls. In the homogeneous streamwise and spanwise directions periodic boundary conditions are applied and the pressure gradient that drives the flow is fixed.
The friction Reynolds number Re τ = u τ h/ν, where u τ = √ τ/ρ is the friction velocity expressed in terms of the average shear stress τ at the wall and the density ρ, for this DNS is Re τ = 360 and the Taylor-based Reynolds number is Re λ = 69 evaluated in the core region of the channel, respectively.   Table 2. Simulation parameters. L x and L z denote the streamwise and spanwise dimensions of the computational box. y + c is the non-dimension wall-normal resolution in the centre of the channel.
An overview of the numerical parameters and selected mean quantities is presented in table 1. Note that the values for case 5 originate from the channel core region, while for case 6 all y-dependent quantities have been averaged over the inhomogeneous direction.

Dissipation element analysis
The motivation for dissipation elements is the reconstruction of the entire 3D scalar field by means of an adequate description of an element's characteristics. Figure 3 shows for each gradient trajectory of case 2 a scatter plot of the curvilinear length s from all calculated trajectories within an element over the linear length l of the corresponding dissipation element. As expected, the curvilinear length is always larger than the linear length. For small elements  (s 0.25), it can be observed that the curvilinear length follows closely the linear distance, while for larger elements an increasing spread between the curvilinear length of the trajectories of an element and its linear length can be observed. On average s ∼ 1.2l. From this, one can deduce that small elements are quite regular in shape, while larger elements are corrugated and possess a highly complicated geometry, a finding that is in good agreement with previous observations; cf [8,9]. As the linear length is a more easily accessible parameter than an appropriate average of the curvilinear length and, in addition, it is unique for each dissipation element, l and k have been chosen as the statistical parameters for describing dissipation elements.
The corresponding joint probability density function (jpdf) P(l, k) is expected to contain most of the information needed for a statistical reconstruction. Based on a trajectory search algorithm, the field of the turbulent kinetic energy has been analyzed for the different DNS cases and the resulting jpdf for case 1 is shown in figure 4 (to relate the values of l and k given in this figure to other characteristic flow quantities, see table 1). In this illustration, different physical effects are illustrated. In addition to a distinct maximum, one observes a decrease at the origin, corresponding to the annihilation of small elements due to molecular diffusion. The 8 region in the upper right-hand area of the jpdf is dominated by extensive strain, as large elements are exposed to large velocity differences.
Additionally, a dotted line corresponding to the conditional mean k|l is included in figure 4 (left), revealing a scaling with Kolmogorov's 2/3. This finding will be discussed in more detail later; see figure 7. However, the spread of the jpdf around this mean value is obviously not symmetric and can therefore not be described by a Gaussian distribution; see figure 4 (right).
The jpdf can be described by a model equation. Based on the Bayes theorem, it is decomposed into a marginal pdf P(l) of the linear distance and a conditional pdf P( k|l) of the scalar difference, yielding where the marginal pdf P(l) is defined by For this pdf in its normalized formP(l), withP = P l m andl = l/l m , the following model equation was derived in [9]: In this equation,ã represents the conditional mean strain rate a of the elements of length l a = u n |l l (6) normalized by its asymptotic value a ∞ , which is approached for l → ∞.
In equation (6), u n denotes the velocity difference between an element's maximum and minimum projected in the trajectory direction with n = ∇k/ |∇k| .
The results for a(l) have been calculated for all DNS cases and its conditional mean is compared with the model proposed in [9], where the value of a ∞ was also obtained from the DNS. Figure 5 depicts the comparison of equation (9) with DNS data. One observes a qualitatively similar shape for all DNS cases, as well as good agreement between DNS data and the model equation. The value of a ∞ ranges from 0.68 for the Kolmogorov flow to 4.10 for the shear turbulence. Finally denotes the normalized drift velocity due to molecular diffusion in equation (5) and is responsible for the linear decrease ofP(l, t) forl → 0 as will be shown below. The constant c in equation (10) is determined from the condition that the total length of the array must not change, cf [9], and D is the molecular diffusion coefficient. In addition, in equation (5) the two non-dimensionalized numbers s and a appear. These describe the splitting (respectively reconnection) of larger (smaller) elements into smaller (larger) ones and are determined from the normalization and the first moment during the solution of the equation as eigenvalues of the problem, cf [12]. Equation (5) can be solved numerically yielding the steady solution depicted by the solid curve in figure 6. Furthermore, the results for the normalized pdf of the length distribution obtained from the different DNS cases are shown. One observes very good agreement with the model solution. Slight differences can be identified in cases 4 and 6, where the model marginally underpredicts the maximal value of the pdf. The linear increase at the origin as well as the exponential tail however, see especially the log insets, follow more or less closely the predicted solution. For cases 1 and 5, the location and the overall shape of the pdf are tilted slightly to the left, resulting in a small deviation from the model, although the branches left and right of the maximum qualitatively agree nicely. (Note: all data presented for case 5 here and in the following just stem from an evaluation in the core region of the channel.) However, deviations in the exponential tails (see the slopes of the pdf in the insets) of the pdf are hard to interpret due to the limited number of sample points. Nevertheless, figure 6 illustrates that the equation forP(l) seems not to be a function of the Reynolds number as the values of R λ vary from 69 for the case of decaying turbulence to 295 for the shear turbulence. Overall, no conclusive influences due to the Reynolds number or based on the type of turbulent flow can be identified, so that the shape of the non-dimensional marginal pdfP(l) and its model equation may be considered independent of inhomogeneities and anisotropies, a finding that is illustrated in particular by the good agreement between the model and the Kolmogorov flow. For the second term in equation (3), the conditional pdf of k, no conclusive model equation was derived so far. In a first step, we will therefore analyse the conditional mean of k conditioned on the length of the respective dissipation element, defined by In unconditioned statistics the first moment is equal to zero. For statistics based on gradient trajectories however, this is not the case, as the value of the turbulent kinetic energy increases per definition monotonically along a trajectory from the minimum to the maximum point. Consequently, we will study the first order conditional moment to examine its scaling. In [8]- [10] this conditional difference based on various scalar fields, such as for instance a passive scalar φ or the instantaneous dissipation ε, was found to scale with Kolmogorov's 1/3 exponential dependence. In the present case for k however, one expects a scaling with 2/3 based on dimensional grounds. For the different DNS cases presented in chapter 2, the results for k|l are shown in figure 7.
One observes different absolute values for the difference of the instantaneous turbulent kinetic energy between the extreme points of a dissipation element conditioned on its length. (Note: for a better graphical illustration the channel flow data (case 5) have been shifted by four orders of magnitude, the shear flow data (case 1) have been multiplied by 5 and the decaying flow data (case 6) have been multiplied by 5.) The symbols of the different DNS cases included in figure 7, however, show a scaling with 2/3 as indicated by the solid line with a slope of 2/3, which is more or less accurate for the different flow types. The best agreement is obtained in the case of decaying turbulence, while for example in shear turbulence a less pronounced 2/3 scaling is observed. In addition, it is obvious that Kolmogorov's scaling is also valid for dissipation element analysis based on the channel flow and the Kolmogorov flow, which are inhomogeneous and anisotropic flows.

Scaling of the velocity difference along gradient trajectories
In the course of this section, we will study the velocity difference at two points along gradient trajectories, motivated by the findings shown in figure 5. Conditioning on dissipation elements and more specifically on points along one trajectory introduces major differences compared with standard statistics in Cartesian correlation space. Due to the limitation of the two-point statistics on points on the same trajectory, only the correlation of specific points is studied. The conditional statistics along gradient trajectories also introduces differences with respect to the correlation coordinate. In contrast to statistics in correlation space of a Cartesian grid, where the correlation coordinate usually denotes the linear distance between the two points under consideration, the correlation coordinate s is defined as the curvilinear distance of the two points along their trajectory. This obviously also introduces a flow-dependent restriction to the length of the new correlation coordinate as the mean distance of two points is the one between a maximum and a minimum point, which is of the order of the Taylor scale, as discussed in the previous section in the context of the dissipation element length and its pdf.
Compared to statistics in Cartesian coordinates, caution has to be exercised when the concepts of isotropy and homogeneity are used. In homogeneous turbulence, the mean of the fluctuating component of the instantaneous velocity is by definition equal to zero. The velocity along trajectories u n , however, is projected in trajectory direction n, cf (8), and therefore time and space dependent, yielding u n = u · n. The result of this product is a scalar, whose mean value is not by definition equal to zero. This difference also arises when the first conditional moment is studied, which is of particular interest for the present work in trajectory coordinates, while it is equal to zero in the Cartesian system.
In the following, we will study the scaling of the conditioned mean velocity increment u n along gradient trajectories. Based on the governing equation for a passive scalar φ in gradient coordinates, a relation for the inertial range was derived in [14]. Neglecting the viscous term and assuming a decorrelation of the velocity difference from the two-point correlation of the scalar gradient, it is concluded that where the integral time scale is defined as τ = k/ε. As this proportionality has only been validated in [14,29] for homogeneous shear turbulence, we will examine its validity in the following for other flows.
In figure 8, the non-dimensional product u n (τ/λ) is shown as a function of s/λ. In this figure, the Taylor scale is used for the normalization rather than the Kolmogorov or an integral scale, as it is the representative length scale for gradient trajectories. One observes that both the normalized velocity difference and the slope for all DNS cases are negative up to roughly 0.5λ, while the zero-crossing is always close to s/λ = 1. One furthermore finds a distinct quasilinear increase with the separation arclength for all DNS cases beyond the minimum. The slope, however, varies for the different flows. The deviations from the linear scaling for large separation distances s > 3.5λ may be attributed to the small number of gradient trajectories of such a length, an obstacle that is overcome in the DNS of shear turbulence, which employs 2048 3 grid points, as in this case the increase is linear up to s ≈ 8λ. As illustrated in figure 4 for the length distribution of dissipation elements, the majority of elements have a linear length of order λ. The probability of finding an element longer than s > 3.5λ is consequently very small. Even though for the conditioned mean velocity increment the curvilinear arclength, which obviously Non-dimensional first-order velocity structure function u n τ/λ along a gradient trajectory for cases 1-6. is longer than the linear distance, is used, the difference between the linear distance of two extreme points and their curvilinear distance is relatively small, cf figure 3. Summarizing, it can be concluded that while the linear increase of u n with s/λ in the inertial range can be considered as well established, the proportionality constant depends on the flow. This is not surprising since u n is normalized in figure 8 with the integral time scale of the respective DNS case. On the other hand, large dissipation elements are strained by a ∞ such that the projected velocity difference u n at the extreme points of an element in the direction of the linear connecting line is proportional to la ∞ . Therefore, one may expect that the strain rate of individual trajectories within an element, especially for large separation distances s, will be exposed to the extensive strain a ∞ .
To normalize the slope of the profile of u n for large s, we therefore propose to nondimensionalize it by using the strain a ∞ , cf table 3, rather than the integral time scale τ . This is shown in figure 9. Compared with figure 8, the new scaling u n /(a ∞ λ) naturally retains the zero-crossing at approximately s = λ, but in addition collapses the profiles so that the results of the different DNS cases now lie close to each other with a slope of approximately 0.5.
This normalization of u n using the strain rate a ∞ connects the two-point correlation along gradient trajectories with the concept of dissipation elements. Although all gradient trajectories connecting the same two extreme points are only described by the parameters of their dissipation element, this generalized 3D information is still valid when returning to a 1D trajectory.  The scaling of l m may be inferred from equations (5)- (10). Since only the diffusion D and the strain rate a ∞ appear as dimensional parameters in these equations and since ν = D in the simulations, we obtain On the other hand, the Taylor scale λ is proportional to leading to as displayed in figure 10.

Conclusion
Dissipation element analysis based on the instantaneous field of the kinetic energy has been extended from homogenous shear turbulence to various other types of turbulent flows exhibiting a wide range of Taylor-based Reynolds numbers. For all cases the mean length of the dissipation elements was found to be of the order of the Taylor scale. The normalized pdf of the length of the dissipation elements was found to be Reynolds number independent, in agreement with the previous findings for homogeneous shear turbulence. In addition, it was found that the flow configuration does not have a noticeable effect on the shapes of the normalized pdfs. This leads to the conclusion that the relevant physical effects that determine the pdf are independent of residual inhomogeneities or anisotropies. The model equation for the pdf was found to be in good agreement for all the cases discussed. The first-order structure function of the scalar difference at the ending points of the dissipation elements was found to scale approximately as l 2/3 in accordance with dimensional analysis. The velocity difference at the ending points, however, revealed a linear scaling with the element's linear length for all DNS cases. This result was further analyzed in terms of the first-order velocity structure function along gradient trajectories.
In a first step, Wang's scaling was examined for the different turbulent flows, which revealed a negative u n (τ/λ) for small separation distance, followed by a linear increase with a zero-crossing around λ. As shown in [14], the slope of this linear scaling varies from case to case and is therefore not universal. In a second step however, u n was normalized by the asymptotic strain rate calculated for dissipation elements, as strain is the dominating physical mechanism acting on long gradient trajectories. The first order structure function along gradient trajectories normalized with the strain rate u n /(a ∞ λ) consequently exhibited the same shape for all turbulent flows. The successful normalization of statistics along gradient trajectories using the strain rate of dissipation elements, in addition, illustrated the physically meaningful generalization used in dissipation element analysis.