Lagrangian network analysis of turbulent mixing

TL;DR

This paper introduces a novel complex network approach to analyze turbulent mixing from a Lagrangian perspective, effectively capturing flow dynamics and regimes in turbulent channel flow.

Contribution

It presents a new network-based formulation for Lagrangian turbulence analysis, enabling robust and computationally efficient identification of flow regimes and particle interactions.

Findings

Network formalism clearly describes particle swarm dynamics

Identifies transient and long-term flow regimes

Highlights flow features like turbulent mixing and advection

Abstract

A temporal complex network-based approach is proposed as a novel formulation to investigate turbulent mixing from a Lagrangian viewpoint. By exploiting a spatial proximity criterion, the dynamics of a set of fluid particles is geometrized into a time-varying weighted network. Specifically, a numerically solved turbulent channel flow is employed as an exemplifying case. We show that the time-varying network is able to clearly describe the particle swarm dynamics, in a parametrically robust and computationally inexpensive way. The network formalism enables to straightforwardly identify transient and long-term flow regimes, the interplay between turbulent mixing and mean flow advection, and the occurrence of proximity events among particles. Thanks to their versatility and ability to highlight significant flow features, complex networks represent a suitable tool for Lagrangian investigations of turbulence mixing. The present application of complex networks offers a powerful resource for Lagrangian analysis of turbulent flows, thus providing a further step in building bridges between turbulence research and network science.

Equations

$$W_{i,j}(\mathbf{a},\mathbf{\lambda },t^{+})=\sum _{p\in l_i}\sum _{q\in l_j}{I}_{p,q}(\mathbf{a},t^{+})K\left({\mathbf{\lambda }}_{p,q}(t^{+})\right),i,j=1,...,N_y,$$

$$K\equiv K_{p,q}^w=\left(\frac{1}{d_{p,q}+1}-\frac{1}{D_{p,q}+1}\right)/\left(1-\frac{1}{D_{p,q}+1}\right),$$

$$\nabla \cdot \mathbf{u}=0,\frac{\partial \mathbf{u}}{\partial t}+\frac{1}{\rho }\nabla p=\mathbf{f}-\mathbf{\omega }\times \mathbf{u}+\nu {\nabla }^2\mathbf{u},$$

Full text

1 Introduction

One of the most remarkable features of turbulent flows is their ability to strongly enhance transport and mixing processes. This outstanding property plays a crucial role in many natural phenomena and engineering applications, ranging from chemical reactions and combustion mechanisms (e.g., see Warnatz et al., 1996; Nguyen & Papavassiliou, 2018) to biophysics (Seuront & Schmitt, 2004), as well as atmospheric dispersion and geophysical phenomena (e.g., see Pasquill & Smith, 1983; Fernando, 2012). In spite of many efforts, several issues regarding the understanding and modelling of turbulent mixing – e.g., anomalous scaling of statistics and signal intermittency – need to be addressed (Warhaft, 2000; Sawford, 2001; Dimotakis, 2005; Toschi & Bodenschatz, 2009). The Lagrangian viewpoint of turbulent mixing – in which the focus is on particles that are advected by the flow – has turned out to be better suited than the Eulerian view in many cases (Falkovich et al., 2001; Sreenivasan & Schumacher, 2010). To this aim, the time evolution of pair- and multi-particles has usually been explored in terms of geometrical features – such as pairwise mean-square separation or multi-particle shape evolution –, and results are typically ensemble averaged (e.g., see Salazar & Collins, 2009; Pitton et al., 2012; Bianchi et al., 2016; Polanco et al., 2018).

In this work, turbulent mixing is investigated in a Lagrangian way by exploiting the recent advances in complex networks, aiming to extend the level of information of classical statistics. Differently to previous approaches, the geometrical representation into complex networks of the particle dynamics offers a twofold tool. On one hand, the network topology systematically geometries the particle dynamics, thus introducing a schematic and synthetic representation of the particle swarm in time. On the other hand, the network formalism provides a robust and well-established framework for studying non-trivial spatio-temporal dynamics of a discrete set of interacting elements (Boccaletti et al., 2006). Here, we propose a network-based geometrization of particle dynamics in which nodes correspond to fluid particles and links are active by means of a spatial proximity criterion, such that the relative position between particles at any time is enclosed into the network structure. In this way, we obtain a time-varying network which is able to capture in a synthetic way both transient and long-term effects of turbulent mixing, as well as the extent to which particles interact with each other.

Network science has emerged in last two decades as an effective technique to study real-world complex systems (Newman, 2018), and growing attention is given to the application to turbulent flows (e.g., see Charakopoulos et al., 2014; Murugesan & Sujith, 2015; Scarsoglio et al., 2016; Taira et al., 2016; Iacobello et al., 2018). Very recently, network-based analysis of fluid flows both in a Eulerian and in a Lagrangian frame has also been carried out. In general, the focus of Lagrangian approaches has been on fluid transport and coherent structure identification, by exploiting the discrete transfer operator (Ser-Giacomi et al., 2015) and spectral-graph procedures (Hadjighasem et al., 2016; Schlueter-Kuck & Dabiri, 2017; Padberg-Gehle & Schneide, 2017; Schneide et al., 2018). Specifically, the linking criterion is usually based on a similarity measure – e.g., Euclidean distance, correlation coefficient or Granger causality (see Donner et al., 2017, for an overview) – between particle trajectories evaluated in a given time interval; by doing so, the temporal details of particle trajectories do not explicitly emerge. Differently, in this work the spatio-temporal evolution of particles is captured at each time step, providing a rich and detailed time-dependent picture of turbulent mixing. Therefore, the time-varying network formulation intrinsically highlights the temporal development of particle dynamics due to the turbulent motion.

As a paradigm of possible applications, the proposed approach is shown for a numerically solved turbulent channel flow. In this way, we are able to provide physical insight into the interplay between the mean flow advection and wall-normal turbulent mixing on particles at different times, thus highlighting the key role of the spatial inhomogeneity.

2 Flow description and networks building

A DNS of a three-dimensional fully-developed incompressible turbulent channel flow was performed at $\Rey_{\tau}=Hu_{\tau}/\nu=950$, where $H$ is half the channel height, $u_{\tau}$ the friction velocity and $\nu$ the kinematic viscosity (Kuerten & Brouwers, 2013). The length of the domain equals $2\pi H$ in streamwise direction, and $\pi H$ in the spanwise direction; in these two directions, periodic boundary conditions are applied to velocity and fluctuating part of the pressure. A constant mean pressure gradient in the streamwise direction is applied to drive the flow. The grid spacing in wall units equals $7.8$ and $3.9$ in the streamwise and spanwise direction, respectively, while the maximum wall-normal grid spacing (at the channel centre) is $7.8$. The total simulation time is $T^{+}=15200$, while the time-step is $\Delta t^{+}=4.75$, where the $+$ superscript indicates wall-units normalization. The $T^{+}$ value was set greater than the wall-normal (Eulerian) mixing time, $T^{+}_{\epsilon}$, namely the integral time-scale after which the Taylor dispersion analysis is asymptotically valid (Fischer, 1973). By definition, $T^{+}_{\epsilon}\nu/u_{\tau}^{2}=H^{2}/\epsilon_{y}$ (Fernando, 2012), where $\epsilon_{y}=0.067Hu_{\tau}$ is Elder’s vertical mixing coefficient (Elder, 1959); as a result, $T^{+}_{\epsilon}\approx 14180$.

The numerical simulation is run until the statistically stationary condition of the flow is reached. After this condition is fulfilled, particles are seeded in the fully-developed turbulent flow and tracked in time (for more simulation details, see Appendix A). At this initial time, a set of $\left(N_{y}\times N_{z}\right)$ fluid particles was arranged as a uniformly distributed grid in the plane $\left(y^{+},z^{+}\right)$ at $x^{+}=0$, where $(x^{+},y^{+},z^{+})$ are the streamwise, wall-normal and spanwise coordinates, respectively. Although in this work we considered fluid tracers, inertial particles could alternatively be employed. The time-dependent particle positions are then obtained from $\mathrm{d}\boldsymbol{x^{+}}/\mathrm{d}t^{+}=\boldsymbol{u}\left(\boldsymbol{x^{+}}(t^{+}),t^{+}\right)$, where $\boldsymbol{x^{+}}$ is the position of a tracer particle. Since the focus is on the wall-normal mixing process, particles were grouped into $N_{y}$ wall-normal levels, $l_{i}$, where each level includes a spanwise-row of fluid particles at the initial time, i.e. $l_{i}$ depends only on the $y^{+}$ value of particles at $t^{+}=0$ (e.g., figure 1(a) shows five levels highlighted with different colors). In this work, we set $N_{y}=N_{z}=100$, where particles closest to the walls are located at $y^{+}=9.5$ and the remaining ones are separated by $\Delta y^{+}=2\Rey_{\tau}/N_{y}=19$, while the initial spanwise separation is $\Delta z^{+}=29.85$. The initial particle grouping into levels naturally emerges here as the presence of the walls introduces an inhomogeneous direction, $y^{+}$. The proposed initial particle arrangement is suitable to numerically represent, for example, a uniformly distributed release of contaminant in the domain. In general, initial particle grouping can be not easily identifiable, resulting in a non-trivial task which requires ad hoc partition techniques (e.g., see Hadjighasem et al., 2017; Balasuriya et al., 2018).

To investigate the turbulent dispersion through the network formalism, the interactions between fluid particles were determined based on mutual spatial proximity. Accordingly, if two particles come sufficiently close in space during their motion, a connection is established between them. To specify the particle proximity, we assume that a particle $i$ is connected to a particle $j$ if $i$ lies inside a reference ellipsoid centred at $j$, and vice versa (by symmetry). For instance, figure 1(a) illustrates the temporal evolution of a particle along its trajectory, and the connections with other particles enclosed in its reference ellipsoid. The ellipsoid was geometrically anchored to each particle location and it was chosen as reference geometry to take into account the anisotropy of the flow. Each ellipsoid is determined by means of its semi-axis lengths, $\boldsymbol{a}=\left(a_{x},a_{y},a_{z}\right)$, representing spatial scales of turbulent motion along each Cartesian direction. Accordingly, if the Euclidean distance between a pair of particles is less than $a_{i}$ in each direction $i=x,y,z$ (namely they are connected), then the two particles share turbulent length scales greater (or equal) than $a_{i}$. The choice of $\boldsymbol{a}$ is generally a non-trivial task which depends on the specific problem under study, such as the presence of inhomogeneities in the flow or if particles are involved in chemical or biological process (in which specific interactions occur when particles are sufficiently close). When characteristic scales are not known a priori, the issue of setting $\boldsymbol{a}$ – namely to assess the typical length scales in the flow – can be faced, for instance, by relying on turbulence spectra, correlation functions, as well as coherent structure identification techniques.

In this work, in order to illustrate the potential of the proposed approach, the semi-axes of the ellipsoid were set proportional to the average pairwise distances in the corresponding Cartesian directions. In this way, the increase of the average mutual distance between particles with time – that is mainly due to the streamwise dispersion and partly to the spanwise mixing – is taken into account. By indicating the average Euclidean distance as $\left\langle\boldsymbol{d}\right\rangle$, then $a_{i}(t^{+})=\alpha_{i}\left\langle d_{i}(t^{+})\right\rangle$, where angular brackets indicate the average over all particle pairs, $i=x,y,z$ and $\boldsymbol{\alpha}=\left(\alpha_{x},\alpha_{y},\alpha_{z}\right)$ is a set of parameters. Although $\boldsymbol{\alpha}$ may also explicitly depend on $t^{+}$, we considered $\alpha_{i}$ constant in time in order to focus only on the temporal dependency of the semi-axes $\boldsymbol{a}$, as retained in $\left\langle\boldsymbol{d}(t^{+})\right\rangle$. In particular, as the simplest case, we selected $\boldsymbol{\alpha}$ as a constant in the range $\alpha\in\left(0,1\right]$, so that the ellipsoid size in each Cartesian direction does not exceed the average distance, $\left\langle\boldsymbol{d}\right\rangle$ (this is crucial along $y^{+}$ since particles can not exceed the inter-wall height, $2H$). As a reference case, we set $\alpha=0.5$ so that the ellipsoid size in each Cartesian direction is equal to the average distance; however, similar results are obtained with other $\alpha$ values.

Since particles follow different trajectories due to turbulent motion, the spatial proximity approach results in a non-trivial time-sequence of connections. An example of a geometrical representation of the particle dynamics is shown in figure 1(b) as a contact sequence: each connection between pairs of particles (depicted as small coloured circles) is indicated as a black arc, meaning that particles are sufficiently close in space at that time. Therefore, we employed a complex network-based approach to geometrically represent and investigate particle dynamics. Complex networks are defined as graphs – namely they are made up of entities called nodes which are interconnected by links – that show non-trivial topological features (Boccaletti et al., 2006; Newman, 2018). In this work, nodes correspond to levels, $l_{i}$, i.e. spanwise-groups of particles initially at the same $y^{+}$. By doing so, information of particles starting at the same wall-normal positions is enclosed in each node. Therefore, the metrics extracted from the network increase their statistical significance because they do not represent the dynamics of a single particle.

For any time, a weight is associated to each link between a pair of nodes $\left(i,j\right)$, which takes into account the total number of connections shared by particles belonging to levels $l_{i}$ and $l_{j}$. As a result, particle dynamics is modelled by means of a time-varying network (Barrat et al., 2004; Holme & Saramäki, 2012), that is a sequence of $N_{T}=T^{+}/\Delta t^{+}$ weight matrices, defined as

[TABLE]

where $W_{i,j}=W_{j,i}$, and the binary indicator function $I_{p,q}$ is equal to $1$ if a particle $p$ lies inside the ellipsoid of a particle $q$ (or vice versa) at time $t^{+}$, and [math] otherwise. The window function $K$ is a weighting function taking into account the interaction strength, $\boldsymbol{\lambda}_{p,q}=\left(\lambda_{x},\lambda_{y},\lambda_{z}\right)_{p,q}$, between particle pairs $\left(p,q\right)$. In this work, $\boldsymbol{\lambda}$ corresponds to the pairwise Euclidean distance (evaluated in wall units) between particles (i.e., $\boldsymbol{\lambda}_{p,q}\equiv\boldsymbol{d}_{p,q}$), but other similarity or distance functions can be adopted. As discussed for the characterization of $\boldsymbol{a}$, the choice of $K$ is conditioned to the specific problem under study and the corresponding aim. For instance, when $K=1$ each link-weight exactly counts the total number of connections shared by particles belonging to pairs of nodes. An example of a time-varying network for $K=1$ is shown in figure 1(c): the particle dynamics represented in figure 1(b) is geometrized into time-evolving networks, where nodes correspond to levels and thickness of each link is proportional to its weight.

Due to its versatility (i.e., by properly setting $\boldsymbol{a}$, $\boldsymbol{\lambda}$ and $K$, as well as the number and the initial arrangement of particles), the proposed time-varying network approach is a powerful tool for turbulence analysis. In fact, the weight matrices, $W_{i,j}(t^{+})$, capture and inherit the turbulent mixing information of particles initially located at different $y^{+}$, at time any $t^{+}$. In particular, a pair of nodes that is not linked at a given time (i.e., $W_{i,j}=0$) consists of two levels whose particles are not sufficiently close in space (namely, a link is absent in the network topology). On the other hand, non-zero $W_{i,j}$ values quantify the intensity of the connection between levels, namely the extent to which particles are close in space. By doing so, nodes do not vary with time (i.e., they still represent the same levels), while link weights depend on time.

Finally, it should be noted that the network-based approach recently proposed by Padberg-Gehle & Schneide (2017) for the study of Lagrangian transport and mixing (which follows the idea by Rypina & Pratt (2017)) can be obtained from equation 1 by setting $\boldsymbol{\lambda}=\boldsymbol{d}$, $a_{x}=a_{y}=a_{z}$, $K=1$ and by checking whether a pair of particles comes sufficiently close in space at least once in the time window considered. Additionally, equation 1 can be generalized by removing the dependence on $\boldsymbol{a}$ (that is, by setting $a_{i}=\infty$), thus only keeping the window function $K$ as a similarity measure between particle pairs (e.g., as adopted in Lagrangian coherent structure approaches, see also Hadjighasem et al. (2016); Schlueter-Kuck & Dabiri (2017); Schneide et al. (2018)).

3 Results and discussion

In this section, we present the results of the time-varying weighted network approach to investigate turbulent mixing. Two configurations of window function $K$ are explored.

First, in order to show the main features of the time-varying networks in a simple case, we set $K\equiv K^{u}=1$, where superscript $u$ indicates a uniform window function (see Sections 3.1-3.3); namely, particle interaction is unweighted (binary representation). In this way, each connection is equally weighted, and the connection criterion is only driven by $\boldsymbol{a}$ – i.e., the main criterion is to check if particles come sufficiently close or not, at any time. In the second configuration (see Section 3.4), we set

[TABLE]

where $d_{p,q}$ is the Euclidean distance between particles $p$ and $q$, given that $p$ and $q$ are connected (i.e., $I_{p,q}=1$, see equation 1). If the ellipsoid is centred on a particle $p$, $D_{p,q}$ is the distance from $p$ and the point of intersection between the ellipsoid border and the straight line between $p$ and $q$ (the same criterion holds if the ellipsoid is centred on the particle $q$). In this way, the window function is bounded as $K^{w}\in\left[0,1\right]$, where the ellipsoid border represent an iso-value $K^{w}=0$, while $K^{w}=1$ is obtained if particle positions coincide (see also Hadjighasem et al. (2016); Schneide et al. (2018) for other works following an inverse distance-based weighting function). In this case, the anisotropy of the flow is explicitly considered as $D_{p,q}$ depends on the length of the ellipsoid semi-axes. By using a monotonically decreasing function for $K$, the smallest spatial scales are more weighted than the largest ones. This corresponds to the assumption that the smallest turbulent scales play a more significant role in particle dynamics.

3.1 Network structure of particle dynamics

To show the networks from the particle dynamics, in figure 2 the weight matrices at six characteristic times and their corresponding network topology are reported (for further visualizations see Movie1). Since particles are initially arranged in a uniform grid at $x^{+}=0$, each level is simply connected to levels close-in-space in the $\left(y^{+},z^{+}\right)$ plane at $t^{+}=0$, resulting in the diagonal weight matrix of figure 2(a1). For small times, the particle dynamics is led by an almost purely advective motion and the particle swarm takes the shape of a bow-like surface (at the very beginning, this surface reproduces the mean velocity profile, $U(y^{+})$, in each $\left(x^{+},y^{+}\right)$ plane). Accordingly, the weight matrix exhibits a predominant diagonal pattern, as shown in figure 2(b1). However, wall-normal turbulent mixing enables out-of-diagonal connections, and an increasing number of weak $W_{i,j}$ values between initially distant levels appears in time (e.g., levels $l_{10}$ and $l_{75}$ in figure 2(c1)). Sufficiently far enough downstream from $x^{+}=0$, all levels are interconnected with each other (e.g., see figure 2(d1)), due to the progressively enhanced transversal mixing. Nevertheless, the mean flow advection is still dominant over mixing at this stage: the initial linear diagonal structure of $W_{i,j}$ evolves into a three-square diagonal pattern, where the central square is bounded between $y^{+}\approx 500$ and $y^{+}\approx 1400$ (namely levels $l_{27}$ and $l_{74}$). This wall-normal coordinate corresponds to the $y^{+}$ value at which the mean shear, $\partial U/\partial y$, sharply decreases towards zero. In fact, while particles initially located far from the walls experience an almost zero mean shear (thus moving downstream at a high mean velocity), particles close to the wall tend to form two long tails due to the large mean shear close to the wall. The advection process makes the tails progressively stretched along the walls as time increases, highlighting the effect of the mean shear on particle swarm. Therefore, the three-square pattern emerges as a consequence of the mean shear on the particle dynamics. Finally, at some time long after, turbulent mixing becomes as effective as streamwise advection. This is first manifested as a smoothing in the three-square pattern of $W_{i,j}$ (figure 1(e1)), and later as a random-like structure (figure 1(f1)). This final state represents the Taylor dispersion regime, in which the streamwise particle distribution approaches a Gaussian distribution (see also Movie2).

The temporal characterization of particle dynamics is also highlighted by a different network topology in figure 2. For short times, the networks show a tree-like elongated structure (figure 2(a2-c2)) in analogy with the diagonal pattern of the corresponding $W_{i,j}$ (figure 2(a1-c1)). On the other hand, turbulence mixing – by enabling links between distant levels – has the effect to induce a clustered topology, that is a network geometry in which nodes tend to aggregate with each other. In fact, as illustrated in figure 2(d2), the nodes of the network at $t^{+}=2327.5$ tend to group based on their wall-normal coordinate: this topology corresponds to the three-square pattern of figure 2(d1). Finally, for large times the turbulent mixing is much more effective and the network topology develops towards a strongly aggregated pattern with a random-like layout (figure 2(e2-f2)).

3.2 Advection-mixing regimes and network robustness

In order to highlight the richness of the information contained in the weighted network as time evolves, we introduce a scalar metric, $E\equiv\sum_{i}{\sum_{j}{W_{i,j}}}/2$, that is the total number of network connections established at each time. Figure 3(a) shows the behaviour of $E$ as a function of $t^{+}$, where we see three temporal regimes. The first regime ranges in $t^{+}\in\left(0,T^{+}_{a}\right]$, where $T^{+}_{a}\approx 400$ is the time-scale in which particle dynamics is primarily led by streamwise advection. Therefore, in this first regime particles are arranged in a bow-like shape, as shown in the top inset in figure 3(a). Since the wall-normal mixing is the main responsible for the activation of distant inter-levels connections, in the first regime the total number of links is almost unchanged and only particles initially close in space are connected with each other (see also figure 2). As time increases, however, mixing progressively strengthens and $E(t^{+})$ decreases up to $t^{+}\approx 5200$. In fact, if two particles are connected (i.e., each particle lies inside the reference ellipsoid of the other one), wall-normal mixing tends to move the two particles apart in the wall-normal direction. Accordingly, particles tend to deviate from the bow-like profile (typical of the first regime) in the $y^{+}$ direction, and they come across a region in which the particles are less dense; as a result, in the second regime $E(t^{+})$ decreases as turbulent mixing progressively intensifies. Therefore, the second regime is an intermediate stage between an advection-dominant and a mixing-dominant regime. At some time long after, the advection process and the transversal mixing are balanced and the particle dynamics approaches the Taylor asymptotic state. Hence, the third regime is characterized by a nearly constant value of $E(t^{+})$, because the particle swarm spreads in the streamwise and spanwise directions without showing any spatial pattern (see the bottom inset in figure 3(a)). According to our analysis, this third regime starts at time $T^{+}_{d}\approx 5200$, corresponding to the time-scale from which Taylor’s dispersion analysis can be applied (Fernando, 2012). By comparing $T^{+}_{d}$ and $T^{+}_{\epsilon}$, we found that Taylor’s analysis holds for $T^{+}_{d}/T^{+}_{\epsilon}\gtrsim 0.37$ which is in excellent agreement with the value $0.4$ reported in literature (e.g., see Fischer, 1973). It is worth noting that, while $T^{+}_{\epsilon}$ (and $T^{+}_{d}\propto T^{+}_{\epsilon}$) can be estimated from a dimensional analysis of the conservation equation of a tracer, the value of $T^{+}_{a}$ is non-trivial. Indeed, $T^{+}_{a}$ is governed by a transient dynamics in which the magnitude of the advection and mixing terms is not easy to quantify. By means of particle geometrization into networks, however, we are able to easily distinguish the onset of both transient (i.e., $T^{+}_{a}$) and long-term regimes (i.e., $T^{+}_{d}$). It should also be noted that the chosen time resolution ($\Delta t^{+}=4.75$) is sufficiently accurate for the present analysis, as the network features smoothly evolve over time (see figure 3(a)).

To summarize, the network structures are directly affected over time by wall-normal turbulent mixing or, in general, by the interplay between mixing and advection. In particular, the effect of turbulent mixing on particle dynamics is captured by the total number of connections, $E$, which is able to reveal to which extent wall-normal mixing breaks the initial particle arrangement towards the Taylor asymptotic state.

For a fixed number of levels, $N_{y}$, the behaviour of $E(t^{+})$ – and in turn of the weighted network – basically depends on two modelling parameters: the number of particles in each level, $N_{z}$, and the constant of proportionality, $\alpha$, between the ellipsoid semi-axes and the average pairwise Euclidean distance. Besides, the DNS spatial resolution could affect the network structure, since a coarse spatial resolution implies an inaccurate velocity field and, in turn, a poor particle position resolution. However, the adopted spatial resolution is sufficient to provide a reliable particle position evaluation (Geurts & Kuerten, 2012; Kuerten & Brouwers, 2013).

To assess the robustness of the proposed approach, we first performed a parametric analysis on $N_{z}$ while keeping $\alpha=0.5$. By decreasing $N_{z}$ (keeping $N_{y}=100$ constant), the total number of particles released in the channel is reduced, thus weakening the statistical significance of the results. However, we found that – by considering a fraction $1/\delta$ of the total number of particles (with $\delta=2,...,10$) – the curve of $E(t^{+})$ scales down from the reference case shown in figure 3(a) proportionally to $1/\delta^{2}$, as expected, with a relative error below $10\%$. The effect of different $\alpha$ values, instead, is to vary the size of the reference volume for the link activation between particles. Hence, a higher (lower) value of $\alpha$ increases (decreases) the possibility that particles connect with each other. By varying $\alpha$ in the range $\left(0,1\right]$, the curve of $E(t^{+})$ is scaled but the values of $T^{+}_{a}$ and $T^{+}_{d}$ do not change (as for the parametric analysis on $N_{z}$). However, in this case, $E(t^{+})$ does not scale as $\alpha^{3}$ at any time, as one would expect, but the scaling depends on the regime. In fact, as shown in figure 3(b), the mean value of $E(t^{+})$ for $t^{+}<T^{+}_{a}$ (namely $E_{a}$) scales as $\alpha^{2}$, because in the first regime only a fraction of each ellipsoid is occupied by particles, resulting from the intersection of a bow-like surface with an ellipsoidal volume. Differently, the mean value of $E(t^{+})$ for $t^{+}>T^{+}_{d}$ (i.e., $E_{d}$) scales as $\alpha^{3}$, because in the third regime particles are spread in all directions and particles occupy the entire volume of each ellipsoid.

3.3 Characterization of turbulent mixing

The investigation of $E(t^{+})$ provides concise insights into the ensemble behaviour of all levels at each time, and it enables to distinguish different advection-mixing regimes. Nevertheless, from the weight matrices $W_{i,j}$ it is possible to extract much more detailed information. In fact – as shown in figure 2(b1-c1) where low $W_{i,j}$ values appear out of the main diagonal – the key effect of mixing is to promote the activation of links between nodes corresponding to distant levels (e.g., $l_{1}$ and $l_{100}$). Since particle geometrization into the network framework is based on the spatial proximity, the appearance of a link between two distant levels represents a peculiar event, which is important information, for instance when particles are involved in chemical reactions. In order to characterize such events, in figure 4(a) we show the temporal behaviour of $W_{i,j}(t^{+})$ for six representative pairs of nodes. If we focus on how connections between particles in the same level change over time (that is the main diagonal of $W_{i,j}$, with $i=j$), link activation starts as expected at the initial time (both for levels close to the wall, $l_{1}$, and at the center, $l_{50}$) and the weight decreases towards an asymptotic average value. For the level pair $\left(50-25\right)$ – namely particles initially started at $y^{+}\simeq 940$ and $y^{+}\simeq 465$, respectively –, link activation starts quickly in time, because particles belonging to $l_{25}$ and $l_{50}$ do not experience a strong velocity gradient and turbulent mixing enables their connection after a short time interval. A link between the pair $\left(1-25\right)$, instead, appears only at $t^{+}\approx 570$, because during the first regime particles in $l_{1}$ experience higher mean shear than particles in $l_{25}$, so only when the mixing is strong enough a link appears between them. Unexpectedly, as shown in figure 4(a), a link between the two furthest levels, $\left(1-100\right)$, appears before a link between levels $\left(1-50\right)$. This can be explained by recalling that if particles in the bow-like swarm are moved apart due to wall-normal mixing, they come across a region less dense of other particles. Therefore, particles in $l_{1}$ are more likely to connect with particles in $l_{100}$ than particles in level $l_{50}$, because when $W_{1,100}>0$ particles in $l_{50}$ are mainly located far downstream so that $W_{1,50}=0$. Therefore, link activation strongly depends on the flow features, namely both mixing and advection, and it reveals non-trivial results. In particular, since each node represents a set of $N_{z}$ particles, the activation of a link between two nodes takes into account all the pairwise connections between the $N_{z}$ particles in each node. By doing so, each link captures and highlights – through its weight – ensemble information about the dynamics of the $N_{z}$ particles in each node. It is worth noting that, for simplicity, we only explored the inter-relation between pairs of nodes, namely $2-$tuples, but the weighted networks comprise the information of all $n-$tuples of nodes.

Finally, we focused on how the presence of different levels inside each ellipsoid varies over time due to the effect of mixing. In fact, any particle belonging to a level $l_{i}$ is connected – due to spatial proximity – with a set of other particles belonging to different levels $l_{j}$, at each time. This is illustrated in figure 1(a) where the presence of different levels inside the ellipsoid of a reference particle is highlighted by different particle colors. To quantify such variability at any time, we evaluated the standard deviation, $\sigma_{l_{i}}$, of the indices, $j=1,...,N_{y}$, of levels $l_{j}$ found inside the ellipsoids of all particles belonging to a reference level, $l_{i}$. In this way, we quantify the efficiency of mixing between different levels. In figure 4(b) we show $\sigma_{l_{i}}$ as a function of time for five representative levels, $l_{i}$: due to the progressive strengthening of turbulent mixing, the values of $\sigma_{l_{i}}$ generally increase with time. Specifically, since in the first regime advection tends to move particles apart in the streamwise direction initially started at different $y^{+}$, levels close to the walls display lower values of $\sigma_{l_{i}}$ because they are unlikely to connect with other levels at higher $y^{+}$ (see top inset in figure 3(a)). For long times, instead, all the $\sigma_{l_{i}}$ approach the value of a uniform distribution, $\sigma_{l,u}$, because of strong turbulent mixing. However, only levels close to the wall (i.e., $l_{1}$ and $l_{1}00$) show $\sigma_{l_{i}}$ values above $\sigma_{l,u}$, as a consequence of their preferential connection with very distant levels (as illustrated in figure 4(a)). It should also be noted that pairs of levels at a similar distance from the wall (i.e., pairs $1-100$ and $20-80$) show analogous behaviour, since they experience similar dynamics.

3.4 Time-varying network: weighting connections, $K^{w}$

The time-varying network built by equally weighting each particle connection (i.e., $K\equiv K^{u}=1$) shows several remarkable features of the turbulent dispersion over time, as well as its effect on particle dynamics. However, in order to explicitly account for the distance between particles in a particle pair, in the present section we show the results for a non-uniform weighting function, which monotonically decreases with the Euclidean distance, as defined in equation 2. The ellipsoid is still used as a spatial proximity limit for particle connections – so that particles outside each ellipsoid are not considered since $I_{p,q}=0$ – but now particles at smaller spacing inside each ellipsoid are more weighted. It should be noted that the (binary) structure of networks for $K^{w}$ at any time is the same as shown for $K^{u}$ in figure 2. In fact, the diagonal and three-square patterns are also found for $K^{w}$ in the first and second regime, respectively, while for large times patterns do not emerge. However, by using a non-uniform window function, the values of the link weights (i.e., the link colors in figure 2) for the $K^{w}$ case tend to be more intense along the diagonal, as shortest connections are weighted more. In figure 5 we show four network metrics to further characterize the time-varying weighted networks in both configurations, $K^{u}$ and $K^{w}$. The four metrics are selected in order to progressively highlight the main network features, ranging from a local (i.e., single nodes) to a global (i.e., the entire network) point of view: the node centrality by evaluating the strength; the node pairs by evaluating the assortativity coefficient; the node triples by evaluating the clustering coefficient; and the node $n-$tuples (of higher order) by computing all the shortest paths.

Figure 5(a) shows the average strength, $\left\langle S(t^{+})\right\rangle$, of the networks as a function of time (angular brackets indicates average over all nodes). The strength of a node $i$ in a weighted network is defined as $S_{i}(t^{+})\equiv\sum_{j}{W_{i,j}(t^{+})}$, and quantifies the intensity of the relation between node $i$ and all other nodes (Barrat et al., 2004). In the case $K^{u}=1$, $\left\langle S\right\rangle$ is related to the number of particle connections, $E$, as $\left\langle S\right\rangle=2E/N_{y}$. The behaviour of $\left\langle S(t^{+})\right\rangle$ as a function of time for the two cases of $K$ is essentially the same, because (as discussed in section 3.2) the effect of turbulent dispersion on particle dynamics is fully captured by the network structure corresponding to $K^{u}=1$. However, the values of $\left\langle S(t^{+})\right\rangle$ for the case $K^{u}$ are globally higher (about two order of magnitude) than the corresponding values for $K^{w}$ (note that in figure 5(a) there are two ordinate axes), as $K^{w}$ ranges in the interval $\left[0,1\right]$. The difference in the temporal behaviour between the two $K$ configurations is more evident in the first regime, where $\left\langle S\right\rangle$ for the $K^{w}$ case (red curve in figure 5(a)) decreases more rapidly. In the first regime, advection is dominant over mixing and most of the links are present between nodes initially close in space (i.e., along the diagonal of $W_{i,j}$). However, as mentioned in Section 3.2, wall-normal turbulent mixing promotes the activation of links between initially distant nodes, as emerges from the appearance of out-of-diagonal links in figure 2(a1-c1) and the increase of standard deviation in figure 4(b). These out-of-diagonal links correspond to connections between particles that are distant in space within an ellipsoid, thus corresponding to very low $K^{w}$ values. In fact, if a particle $q$ enters inside the ellipsoid centred in a particle $p$ at a given time (with $p$ and $q$ belonging to distant initial levels), it is very likely that $q$ is close to the border of the ellipsoid of $p$ (where $K^{w}_{p,q}=0$). On the other hand, for the $K^{u}$ case, connections between particles either close or distant in space are equally weighted. Therefore, the activation of distant links (i.e., out-of-diagonal points in figure 2(b1-c1)) makes the values of $\left\langle S\right\rangle$ decrease faster for the $K^{w}$ case than for the $K^{u}$ case.

Following the concept of strength, in figure 5(b), we show the assortativity coefficient, $r(t^{+})$, which is defined as the Pearson correlation coefficient of the strength of the nodes at the ends of each link (Boccaletti et al., 2006). A network is named assortative ($0<r\leq 1$) or disassortative ($-1\leq r<0$) if nodes tend to link with other similar or dissimilar nodes (this similarity is here measured through the strength), respectively; otherwise the network is said non-assortative (i.e., $r\approx 0$). In both $K$ configurations, $r(t^{+})$ similarly decreases from approximatively $0.8$ at small times to zero at large times. In the first regime (especially at short times) nodes are primarily linked with other similar nodes (e.g., see figure 2(a1-c1)), which implies that the networks are assortative; on the other hand, for long times, the networks lose any pattern (e.g., see figure 2(e1-f1)), thus showing a non-assortative behaviour.

Figure 5(c) shows the average clustering coefficient, $\left\langle C(t^{+})\right\rangle$, for the configurations $K^{u}$ (black curve) and $K^{w}$ (red curve). The clustering coefficient quantifies the probability that two randomly chosen neighbours of a node $i$ (i.e., two nodes linked to $i$) are also neighbours, thus ranging in the interval $\left[0,1\right]$. It is formally defined as $C_{i}=N_{\Delta}(i)/N_{3}(i)$, where $N_{\Delta}(i)$ is the number of triangles and $N_{3}(i)$ is the number of triples involving node $i$, respectively (Newman, 2018). A weighted clustering coefficient takes into account the interaction intensity between nodes comprising triplets (Barrat et al., 2004). As shown in figure 5(c), $\left\langle C\right\rangle$ increases to one for both configurations as time increases: indeed, for large times particles are well mixed with each other and it is very likely that nodes form triangles, namely nodes are locally very well inter-connected. The effect of a weighted network is to enhance the local cohesiveness especially in the first regime, thus the values of $\left\langle C\right\rangle$ are the highest in the $K^{w}$ case, for which spatial proximity plays a more significant role.

Another feature of complex network structure is the concept of shortest path, namely the path with minimal cost between two nodes, where a path is an alternating sequence of nodes (and edges) considered only once (Boccaletti et al., 2006). We investigate the average path length, $L(t^{+})$, that is the average of all shortest path lengths in the network. Specifically, the average path length is $L=\sum_{i\neq j}{d^{G}_{i,j}}/(N_{y}^{2}-N_{y})$, where $d^{G}_{i,j}$ is the length (cost) of the shortest path between node $i$ and $j$ in a graph $G$ (Boccaletti et al., 2006). In this work, since the weight associated to each link, $W_{i,j}$, represents the intensity of spatial proximity between nodes, higher $W_{i,j}$ implies closer distances. Accordingly, the shortest paths for the evaluation of $L(t^{+})$ are computed by using a weighting cost $\left\langle W\right\rangle/W_{i,j}$, where $\left\langle W\right\rangle$ is the average link weight of the networks at any time (Opsahl et al., 2010). The entries of the weighted matrix, $W_{i,j}$, are normalized through $\left\langle W\right\rangle$ in order to consider the change of $W_{i,j}$ range as time increases (e.g., see the colorbar ranges in figure 2(a1-f1)). Therefore, a high weight means a low cost of the path. Figure 5(d) shows the behaviour of $L$ as a function of $t^{+}$. The networks during the first regime are able to display high $L$ values, since two nodes are reachable via a high-cost path. By inspecting the values of $d^{G}_{i,j}$ we found that outliers in $L$ for $t^{+}\in\left[20,50\right]$ are due to the nodes starting very close to the walls (namely, $l_{i}$ with $i=\left\{1,100\right\}$). In the first regime, nodes starting very close to the walls are disconnected from the other nodes in the network (see blue-dotted curve in figure 5(d)) – namely $d^{G}_{1,j}=d^{G}_{100,j}=\infty$ – because of a strong effect of mean shear on particle positions (see also figure 2(b1) or the visualizations in Movie1). Since disconnected nodes are conventionally excluded from the computation of the average path length (Boccaletti et al., 2006), the peaks in $L$ appear when nodes $1$ and $100$ are linked to the other nodes, that is at $t^{+}\approx 25$ and $t^{+}\approx 40$, respectively. From the point of view of the network topology, the re-attachment of the disconnected nodes implies high-cost shortest paths, as the values of $d^{G}_{i,j}$ for nodes $i=\left\{1,100\right\}$ are two orders of magnitude larger than $d^{G}_{i,j}$ for all the other nodes. As a consequence, very high values of $L$ are detected. It should be noted that, although the effect of the mean shear is evident in the first regime, disconnected nodes are not expected for $t^{+}\rightarrow 0$, as follows from the initial particle arrangement. For long times, instead, $L$ tends to decrease as an effect of turbulent mixing. By comparing the two weighted configurations, $K^{u}$ and $K^{w}$, the values of $L$ in the first two regimes are higher for the $K^{w}$ case because (as mentioned for the assortativity and the clustering coefficient) the window function highlights the local mixing by enhancing the spatial proximity. In the third regime, instead, both $K$ configurations approach the same constant $L$ value.

In conclusion, the analysis of the metrics at different network levels (from the single node, to node pairs and triples, as well as the shortest path of variable length) of the time-varying weighted network, is able to unveil both the main general features and the presence of extreme cases in the particle dynamics.

4 Conclusions

The proposed Lagrangian network-based approach is exemplified by means of DNS of a turbulent channel flow, where the dynamics of fluid particles is characterized by a spatial proximity criterion. The resulting time-varying weighted network is fully able to inherit the non-trivial time sequence of connections between (groups of) particles, which emerges due to turbulent motion. Indeed, we can identify in a straightforward way the characteristic regimes of particle dynamics, the appearance of peculiar events (e.g., the time of first contact between initially distant levels), as well as the intensity of wall-normal turbulent mixing (quantified by the total number of links). Accordingly, the potential of the Lagrangian-based networks is twofold, since the time-varying weighted network – represented by the weight matrices – captures in a unique framework both the qualitative spatial features of the particle swarm and the strength of turbulent mixing. The proposed complex network geometrization reveals to be robust and frame-invariant (as the Euclidean distance is used). Moreover, this approach is computationally affordable for a typical number of tracers of the order of $10^{4}-10^{5}$, and it is thus suitable for experimental techniques, such as particle tracking velocimetry, where the number of tracers is usually of the order of $10^{2}-10^{3}$ (Kim et al., 2013). Due to its versatility, particle geometrical representation into time-varying networks can easily be extended to other flows and other tracers, such as inertial particles or passive scalars. Based on present findings, Lagrangian-based networks can pave the way for a systematic network-based investigation of turbulent mixing, especially in the context of dispersion modelling.

Acknowledgements

This work was sponsored by NWO Exacte en Natuurwetenschappen (Physical Sciences) for the use of supercomputer facilities, with financial support from the Netherlands Organization for Scientific Research, NWO.

Appendix A Numerical simulation details

A DNS of a three-dimensional fully-developed incompressible turbulent channel flow was performed at $\Rey_{\tau}=950$. To this end, the continuity equation and the Navier-Stokes (NS) equations

[TABLE]

are solved. In these equations $\rho$ and $\boldsymbol{u}$ denote mass density and velocity of the fluid, $p$ the total pressure, $\boldsymbol{\omega}$ the vorticity and $t$ is time. The friction Reynolds number of the flow is kept fixed by prescribing the mean driving force per unit mass, ${\boldsymbol{f}}$, in the streamwise direction parallel to the plates. The velocity, $\boldsymbol{u}$, satisfies no-slip conditions at the two plates, i.e. $y^{+}=0$ and $y^{+}=1900$. In the other two directions, $x^{+}$ and $z^{+}$, periodic boundary conditions are applied for the velocity and the fluctuating part of the pressure. In the two periodic directions of the domain a Fourier-Galerkin approach is applied, and in the wall-normal direction $y^{+}$ a Chebyshev-tau method. The continuity equation is exactly satisfied by using the wall-normal component of the vorticity vector and the Laplacian of the wall-normal velocity component as dependent variables, instead of the three velocity components. The non-linear terms in the Navier-Stokes equations are calculated in the physical space by fast Fourier transform, with application of the $3/2$ rule in both periodic directions. The equations are integrated in time by a combination of a three-stage explicit Runge-Kutta method and the Crank-Nicolson method. In the two periodic directions $768$ Fourier modes are used, while in the wall-normal direction $385$ Chebyshev polynomials are employed. The time step used in the simulation equals $0.095\nu u_{\tau}^{-2}$, while the grid spacing are $\Delta x^{+}=7.8$, $\Delta y^{+}_{max}\approx 7.8$ (at the channel centre), and $\Delta z^{+}=3.9$.

In order to extract particle trajectories, $\boldsymbol{x^{+}}(t^{+})$, the equation $\mathrm{d}\boldsymbol{x^{+}}/\mathrm{d}t^{+}=\boldsymbol{u}\left(\boldsymbol{x^{+}}(t^{+}),t^{+}\right)$, is solved with the same explicit second-order accurate Runge-Kutta method as used in the solution of the NS equations. The fluid velocity is interpolated to the particle location by tri-linear interpolation. The accuracy of the numerical method has been assessed in previous papers (Geurts & Kuerten, 2012; Kuerten & Brouwers, 2013).