Approximate streamsurfaces for flow visualization

Abstract

Instantaneous features of three-dimensional velocity fields are most directly visualized via streamsurfaces. It is generally unclear, however, which streamsurfaces one should pick for this purpose, given that infinitely many such surfaces pass through each point of the flow domain. Exceptions to this rule are vector fields with a non-degenerate first integral whose level surfaces globally define a continuous, one-parameter family of streamsurfaces. While generic vector fields have no first integrals, their vortical regions may admit local first integrals over a discrete set of streamtubes, as Hamiltonian systems are known to do over Cantor sets of invariant tori. Here we introduce a method to construct such first integrals approximately from velocity data, and show that their level sets indeed frame vortical features of the velocity field in examples in which those features are known from Lagrangian analysis. Moreover, we test our method in numerical datasets, including a flow inside a V-junction and a turbulent channel flow. For the latter, we propound an algorithm to pin down the most salient barriers to momentum transport up to a given scale providing a way out of the occlusion conundrum that typically accompanies other vortex visualization methods.

1. Introduction

Streamsurfaces of a three-dimensional flow are two-dimensional surfaces composed of streamlines. Even a small but judiciously chosen set of such surfaces can give an effective characterization of the global topology of the velocity field. In steady flows, streamsurfaces are also invariant manifolds for the particle motion and hence frame the Lagrangian particle dynamics. For these reasons, streamsurfaces should be, in principle, the simplest tool for illustrating instantaneous features of a velocity field.

Analogously, to identify and visualize vortical features of a velocity field, Yang & Pullin (2010) generalized the notion of vortex tubes and sheets (Batchelor 2000) by defining the vortex-surface field (VSF), i.e. a smooth scalar field whose isosurfaces act as two-dimensional (2-D) invariant manifolds of the vorticity field. Initially developed for symmetric, inviscid flows (Yang & Pullin 2010), this work has been extended to capture approximations to the Lagrangian (material) evolution of VSFs in analytic viscous flows (Yang & Pullin 2011), shear flows (Xiong & Yang 2017), compressible flows (Peng & Yang 2018), transitional wall flows (Zhao, Yang & Chen 2016a,b) and in homogeneous isotropic turbulence (Xiong & Yang 2019). A different approach, based on the spherical Clebsch maps, was used by Chern et al. (2017) to visualize vortex lines and surfaces in computer graphics. However, neither VSFs nor streamsurfaces are objective (frame-indifferent; Haller 2005), which renders their experimental detection ambiguous.

In contrast, recent results on the transport of dynamically active vector fields (such as the momentum and vorticity fields) have yielded objectively defined barrier vector fields whose distinguished invariant surfaces turn out to be frame-indifferent material barriers to active transport (Haller et al. 2020). Specifically, the barrier vector fields are the velocity and vorticity Laplacian in the instantaneous (Eulerian) limit and the time-averaged pullbacks of these Laplacian fields in the Lagrangian case. The invariant manifolds of these barrier vector fields have been shown to highlight vortical features of the velocity field with increased accuracy in several 2-D and three-dimensional (3-D) flows (Haller et al. 2020; Aksamit & Haller 2022).

Irrespective of the underlying vector field, no general methodology is available for its efficient visualization via a well-placed set of invariant surfaces. This is because streamlines launched from any smooth curve form a streamsurface by definition. As a consequence, infinitely many streamsurfaces cross through any point of the flow domain. It is thus unclear which (if any) of these surfaces should be selected even locally as a representative of the vector field topology. As a consequence, flow visualization packages tend to rely on user input for seed points of streamlines and streamsurfaces.

Most related ongoing research in the scientific visualization community focuses either on the more accurate computation of streamsurfaces from select seed points (Hultquist 1992) or streamline selection for streamsurfaces based on visual optimization (Born et al. 2010; Schulze et al. 2014). A relatively recent realization is that the streamsurfaces framing the instantaneous flow behaviour most efficiently are the key invariant manifolds of the instantaneous velocity field. Local stable and unstable manifolds near stagnation points and closed streamlines indeed well illustrate the instantaneous local velocity field geometry (see, e.g. Peikert & Sadlo 2009) but generally stretch and fold globally. These globally filamented streamsurfaces then lose their ability to demarcate different flow regions efficiently, with unavoidable inaccuracies also arising in their computation (Sadlo & Peikert 2007). Closest in spirit to our work is the observation of Van Wijk (1993), who seeks streamsurfaces as level sets of a scalar function. After defining a scalar distribution within an inflow boundary, the level curves of this scalar field are propagated along streamlines into the flow. As a consequence, the resulting surfaces will generally stretch, fold and self-accumulate, resulting in filamenting streamsurfaces that also depend on the choice of the initial scalar distribution. Reviews of all these approaches in the scientific visualization community can be found, e.g. in Peikert & Sadlo (2009) and Martinez–Esturo et al. (2013).

A representative set of streamsurfaces for flow illustration should arguably include at least one surface for each typically observed streamline topology, as well as the (generally unobserved) surfaces separating these different topological classes. Finding such a set of streamlines is straightforward for 2-D incompressible flows by the existence of a first integral (conserved quantity) for the equation of instantaneous streamlines. This first integral is the streamfunction, whose set of level curves contains all typically observed streamline families as well as separatrices among them. Therefore, one can systematically scan through the one-parameter family of level curves and select those that stand out to be included in the visualization.

For the lack of a streamfunction in 3-D flows, the above program can only be carried out for streamsurfaces of integrable 3-D flows, such as steady Euler flows that do not satisfy the Beltrami property. For such flows, the Bernoulli function provides a non-degenerate first integral (Arnold & Khesin 1999) whose one-parameter family of level surfaces can be scanned and filtered to obtain the required representative set of streamsurfaces. A similar result is available for incompressible velocity fields with a volume-preserving symmetry group (Haller & Mezić 1998). For general 3-D incompressible flows, however, no first integrals will exist. This is because even steady 3-D flows can be chaotic (see, e.g. Dombre et al. 1986) and hence cannot have smooth non-trivial conserved quantities. This is equally valid for the autonomous 3-D differential equation defining the instantaneous streamlines of a generic 3-D unsteady flow, excluding the existence of a regular foliation of the flow domain by streamsurfaces that are level sets of a smooth scalar function.

Often, however, the most important parts of a flow turn out to be vortical regions filled with tubular or toroidal streamsurfaces. One cannot expect these surfaces to necessarily form a continuous family, especially the toroidal ones. Indeed, families of 2-D tori composed of streamlines typically form Cantor sets (as opposed to a continuous family) in the 3-D set of differential equations generating the streamlines (Cheng & Sun 1989). While this is also the case for classic Hamiltonian systems admitting families of invariant tori (Arnold 1989), those systems nevertheless turn out to be integrable restricted to this Cantor set of tori in an appropriate sense (see Chierchia & Gallavotti 1982; Pöschel 1982). Specifically, smooth functions can be constructed that act as first integrals over the Cantor family of tori but not in the small gaps among those tori.

Motivated by this idea of integrability over Cantor sets in Hamiltonian systems, we seek here smooth scalar functions that serve as approximate first integrals over a set of streamsurfaces forming vortical (elliptic) regions of a given vector field. As there is no widely accepted definition of a vortex, we use the term ‘vortical’ loosely to describe families of toroidal or cylindrical surfaces to which either the velocity, vorticity or barrier field is tangent. The approximate first integrals arising from this procedure will be steady for steady vector fields and time-varying for unsteady vector fields. We construct these (approximate) first integrals by seeking scalar functions whose gradient vector field is as close to being normal to the given vector field as possible. To avoid the trivial solution to this problem, we use a constrained minimization approach that does not allow for globally constant first integrals.

Our method resembles that of Yang & Pullin (2010) for the construction of VSFs in inviscid, analytic and highly symmetric flows. Ours, however, differs crucially in that we work with a grid in the physical space over which we expand the unknown approximate first integral in a Fourier series. We then only use the known values of the vector field at the gridpoints. We find that this approach results in an homogeneous linear system of equations whose unique, unit-norm least-squares solution yields the unknown Fourier coefficients of the approximate first integral. Thus, in contrast to Yang & Pullin (2010), this procedure is free from any symmetry assumptions on the first integral and does not require rewriting the homogeneous system as an inhomogeneous one under further assumptions. As we show in one of the appendices, these features of our approach significantly enhance the quality of the final solution.

Outside elliptic regions (i.e. in hyperbolic streamline domains), the streamlines are generically chaotic and hence will admit only trivial approximate first integrals. Accordingly, we expect the approximate first integrals to be nearly constant outside vortical regions, while admitting non-trivial shapes inside such regions. In those elliptic regions, level sets of the approximate first integrals will be close to streamsurfaces that form vortical features.

The result of this approach is an automated numerical visualization method that does not require the user to guess seed points in the flow for streamsurfaces in vortical regions. The simplicity and generality of the proposed method allows us to employ it to complex flows defined either analytically or through numerical data. Our examples include spatially periodic integrable and non-integrable flows, a non-periodic vortex ring flow, a V-junction flow and a fully developed turbulent channel flow. The latter flow exemplifies a case wherein exact streamsurfaces tend to obscure the visualization of the most prominent features of the barrier field. Indeed, for such a flow, vortical regions have been delineated via diagnostic tools such as the active version of the finite-time Lyapunov exponent (aFTLE), as described by Haller et al. (2020). Yet, tracking barrier streamlines originating in the neighbourhood of aFTLE ridges quickly results in the streamsurface falling apart, despite some initial vortical motion. In contrast, the structures based on the approximate first integral are able to follow closely the valleys around the aFTLE ridges allowing for a better visualization.

To obtain these results, we divide the computational domain into smaller subdomains and seek approximate first integrals in each one of them separately. In this way, we can capture the most salient structures up to a given scale without the problem of obstruction by smaller structures. These in turn can be captured by further refining the domain subdivision based on their signatures in the aFTLE field. Further, upon assuming mild convexity for the vortical structures to be extracted, we obtain families of barrier surfaces whose outermost members act as vortex boundaries. Thus, we propound a way out of the isocontour value dilemma that besets the typical vortex identification criteria used in the literature (Kim, Moin & Moser 1987; Jeong & Hussain 1995; Zhou et al. 1999), causing them to produce dissimilar structures for different isocontour values.

2. Set-up of the minimization scheme and reconstruction algorithm

Let $\boldsymbol {v}(\boldsymbol {x},t)$ be a smooth vector field defined on some open subset $U\subseteq \mathbb {R}^{3}$. The associated dynamical system for the instantaneous streamlines of this vector field at time $t$ is given by

(2.1)\begin{equation} {\boldsymbol{x}}^\prime = \boldsymbol{v}(\boldsymbol{x},t),\quad \boldsymbol{x} \in U, \end{equation}

where the prime denotes differentiation with respect to a curve-parameter $s\in {\mathbb {R}}$ along the streamline. A continuously differentiable, scalar function $H(\boldsymbol {x},t)$ is called an (instantaneous) first integral for $\boldsymbol {v}(\boldsymbol {x},t)$ at time $t$ if it is constant along each solution of (2.1), i.e. $({\partial }/{\partial s}) H(\boldsymbol {x}(s),t) = 0$. This condition implies that

(2.2)\begin{equation} \boldsymbol{\nabla} H(\boldsymbol{x},t)\boldsymbol{\cdot} \boldsymbol{v}(\boldsymbol{x},t)=0, \end{equation}

for all $\boldsymbol {x}\in U$, where $\boldsymbol {\nabla }$ denotes the gradient with respect to $\boldsymbol {x}$.

For an arbitrary vector field $\boldsymbol {v}$, no first integral will exist, in general, pointwise, i.e. for all $\boldsymbol {x}\in U$. We may, however, relax the constraint (2.2) by seeking a function $H(\boldsymbol {x},t)$ that minimizes the functional

(2.3)\begin{equation} J[H]=\frac{1}{2}\int_{U}|\boldsymbol{\nabla} H \boldsymbol{\cdot} \boldsymbol{v}|^{2}\,{\rm d}V, \end{equation}

which measures the average deviation of $H(\boldsymbol {x},t)$ from being an exact, pointwise first integral in the domain $U$ at time $t$. Any minimizer of (2.3) is called an approximate first integral.

First, let us assume that the domain $U$ is triply periodic. This allows us to expand $H$ in a Fourier series and define its modal truncation of order $N$ as

(2.4)\begin{equation} H^{{\leq} N}(\boldsymbol{x})=\sum_{0<|\boldsymbol{k}|\leq N}\hat{H}_{\boldsymbol{k}}{\rm e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}. \end{equation}

For the application to the examples in later sections, let us remark that the truncated expansion (2.4) can also be used locally on non-periodic domains, as long as we stay away from the domain boundaries. The integrand of the functional (2.3) in Fourier basis takes the form

(2.5)\begin{equation} \boldsymbol{\nabla} H\boldsymbol{\cdot}\boldsymbol{v}=\sum_{0<|\boldsymbol{k}|\leq N}\hat{H}_{\boldsymbol{k}} {\rm e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}{\rm i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{v}. \end{equation}

Because of the linearity of the gradient operator, the expression in (2.5) is linear in the unknown Fourier coefficients.

We assume that the vector field $\boldsymbol {v}$ is known on a discrete, three-dimensional grid of points which we enumerate from $1$ through $m$, where $m$ is the total number of points. Also, to target data-driven applications specifically, we work with the discretized version of (2.3) using the vector of pointwise inner products defined as

(2.6)\begin{equation} \left[\begin{array}{cccc} \boldsymbol{\nabla} H_{1}\boldsymbol{\cdot}\boldsymbol{v}_{1} & \boldsymbol{\nabla} H_{2}\boldsymbol{\cdot}\boldsymbol{v}_{2} & \cdots & \boldsymbol{\nabla} H_{m}\boldsymbol{\cdot}\boldsymbol{v}_{m}\end{array}\right]^{{\rm T}}=\boldsymbol{\mathsf{C}} \boldsymbol{\boldsymbol{h}}, \end{equation}

with $\boldsymbol{\mathsf{C}}\in {\mathbb {C}}^{m\times n}$. Here, $n$ is the number of modes used, ${\mathsf{C}}_{ij}=\exp ({{\rm i}\boldsymbol {k}_{j}\boldsymbol {\cdot}\boldsymbol {x}_{i}})\boldsymbol {k}_{j}\boldsymbol {\cdot}\boldsymbol {v}_{i}$ is the $(i,j)$ entry of $\boldsymbol{\mathsf{C}}$ and $\boldsymbol {\boldsymbol {h}}=\{ \hat {H}_{\boldsymbol {k}}|\boldsymbol {k}\in \mathbb {Z}^{3},\, 0<|\boldsymbol {k}|\leq N\}$ is the vector comprising the Fourier coefficients to be determined. Approximating a first integral through the functional (2.3) then amounts to minimizing the squared vector norm $|\boldsymbol{\mathsf{C}}\boldsymbol {\boldsymbol {h}}|^2$. To exclude the trivial solution $H\equiv {\rm const.}$ from our analysis, we add the constraint $|\boldsymbol {\boldsymbol {h}}|^2=1$.

Solving this optimization problem is equivalent to finding the eigenvector corresponding to the minimal eigenvalue of the symmetric matrix $\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{C}}^{*}\boldsymbol{\mathsf{C}}$, where $\boldsymbol{\mathsf{C}}^*$ denotes the conjugate transpose of $\boldsymbol{\mathsf{C}}$ (see Appendix B). Since the eigenvectors of $\boldsymbol{\mathsf{A}}$ are the right-singular vectors of $\boldsymbol{\mathsf{C}}$, the solution to our algorithm can also be calculated from the singular-value decomposition (SVD) of $\boldsymbol{\mathsf{C}}$.

We refrain from expanding the known vector field $\boldsymbol {v}$ into a Fourier series and we write down (2.2) explicitly for every point in the computational grid without working with the coefficients of each Fourier mode. These features distinguish our method from the one presented by Yang & Pullin (2010), as already noted in § 1. This difference will allow us to obtain unique solutions as well as apply our approach even to turbulent flows without additional assumptions.

In the absence of further constraints, the resulting approximate first integral $H$ will, generally, be a complex-valued scalar field. Denoting by $H_{r}$ and $H_{i}$ its real and imaginary parts, respectively, we have $\langle \boldsymbol {\nabla } H,\boldsymbol {v} \rangle _{\mathbb {C}} = \langle \boldsymbol {\nabla } H_{{r}},\boldsymbol {v} \rangle - \langle \boldsymbol {\nabla } H_{{\rm i}},\boldsymbol {v} \rangle \mkern 3mu {\rm i}$. This implies that by considering the expression in (2.6), we essentially optimize both the real and imaginary parts of $H$. This allows us to use $| H |$ as the approximate first integral. Alternatively, we can require $H$ to be real a priori. We discuss the latter procedure in Appendix C and show that its results are similar to those obtained without imposing this constraint.

As already noted, for generic flows, a non-trivial approximate first integral is only expected to exist in vortical regions. Outside such regions, we expect our algorithm to yield almost flat first integrals. Such a function, locally supported on several vortical regions, would generally require a large number of Fourier basis functions, which in turn would lead to numerical inefficiencies and cost. To avoid this, we work with a comparably small number of basis functions and only use the level surfaces of the emerging approximate first integral in regions where those surfaces are indeed nearly tangent to the given vector field. To identify such regions, we introduce the invariance error as

(2.7)\begin{equation} E_{A} =\frac{1}{A} \sum_{i=1}^{p} \left| \frac{\boldsymbol{\nabla} \left| H_{i} \right| \boldsymbol{\cdot} \boldsymbol{v}_{i}}{\left|\boldsymbol{\nabla} \left|H_{i}\right| \right| \left| \boldsymbol{v}_{i} \right|} \right|, \end{equation}

where $A$ is the surface area of a level set and $p$ the number of points in the level set. This type of error estimate was first introduced by Yang & Pullin (2010). In our visualizations, we will discard streamsurface candidates with invariance errors exceeding a certain threshold value.

Finally, we note that our minimization procedure can also be viewed as finding for $\boldsymbol {v}$, in the appropriate norm, the closest member of the integrable, 3-D incompressible vector field family

(2.8)\begin{equation} {\boldsymbol{x}}^\prime = \boldsymbol{\mathsf{J}}(\boldsymbol{x})\boldsymbol{\nabla} H(\boldsymbol{x},t),\quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{J}}\equiv\boldsymbol{0}, \end{equation}

with $\boldsymbol{\mathsf{J}}=-\boldsymbol{\mathsf{J}}^{\rm T}$ and $\boldsymbol {\nabla }\boldsymbol {\cdot} \boldsymbol{\mathsf{J}}$ denoting the divergence of the tensor field $\boldsymbol{\mathsf{J}}$. Indeed, all these vector fields in (2.8) share the same streamsurfaces, the level sets of $H$. Working with (2.8) directly, however, is much more demanding numerically in our experience and would also require specific assumptions on the form of $\boldsymbol{\mathsf{J}}$.

Before moving to specific examples, we note that finding an exact, pointwise first integral in (2.2) is not a well-posed problem by itself. Indeed, if $H(\boldsymbol {x},t)$ is a solution, then, for any sufficiently smooth function $F$, $F(H(\boldsymbol {x},t))$ will also be a solution due to the homogeneity of (2.2). This would not be an issue for the detection of streamsurfaces if the isocontours of $H$ and $F(H)$, expressed through the gradient of these fields, represent the same geometric and topological features for the given vector field $\boldsymbol {v}$. Unfortunately, we can construct simple counter-examples where this is not the case. For instance, if we denote by $v_x, v_y$ and $v_z$ the three components of $\boldsymbol {v}$ and assume that $v_x = 0$, then the function $G(x) H(y,z)$ will be an exact, pointwise first integral as long as $v_y \theta _y H + v_z \theta _z H = 0$. By tweaking $G$, one can easily obtain markedly different topological features resulting from the corresponding streamlines. We also refer to Pullin & Yang (2014) for more examples of first integrals with different topology for the vorticity field of the Taylor–Green flow stemming from the superposition of independent solutions to (2.2).

Consequently, our variational principle (2.3) will exhibit the same non-uniqueness issues whenever an exact, pointwise first integral is admitted by the underlying vector field. To resolve this, we will only consider approximate first integrals for which the weakest eigenvalue of $\boldsymbol{\mathsf{A}}=\boldsymbol{\mathsf{C}}^*\boldsymbol{\mathsf{C}}$ is not considered (numerically) zero. In addition, we will only retain the streamsurfaces whose topology remains unaltered when, for the same set of grid points, a larger number of Fourier modes is used allowing for small geometric refinements due to the increased accuracy of the solution. If these two conditions are met, we will consider the resulting structures as robust and they will be included in the visualization. Furthermore, we will see that the more complex a flow is, the larger the spectral gap to the second-smallest eigenvalue will be. Irrespective of this gap, however, in all the examples that follow, we will build our solution based only on the eigenvector corresponding to the smallest eigenvalue of $\boldsymbol{\mathsf{A}}$. We close this section by noting that our approach is in agreement with the findings of Xiong & Yang (2017, 2019) who showed that the construction of VSFs in turbulent flows through the use of partial differential equations leads to robust structures despite the non-uniqueness issues.

3. Approximate first integrals for explicit solutions of the Euler equations

In this section, we illustrate our minimization algorithm to construct approximate first integrals and use their level sets as approximate streamsurfaces in analytic examples.

3.1. ABC flow

As a first test case, we investigate the ABC (Arnold–Beltrami–Childress) class of flows (Dombre et al. 1986; Henon 1966), defined as

(3.1)\begin{equation} \left. \begin{gathered} \displaystyle \dot{x} = A \sin z + C \cos y, \\ \displaystyle \dot{y} = B \sin x + A \cos z, \\ \displaystyle \dot{z} = C \sin y + B \cos x, \end{gathered} \right\} \end{equation}

for $A,B,C\in \mathbb {R}$ and $(x,y,z) \in [0,2{\rm \pi} ]^3$, together with periodic boundary conditions. The right-hand side of (3.1) defines an exact steady solution to the incompressible Euler equations. For $ABC \neq 0$, the flow exhibits chaotic behaviour (Dombre et al. 1986; Henon 1966), whereas some analytic non-integrability results have been reported by Ziglin (1988, 1998).

3.1.1. Integrable case

We first analyse the ABC flow with $A=0$ for which (3.1) is completely integrable. For $BC \neq 0$, an exact, pointwise first integral is given by $H_1(x,y) = C \sin y + B \cos x$, while another independent first integral can be constructed through the use of elliptic functions (Llibre & Valls 2012). The level sets of $H_1$ are depicted in figure 1(a) on one cross-section of the computational domain for $B=\sqrt {2}$ and $C = 1$. These curves, therefore, represent the intersections of a representative set of streamsurfaces with the $z=0$ plane. These streamsurfaces are also exact invariant manifolds for the Lagrangian particle motions in this steady flow.

Figure 1. Analysis of the integrable ABC flow using a computational grid of $100^3$ points and approximately $9000$ Fourier modes. (a) Intersections of the level surfaces of $H_1$ with the $z=0$ plane. (b) Intersections of the level surfaces of the approximate first integral $H$ with the $z=0$ plane. (c) Same as panel (b) but after the removal of small-scale structures of panel (b) as well as the structures with $E_{l} > 10^{-5}$.

To test our proposed algorithm for constructing approximate first integrals, we start from a discretized version of the full 3-D velocity field (3.1) along $100$ points per direction, using approximately $9000$ Fourier modes in (2.4)–(2.6). Solving the underlying optimization problem using the algorithm in Appendix B, we obtain the results shown in figure 1(b) at the same cross-section as in figure 1(a). The numerically constructed approximate level sets match the analytic first integral closely. To measure the proximity of streamsurfaces and approximate streamsurfaces along the $z=0$ plane, we introduce a planar version of the general invariance error (2.7) by defining

(3.2)\begin{equation} E_{l} =\frac{1}{l} \sum_{j=1}^{p} \left| \frac{\boldsymbol{\nabla} \left| H_{j} \right| \boldsymbol{\cdot} \boldsymbol{\nabla} H_{1j}}{\left| \boldsymbol{\nabla} \left|H_{j}\right| \right| \left| \boldsymbol{\nabla} H_{1j} \right|} \right|, \end{equation}

where $l$ is the length of the level and $p$ is the number of points for each level set. We use this metric to remove level curves with fewer than $30$ points and those with invariance errors $E_{l} > 10^{-5}$. The results shown in figure 1(c) confirm that choosing these thresholds removes small-scale artefacts arising from numerical inaccuracies.

We perform the same analysis on the $z=0$ plane using $150$ points per direction but only approximately $1200$ Fourier modes. The results are depicted in figure 2. Again, we observe close agreement between the known first integral and the reconstructed curves. At some points, the agreement is even closer when compared to figure 1(c) due to the higher spatial resolution, even though the number of Fourier modes is significantly smaller.

Figure 2. Same as figure 1 but with a computational grid of $150^3$ points and $1200$ Fourier modes.

3.1.2. Non-integrable case

For a different set of parameter values ($A = \sqrt {3}$, $B = \sqrt {2}$ and $C = 1$), the ABC flow (3.1) is non-integrable and shows chaotic behaviour in some regions. The dynamic behaviour of trajectories for this set of parameter values is well studied, including the KAM-type tori highlighted by Poincaré maps (Dombre et al. 1986) and elliptic Lagrangian coherent structure (LCS) techniques (Blazevski & Haller 2014; Oettinger, Blazevski & Haller 2016). We use this velocity field as a benchmark to test different solution algorithms for finding an approximate first integral for a non-integrable flow. Also, this will serve as a proof of concept for finding elliptical regions.

We have already noted that the unit-norm least-squares solution to the homogeneous system of (2.6) coincides with the right-singular vector of $\boldsymbol{\mathsf{C}}$ associated with its smallest singular value or, equivalently, with the eigenvector associated with the smallest eigenvalue of $\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{C}}^{*}\boldsymbol{\mathsf{C}}$. To improve numerical stability, the SVD-based solution is preferred (Golub & Pereyra 1973). Indeed, the eigenvalue calculation requires a matrix multiplication to form $\boldsymbol{\mathsf{A}}$, which invariably squares the condition number of $\boldsymbol{\mathsf{C}}$. For comparison, we calculate both the singular-vector-based and the eigenvector-based solutions on the triply periodic box $[0,2 {\rm \pi}]^3$ with $100$ points per direction and $N=13$ (or $9170$ Fourier modes). In this setting, running the SVD algorithm of MATLAB on $\boldsymbol{\mathsf{C}}$, which is a $(10^2)^3 \times 9170$ matrix, would require an exorbitant amount of memory (more than $400$ GB), indicating that the classical SVD algorithm is not optimized for tall-skinny matrices (such as our coefficient matrix). To proceed with our comparison test, we instead follow the modified SVD method discussed in Appendix D for tall-skinny matrices.

With this modification to the SVD-based solution, the results from the two approaches for the non-integrable ABC flow are presented in figure 3 on the $z=0$, $y=0$ and $x=0$ planes. We observe that the differences between the eigenvector-based and singular-vector-based computations are marginal, indicating that the larger condition number of $\boldsymbol{\mathsf{A}}$ does not affect the results. Furthermore, we use the same three planes as Poincaré sections to integrate trajectories up to an arclength of $10^4$ from a uniform grid of $20 \times 20$ initial conditions on each section. We observe a very good agreement between the predicted structures and the intersections of the KAM surfaces with each of these sections. This is highlighted perhaps even better by the reconstructed KAM surfaces as approximate streamsurfaces in figure 4, which are to be contrasted with the $\lambda _2$-based structures in Appendix A. Since the ABC flow is a Beltrami flow, its velocity $\boldsymbol {v}$ is parallel to the vorticity $\boldsymbol {\omega } = \boldsymbol {\boldsymbol {\nabla }} \times \boldsymbol {v}$ (in fact, $\boldsymbol {v}=\boldsymbol {\omega }$); consequently, the approximate-first-integral-based tori we have constructed are also VSFs. This illustrates that our algorithm can identify VSFs in flows where the methodology of Yang & Pullin (2010) is inapplicable. Indeed, as already noted, the non-integrable ABC flow has chaotic streamlines and, thus, no symmetry assumptions regarding these streamlines can be used to accelerate the convergence rate for the optimization technique presented by Yang & Pullin (2010). Even if this rate was irrelevant, however, expanding the known velocity field in a Fourier series would result in an optimization problem with many (numerically) zero eigenvalues and, thus, infinitely many possible minimizers.

Figure 3. Analysis of the non-integrable ABC flow using a computational grid of $100^3$ points and $9170$ Fourier modes. Level sets of the approximate first integral at (a,d) $z=0$, (b,e) $y=0$ and (c,f) $x=0$. Panels (a–c) are constructed from the eigenvector of $\boldsymbol{\mathsf{A}}$ corresponding to the smallest eigenvalue, whereas panels (d–f) are produced using the SVD of $\boldsymbol{\mathsf{C}}$. The overlaid Poincaré map (black dots) on each section is based on a uniform grid of $20 \times 20$ initial conditions.

Figure 4. Two different views of the approximate streamsurfaces (level sets of the approximate first integral) closely approximate the KAM-type surfaces of the non-integrable ABC flow in elliptic regions. Also shown are iterations of the Poincaré map (black dots) on three orthogonal planes. The results were obtained using the weakest eigenvector of the positive definite matrix $\boldsymbol{\mathsf{A}}$. See also the supplementary movie 1 available at https://doi.org/10.1017/jfm.2022.992.

Upon taking a closer look at the results of figure 3, we notice that the reconstructed level sets attain their values in a longer range (i.e. $[0,3.5]$) for the larger KAM surfaces, whereas, in the vicinity of the smaller structures, they are confined to a narrow band (i.e. $[0.25,0.35]$). Here the adjectives larger and smaller are used to refer to either the area (figure 3) or the volume (figure 4) these families enclose. This is a type of overfitting that we would like to mitigate. One way to achieve this is by considering a slightly different optimization problem (see Appendix E) which resembles the one put forward by Yang & Pullin (2010). This different approach, however, turns out to be computationally intense, posing severe limitations on its use for typical grid sizes while its results are arguably of inferior quality. All in all, obtaining the solution to the proposed algorithm as the eigenvector corresponding to the smallest eigenvalue of $\boldsymbol{\mathsf{A}}$ is computationally superior to all the other techniques used and, thus, it is the one that we will follow in the rest of this article.

We conclude this section by performing a convergence analysis for different numbers of modes in figure 5. We observe that the least-squares error (as the smallest eigenvalue of $\boldsymbol{\mathsf{A}}$) approaches zero as the number of modes used in the analysis increases. Furthermore, the second smallest eigenvalue converges to the smallest one for higher modes. In Appendix B, we show that when the smallest eigenvalues are almost equal, we can construct the solution as a linear combination of the eigenvectors corresponding to these near-identical eigenvalues. Here, however, we can base our solution only on the weakest eigenvector since we have $d_1 = 0.015415$ and $d_2=0.01596$ for $N=13$. This gap between $d_1$ and $d_2$ will prove significantly larger in the following, data-based examples, confirming the uniqueness of solutions to the optimization problem.

Figure 5. Five smallest eigenvalues of $\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{C}}^{*}\boldsymbol{\mathsf{C}}$ for different numbers of modes ($2108, 3070, 4168, 5574, 7152, 9170$ modes for $N=8,9,10,11,12,13$, respectively).

3.2. Further analytic solution to the Euler equations

A set of analytic, unsteady, tri-periodic laminar solutions of the Navier–Stokes equations was put forward by Antuono (2020). Here, we consider only the steady part of these solutions, which is a Beltrami solution to the Euler equations with no known first integral. We therefore expect the streamlines of this velocity field to be chaotic and the overall dynamics to be non-integrable. Representative (approximate) streamsurfaces have not yet been constructed for this flow in the literature.

The velocity field is given by

(3.3)\begin{equation} \boldsymbol{v}=\dfrac{4\sqrt{2}}{3\sqrt{3}}\left(\begin{array}{c} \sin\left(x-\dfrac{5{\rm \pi}}{6}\right)\cos\left(y-\dfrac{\rm \pi}{6}\right)\sin(z)-\cos\left(z-\dfrac{5{\rm \pi}}{6}\right)\sin\left(x-\dfrac{\rm \pi}{6}\right)\sin(y)\\ \sin\left(y-\dfrac{5{\rm \pi}}{6}\right)\cos\left(z-\dfrac{\rm \pi}{6}\right)\sin(x)-\cos\left(x-\dfrac{5{\rm \pi}}{6}\right)\sin\left(y-\dfrac{\rm \pi}{6}\right)\sin(z)\\ \sin\left(z-\dfrac{5{\rm \pi}}{6}\right)\cos\left(x-\dfrac{\rm \pi}{6}\right)\sin(y)-\cos\left(y-\dfrac{5{\rm \pi}}{6}\right)\sin\left(z-\dfrac{\rm \pi}{6}\right)\sin(x) \end{array}\right). \end{equation}

Using a uniform grid of $20 \times 20$ initial conditions on the $y=0$ plane, we integrate trajectories up to an arclength of $10^4$ and, upon retaining their long-term behaviour from the interval $[5 \times 10^3, 10^4]$, we show the resulting Poincaré map in figure 6(a) where a plethora of KAM tori are discernible. Selecting the triply periodic box $[0,2 {\rm \pi}]^3$ and sampling it with $100^3$ points, we run our algorithm for $N=13$. Upon constructing isosurfaces for $10$ different isovalues of the resulting approximate first integral, we locate their intersections with the $y=0$ plane and superimpose them on figure 6(a). We then compute the invariance error based on (2.7) and colour the isocontours of figure 6(a) blue or red depending on whether their average error corresponds to an angle of less or more than $5^{\circ }$. Based on this, contours lying inside the chaotic sea of the Poincaré map show the largest invariance error despite our algorithm not using any knowledge of the long-term, chaotic dynamics.

Figure 6. (a) Comparison of the Poincaré map on the plane $y=0$ for the steady Euler flow (3.3) overlaid on the intersections of the tori obtained from an approximate first integral with the same plane for $N = 13$. The blue (red) isocontours depict tori whose invariance error (see (2.7)) is smaller (larger) than $5^{\circ }$ on average. (b) 95th percentile of the distance (in non-dimensional units) between trajectories emanating from the isocontours of panel (a) with the corresponding 2-D tori inside the computational box $[0,2{\rm \pi} ]^3$. The points to the left (right) of the dash–dotted line correspond to trajectories originating in the blue (red) contours of panel (a).

For every isocontour of figure 6(a), we launch trajectories following the vector field $\boldsymbol {v}$ of (3.3) until they leave the computational domain $[0,2{\rm \pi} ]^3$, and then calculate their distances from the corresponding isosurfaces. Thus, in figure 6(b), we present the 95th percentile of those distances for each streamsurface after separating them according to figure 6(a), i.e. we depict the trajectories corresponding to the blue contours in the left side of the dash–dotted line, whereas the ones for the red contours in the right side. We note the correlation between the retained (blue) isocontours of figure 6(a) and the significantly smaller percentiles of figure 6(b).

Similarly, we present the reconstructed tori for $N = 15$ or $14\,146$ modes in figure 7(a). We observe that for $N=13$, the approximate first integral captures virtually all the KAM surfaces indicated by the Poincaré map. In contrast, for $N=15$, some of the structures are captured more accurately while others are missed completely. The convergence analysis depicted in figure 8(a) shows that the least-squares error follows a declining trend as $N$ grows. This prompts us to consider an error measure similar to the one in (2.7), defined as

(3.4)\begin{equation} E_m=\frac{1}{m}\sum_{i=1}^{m}\left| \frac{\boldsymbol{\nabla} \left| H_{i} \right| \boldsymbol{\cdot} \boldsymbol{v}_{i}}{\left| \boldsymbol{\nabla} \left| H_{i} \right|\right| \left| \boldsymbol{v}_{i}\right| } \right|. \end{equation}

In this expression, the summation is taken over all the grid points, providing an estimate for the mean invariance error of the entire solution. This allows us to make a direct comparison among solutions corresponding to different numbers of modes. Specifically, excluding points that lie in the vicinity of fixed points of either $\boldsymbol {v}$ or $\boldsymbol {\nabla } | H |$, we show the dependence of the invariance error $E_m$ on the number of Fourier modes used in our algorithm in figure 8(b). We observe that the error attains the minimum value for $N=13$, in agreement with what is inferred from figure 6(a).

Figure 7. (a) Same as in figure 6(a) with the tori reconstructed for $N = 15$. (b) Closeup view on a region filled with two families of invariant tori. Overlaid on the Poincaré map are the tori obtained from an approximate first integral for $N = 19$.

Figure 8. Numerical details of the approximate first integral calculation for the steady Euler flow (3.3). (a) Five smallest eigenvalues of $\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{C}}^{*}\boldsymbol{\mathsf{C}}$ and (b) normalized error estimate (3.4) for different numbers of modes.

Moreover, to mitigate minor discrepancies between the reconstructed tori and those outlined by the Poincaré map of figure 6(a), we increase the number of Fourier modes to $28\,670$ ($N=19$), while keeping the same set of grid points. A closeup view of a region filled with invariant tori in figure 7(b) confirms that there is a close agreement with the tori obtained from an approximate first integral. There is, however, a trade-off between the increased accuracy and the ensuing computational burden that an end user must consider.

We also note that, while for $N=15$ we have $d_1 \approx d_2 = 839.8079$, for $N=19$ we have $d_1 = 204.5165$ and $d_2 = 218.2078$, further corroborating that our approach leads to a unique solution. Finally, we conclude this section with figure 9(a) showing 3-D rendered images of the three Poincaré maps on the planes $x=y=z=0$, figure 9(b) showing representative streamsurfaces obtained as level surfaces of an approximate first integral for $N=13$ and figure 9(c) showing the superimposition of panel (a) on panel (b), which confirms the close agreement between the expected and reconstructed structures.

Figure 9. Results for the steady Euler flow (3.3). (a) Poincaré maps on $x=y=z=0$. (b) Streamsurfaces approximating the KAM surfaces of (3.3). (c) Panel (a) superimposed on panel (b).

3.3. Hill's spherical vortex

We now turn to a spatially non-periodic, integrable flow given by Hill's spherical vortex (Hill 1894). The axial-symmetric (approximately the $z$-axis) streamfunction for Hill's solution to the Euler equations is

(3.5)\begin{equation} \psi(h,z)=\left\{ \begin{array}{@{}ll} \dfrac{3}{4}U_{0}h^{2}\left(1-r^{2}\right), & r\leq1,\\ -\dfrac{1}{2}U_{0}h^{2}\left(1-\dfrac{1}{r^{3}}\right), & r>1, \end{array} \right. \end{equation}

where $h$ is the distance from the axis of symmetry ($h^2 = x^2 +y^2$) and $r^{2}=h^{2}+z^{2}$. Using the relations $u_{h}=-({1}/{h})({\partial \psi }/{\partial z})$ and $u_{z}=({1}/{h})({\partial \psi }/{\partial h})$, we obtain the corresponding velocity field in Cartesian coordinates

(3.6)\begin{equation} \displaystyle \boldsymbol{v}(x,y,z)=\left\{ \begin{array}{@{}ll} \dfrac{3}{2}U_{0} \left(xz, \ yz, \ 1-(r^2+h^2)\right) , & r\leq1,\\ \dfrac{3}{2}U_{0} \left(\dfrac{xz}{r^5}, \ \dfrac{yz}{r^5}, \ \dfrac{2}{3} \dfrac{1}{r^3} - \dfrac{h^2}{r^5}\right), & r>1. \end{array} \right. \end{equation}

As for periodic domains, we can compute a Fourier expansion on a bounded subdomain, bearing in mind that the convergence of the partial Fourier sum will be slow in general (Gottlieb & Shu 1997). To mitigate this issue, we will have to consider a sufficiently high number of modes. Furthermore, due to the well-known Gibbs phenomenon, there will be sizeable spurious oscillations in the approximate first integral near the box boundary that do not diminish after an increase in the number of modes. To damp this effect in our algorithm, we discard all the reconstructed surfaces that fall within $5\,\%$ of the box size in all three directions.

The velocity field of (3.6) has toroidal streamsurfaces inside the spherical domain $\{\boldsymbol {x} \in \mathbb {R}: |\boldsymbol {x}| \leq 1\}$ (see figure 10a). We reconstruct these streamsurfaces over the computational domain $[-2,2]^3$ with $60$ points per direction, running our algorithm for different numbers of Fourier modes with $U_{0}=1$. In some cases with $N<15$, the non-periodic nature of the velocity field results in the breakdown of the reconstructed toroidal structures. For $N\ge 15$, however, we obtain fully symmetric solutions. Specifically, for illustration purposes, we work with $N=17$ (or $20\,478$ Fourier modes) which yields a unique solution corresponding to the smallest eigenvalue of $\boldsymbol{\mathsf{A}}$, $d_{1} = 23.5169$ ($d_2 = 25.3118$). Based on this, we create $10$ isosurfaces in the interval $\big[ |H|_{min},|H|_{max} \big]$. Upon locating the intersections of these surfaces with eight radially equidistant planes, we launch trajectories corresponding to these intersections and compute the pointwise distance of the solution curves to the reconstructed surfaces. The results (in terms of percentiles) are presented in figure 10(b). A comparison between the reconstructed streamsurfaces and indicative solution curves is given in figure 11.

Figure 10. (a) Level sets of Hill's streamfunction on the $x=0$ plane. (b) 25th, 50th, 75th and 95th percentile of the pointwise distances (in non-dimensional units) between solution curves and $10$ different reconstructed streamsurfaces for $N=17$.

Figure 11. (a) Streamsurfaces of figure 10 that have a 95th percentile less than $0.05$ for $x>0$. (b) Solution curves of (3.6) for $5$ different points lying on the outer streamsurface of panel (a).

3.4. Flow inside a V-junction

Next we investigate the flow inside a V-junction, as depicted in figure 12. Despite the simple geometry, recent experiments have suggested that pumping a particle-laden fluid into such a configuration allows light particles, such as gas bubbles in water, to be permanently trapped in the junction (Vigolo, Radl & Stone 2014). This phenomenon arises for a wide range of Reynolds numbers and for various junction angles (Ault et al. 2016). Here we consider one such flow with junction angle $70^{\circ }$ and ${Re} = ( \bar {u}L/\nu ) = 230$ (Shin, Ault & Stone 2015), where $\bar {u}$ is the average inlet flow speed, $L$ is the side length of the square channel and $\nu$ is the kinematic viscosity.

Figure 12. Velocity streamlines, colour coded with their magnitude, emanating from the inlet of the V-junction and leaving the domain from the two outlets. For ${Re} = 230$ and a junction angle of $70^{\circ }$, four symmetric vortex-breakdown bubbles are portrayed in grey.

We use a finite-volume solver from the OpenFOAM library to obtain a steady solution to the 3-D incompressible Navier–Stokes equation (Weller et al. 1998). The same numerical solution has recently been analysed using methods from dynamical systems theory, which revealed large, anchor-shaped trapping regions for light particles (Oettinger et al. 2018). These trapping regions, however, have been invariably linked to bubble-type vortex breakdown structures in the fluid flow, which are formed downstream in the junction (Ault et al. 2016). These structures are depicted as dark grey blobs in figure 12, obtained from a careful advection of streamlines from the vicinity of known stagnation points. Their construction, therefore, is by no means automated and assumes a detailed knowledge of the streamline geometry.

Each vortex breakdown bubble is demarcated by the 2-D stable and unstable manifolds of two saddle-type fixed points. In a generic 3-D flow, these two manifolds do not coincide because that configuration would not be structurally stable. Instead, they are expected to intersect transversely. In such a scenario, the original streamsurface breaks down allowing fluid particles from the main stream to be entrained into the bubble and, conversely, particles to return from the bubble to the main stream (Holmes 1984; Peikert & Sadlo 2007). Sotiropoulos, Ventikos & Lackey (2001) showed that a very careful mesh refinement is required to reveal the splitting of the manifolds. Here, we will only be interested in reconstructing a closed streamsurface that manifests minimal fluid exchange (see Peikert & Sadlo (2007) for a discussion of this surface).

We select as computational domain the box depicted in figure 13 with $110, 90, 80$ points in the $x,y$ and $z$ direction, respectively. This box is chosen so that the bubble-like vortical structure lies approximately in its centre. The dimensions of the box are $L_{x} = 1\, \textrm {m}$, $L_{y} = 0.6\, \textrm {m}$ and $L_{z} = 0.4\, \textrm {m}$. This procedure incorporates a priori information approximating the rough location of the streamsurface of interest. We will see in the next section how the same algorithm can be extended to uncover a priori unknown structures in a turbulent flow.

Figure 13. Different views of the computational box used to construct an approximate first integral for one of the bubble-like structures in the V-junction flow.

We use our algorithm with $N=15$ or $14\,146$ Fourier modes and show a 2-D cross-section at the middle of the computational domain in figure 14(a). For the smallest eigenvalue of $\boldsymbol{\mathsf{A}}$, $d_1 = 58.1538$ ($d_2 = 376.9696$), a very pronounced circular structure is clearly visible approximately in the middle of the domain bordering an otherwise flat landscape. To obtain the full, 3-D reconstruction, we use $40$ equidistant values between $| H |_{min}$ and $| H |_{max}$ and extract $28$ different level surfaces of the approximate first integral. We then sort them in decreasing order based on the volume they enclose.

Figure 14. Results for an approximate first integral in the V-junction flow, framing one of the elliptical vortex regions. (a) Approximate first integral distribution on a plane which coincides with the middle, in the $x$ direction, of the computational domain presented in figure 13. (b) Normalized invariance error as a function of the extracted isosurfaces sorted in descending order with respect to their volume.

Plotting this error estimate for the extracted isosurfaces (figure 14b) yields a global minimum for the surface $i=25$ and a local minimum for the surface $i=15$. Rendered depictions of the extracted isosurfaces corresponding to these two minima are shown in figure 15. We note the close agreement between the reconstructed structures and the structure delineated by the judiciously chosen streamlines. This is in stark contrast to the structures suggested by the $Q$-criterion as demonstrated in Appendix A.

Figure 15. Streamsurfaces as level sets of an approximate first integral in the V-junction flow, corresponding to one local $(i = 15)$ and the global minimum $(i = 25)$ of the normalized error $E_{A}$ depicted in figure 14(b). The duct is cut transversely to help the visualization.

4. Momentum transport barriers in a turbulent channel flow

As we have mentioned in § 1, the shape of streamsurfaces is observer-dependent: their geometry changes under general time-dependent rotations and translations of the observer. Therefore, unless the flow has a distinguished frame in which streamsurfaces coincide with material surfaces, streamsurfaces simply highlight features of the velocity field in a given frame, as opposed to intrinsic, observer-indifferent features of the flow of fluid particles. For studies seeking to be consistent with observed flow physics, the latter features are relevant (Haller 2021). This is because flow visualization experiments with dye particles highlight material (and hence objectively defined) transport barriers, which generally differ substantially from streamsurfaces in unsteady flows. At the same time, classic studies focused on momentum, energy or vorticity transport are inherently frame-dependent by the frame-dependence of all these Eulerian fields.

To reconcile these two objectives, Haller et al. (2020) developed a theory of objectively defined barriers to the transport of active vector fields, such as the vorticity and momentum. Aksamit & Haller (2022) used this theory to locate instantaneous (Eulerian) frame-indifferent momentum transport barriers in 3-D turbulent channel flows. Active barriers turn out to be distinguished streamsurfaces of appropriately defined steady, 3-D incompressible vector fields (Haller et al. 2020). Specifically, momentum transport barriers of a Navier–Stokes velocity field $\boldsymbol {v}(\boldsymbol {x},t)$ at time $t$ are streamsurfaces of the barrier equation

(4.1)\begin{equation} \boldsymbol{x}^{\prime}(s)=\boldsymbol{\Delta} \boldsymbol{v}\left(\boldsymbol{x}(s),t\right), \end{equation}

with $s\in {\mathbb {R}}$ denoting a parametrization of streamlines forming the streamsurfaces and $\boldsymbol{\Delta}$ the Laplace operator. For an incompressible velocity field, we further note that $\boldsymbol{\Delta} \boldsymbol {v} = -\boldsymbol {\nabla } \times \boldsymbol {\omega }$ holds.

Haller et al. (2020) and Aksamit & Haller (2022) detected distinguished streamsurfaces of the momentum barrier equation (4.1) using active versions of some of the passive hyperbolic and elliptic LCS diagnostics reviewed by Haller (2015). These Lagrangian calculations involve arrays of trajectories and return diagnostic scalar fields for visual inspection rather than explicit streamsurfaces families. In the following, we will use our approach for finding approximate first integrals to obtain vortical momentum barriers as level surfaces of an approximate first integral for (4.1).

4.1. Numerical dataset

We study the 3-D incompressible, turbulent channel flow which can be found here. The friction Reynolds number is ${Re}_{\tau } = u_{\tau }h/ \nu = 150$, where $u_{\tau }$ is the friction velocity, $h$ is the channel half-height and $\nu$ is the kinematic viscosity. We denote by $x$, $z$ and $y$ the streamwise, spanwise and wall-normal directions, respectively. The computational domain is $L_x = 2.5{\rm \pi} h$ long and $L_z = {\rm \pi}h$ wide.

The number of Fourier modes used in the simulation was $192$ in both the streamwise and the spanwise direction. The number of points in the wall-normal direction is $194$, inhomogeneously spaced so that the grid becomes more refined closer to the walls. No-slip boundary conditions were applied at the two channel walls and the governing equations were integrated forward in time with a constant pressure gradient $\textrm {d}p/{\textrm {d}\kern0.06em x} = -1$ enforced so as to drive the flow through the channel. Ten thousand velocity field snapshots were stored at multiples of the simulation time step $\boldsymbol{\Delta} t = 0.001$ once the flow reached a statistically stationary state. We will identify the $5000$th snapshot of this time series with the time $t = 0$.

4.2. Momentum barrier extraction

We seek to uncover objectively defined instantaneous momentum transport barriers corresponding to the time $t = 0$ as specific streamsurfaces of the vector field $\boldsymbol{\Delta} \boldsymbol {v}(\boldsymbol {x},0)$ in (4.1). First, as illustration of prior results by Haller et al. (2020) on this problem, we show the FTLE field computed for $\boldsymbol{\Delta} \boldsymbol {v}(\boldsymbol {x},0)$ (called active FTLE or aFTLE) from a grid of $800 \times 1000$ initial conditions uniformly placed in the wall-normal and spanwise directions, respectively, up to $s = 10^{-3}$ (see figure 16a). This plot indicates the signatures of two larger barriers to momentum transport, with the first located approximately around $[z/h,y/h] = [0.75,1.25] \times [1.5,1.75]$. The second, a mushroom-type barrier, is located around $[2.25,2.75] \times [1.25,1.75]$. A plethora of other smaller-scale structures are also present in the aFTLE plots.

Figure 16. (a) Active FTLE (aFTLE) for the momentum barrier field (4.1) at $t = 0$ over a 2-D cross-section of the channel at $x/h = 2$. (b–d) Projections of the approximate first integral (black lines) on the same cross-section using as computational domains for the analysis the (red) boxes $[1,3] \times [0,2] \times [0,{\rm \pi} ]$, $[1,3] \times [0,2] \times [0.5,2.5]$ and $[1,3] \times [1.25,2] \times [0.4,1.4]$, respectively, superimposed on the aFTLE landscape.

We use three different computational domains to illustrate the potential of the algorithm presented in § 2: the boxes $[1,3] \times [0,2] \times [0,{\rm \pi} ]$, $[1,3] \times [0,2] \times [0.5,2.5]$ and $[1,3] \times [1.25,2] \times [0.4,1.4]$ with $80 \times 85 \times 80$ points in the streamwise, wall-normal and spanwise direction, respectively. In each of these boxes, the points are evenly spaced and a tri-linear interpolation scheme is employed to obtain the necessary $\boldsymbol{\Delta} \boldsymbol {v}$ values. Moreover, we perform another rescaling of the dummy time $s$ to pointwise normalize the right-hand side of (4.1). This eliminates the high norms of $\boldsymbol{\Delta} \boldsymbol {v}$ near the wall in the expression (2.6). If this normalization is not used, our algorithm tends to miss the turbulent regions close to the wall and touches more on the central part of the channel.

For $N = 13$ $(9170$ Fourier modes$)$, $20$ isovalues of the approximate first integral are depicted in figure 16(b–d) on the $x/h = 2$ plane. We note the progressive improvement in the reconstruction of the first structure discussed above, as the computational domain starts to close in on it. This observation prompts us to use a greedy algorithm compartmentalizing the domain into smaller, overlapping domains. The size of these domains will be dictated by the scale of the structures to be extracted.

Our algorithm so far generates families of streamsurfaces. We now seek to locate outermost members of nested elliptical barrier families. To this end, for each of the overlapping computational boxes, we use $n_{iso}$ different values to produce isosurfaces in the interval $[ | H |_{min},| H |_{max}]$. We then group the resulting isosurfaces into foliations (see figure 20(f) for two such foliations) or, more generally, families by finding the ones that either lie entirely inside others or have an intersection volume above a certain threshold, respectively. Subsequently, we discard the foliations or families that have a number of members less than a fixed percentage of $n_{iso}$. We note that for the isosurface extraction, we use MATLAB's built-in function which is based on the ‘Marching cubes’ algorithm (Lorensen & Cline 1987). Other techniques generating isosurfaces (see, e.g. the concept of contour trees (Carr, Snoeyink & Axen 2003)) are known to produce robust structures in a more efficient fashion and could certainly be used as viable alternatives.

We also follow Haller et al. (2016) and define the convexity deficiency of a closed surface as the ratio of the volume between the surface and its convex hull to the volume enclosed by the surface. Here, however, we will only discard the most non-convex surfaces.

We summarize the main steps of the computations we described previously in Algorithm 1. In the next section, we will use this algorithm to uncover momentum transport barriers tied to different scales in the turbulent channel flow.

Algorithm 1 Extraction of instantaneous barriers to momentum transport

4.3. Results

4.3.1. Channel partition into large subdomains

To illustrate our algorithm, we use the entire computational domain of the direct numerical simulation (DNS) described in § 4.1. First we partition it in the following overlapping subdomains $[x_{1_i},x_{2_i}] \times [y_{1_j},y_{2_j}] \times [z_{1_k},z_{2_k}]$ with $x_1 \in \{-0.5,0.5,\ldots,5.5\}$, $x_2 \in \{1.5,2.5,\ldots,7.5\}$, $y_1 \in \{0,0.5,1\}$, $y_2 \in \{1,1.5,2\}$, $z_1 \in \{0,1,2\}$, $z_2 \in \{1.5,2.5,3.5\}$, $i = 1, \ldots, 7$, $j = 1, 2, 3$, $k = 1, 2, 3$. For each of these domains, we use a grid of $70 \times 75 \times 70$ evenly spaced points on which we compute the normalized velocity Laplacian. For the approximate first integral, we use $N=13$ (or $9170$ Fourier modes). This comes as a necessary trade-off between computational time and surface accuracy after we observed that the reconstructed structures do not vary significantly when we use $N=12$ or $N=14$. With it, we extract $n_{iso} = 40$ isosurfaces in every domain. Of them, we discard the surfaces described in steps $5$ and $6$ of Algorithm  by employing $p_b = 2\,\%$, $p_f = 10\,\%$ and $d_{max} = 20\,\%$.

The results of this computation are shown in figure 17. We observe that the entire channel is populated by a host of different momentum transport barriers, the majority of which have a clear quasi-streamwise direction. These structures are reminiscent of those investigated in studies of wall-bounded flows (Robinson 1991). The structures we obtain, however, are more evenly scattered throughout the channel, appearing not only near the channel walls but also in the more quiescent central part. This is in stark contrast to typical predictions from classic vortex criteria for these flows (Hunt, Wray & Moin 1988; Chong, Perry & Cantwell 1990; Jeong & Hussain 1995). This phenomenon, i.e. structures penetrating into and spanning the bulk flow region, has already been noted in the literature (Haller et al. 2020; Aksamit & Haller 2022) by analysing the 2-D signatures of these structures via diagnostic fields. Here, however, we explicitly construct these structures via the streamsurfaces surrounding them.

Figure 17. Instantaneous barriers to momentum transport in the 3-D turbulent channel flow, derived using the larger partition described in § 4.3.1.

4.3.2. Channel partition into small subdomains

We now employ smaller partitions of the channel, using the overlapping subdomains $[x_{1_i},x_{2_i}] \times [y_{1_j},y_{2_j}] \times [z_{1_k},z_{2_k}]$ with $x_1 \in \{-0.5,0,\ldots,6.5\}$, $x_2 \in \{0.5,1,\ldots,7.5\}$, $y_1 \in \{0,0.5,1\}$, $y_2 \in \{1,1.5,2\}$, $z_1 \in \{0,1,2\}$, $z_2 \in \{1.5,2.5,3.5\}$, $i = 1, \ldots, 15$, $j = 1, 2, 3$, $k = 1, 2, 3$, while keeping all the other parameters the same as in § 4.3.1.

The results for this computation are shown in figure 18. In this case, we again observe momentum transport barriers that tend to align with the streamwise direction but in larger numbers, when compared with figure 17. Another notable difference is that the reconstructed structures show generally smaller scales than before. This demonstrates that our domain partition algorithm acts as a filter of various scales. This provides a natural way out of the occlusion quandary that besets many visualization approaches while retaining the algorithm's capacity to capture the smallest structures, delineated in the aFTLE plots, by using finer partitions. In the next subsection, we will demonstrate this capacity by pinpointing structures that the active diagnostics-based methods in Haller et al. (2020) and Aksamit & Haller (2022) would likely miss. We conclude this section with figure 19 depicting a well-formed spectral gap between the smallest and the second smallest eigenvalue of $\boldsymbol{\mathsf{A}}$ for all the computational boxes. Finally, we carry out an unsteady barrier analysis for times varying over the interval $[t_{0},t_{1}] = [0,0.5]$, as presented in the supplementary movie 2. In it, the fact that some barrier surfaces seem to appear, then disappear and perhaps reappear at a later time is a direct consequence of the flat convexity deficiency threshold we allow for these surfaces (20 %) throughout the channel. We expect this phenomenon to be resolved when the purely Lagrangian approach of Haller et al. (2020) is combined with our algorithm in follow-up work.

Figure 18. Instantaneous barriers to momentum transport in the 3-D turbulent channel flow, derived using the smaller partition described in § 4.3.2. The supplementary movie 2 of the supplementary material portrays the extracted structures for the time interval $[t_{0},t_{1}] = [0,0.5]$.

Figure 19. Distribution of the smallest (blue) and the second smallest (red) eigenvalues of $\boldsymbol{\mathsf{A}}$ over the different computational domains used for the barriers surfaces presented in figure 18.

4.3.3. Uncovering hidden structures

We now focus on the region we highlighted in figure 16(d). The signature of a vortical structure is evident from the aFTLE landscape and it is corroborated by the isolines of an approximate first integral. Zooming in on the vicinity of this structure in figure 18, however, reveals the existence of another structure (see figure 20), showing no imprint on the aFTLE landscape.

Figure 20. Two branches of a mushroom-like, objective vortex (in grey) captured by Algorithm 1 in the 3-D turbulent channel flow. The branches are superimposed on 2-D cross-sections of the aFTLE field at $x/h = 1, 1.4, 1.6, 1.9$ and $2.1$ in panels (a), (b), (c), (d) and (e), respectively. (f) A transverse cut at $x/h = 1.7$ revealing the foliations of structures constituting the two branches of the mushroom-like vortex.

To investigate this, we populate the neighbourhood around these two structures with 2-D cross-sections and compute the aFTLE field for each of them. At $x/h = 1$ (figure 20a), we note two mushroom-like pairs of structures, one larger and one smaller, but neither of the two display prominent closed regions that could signal a vortex. As we move towards larger $x/h$ values, however, the aFTLE topography changes drastically. At $x/h = 1.4$ (figure 20b), the larger structure starts to fall apart before it completely disintegrates in the following cross-sections. In contrast, the smaller structure develops a salient closed loop, delineated by an aFTLE ridge, which is coincident with the largest of the reconstructed surfaces. In the next two cross-sections (figure 20c,d), the smaller of the two reconstructed surfaces comes into play as we further observe both parts of the vortex pair following closely the aFTLE ridges downstream. Finally, in the last cross-section (figure 20e), the vortex pair pattern is broken up in agreement with the reconstructed vortex pair coming to a halt.

This example epitomizes two notable features of Algorithm 1. First, using finer partitions results in the reconstruction of smaller structures consistent with the well-known hierarchy of coherent structures in turbulence. Second, this reconstruction takes place in an almost automated fashion with minimal reliance on user-defined parameters and with no need to advect trajectories of the barrier equation.

5. Conclusions

We have introduced a minimization principle to find an instantaneous approximate first integral of a given 3-D vector field. Level surfaces of this approximate first integral are expected to provide automated approximations to vortical streamsurfaces highlighting the instantaneous elliptic regions of the vector field. We have also proposed and tested a numerical algorithm for solving this minimization problem by means of total least squares and ridge regression. This algorithm has performed well on various solutions of the steady 3-D Euler equations on triply periodic and non-periodic domains. For steady velocity fields, the parametrized family of streamsurfaces rendered by our approach approximate elliptic (toroidal or cylindrical) LCSs without the need to advect a large number of trajectories required by Lagrangian methods (Haller et al. 2016). We have also illustrated on a 3-D junction flow the applicability of our algorithm for velocity fields defined as numerical datasets produced from computational fluid dynamics simulations.

We have additionally used the described algorithm to extract objectively defined instantaneous momentum transport barriers in a 3-D turbulent channel flow in an almost-automated fashion. Such barriers are streamsurfaces of an incompressible barrier equation defined by the Laplacian of the velocity field (Haller et al. 2020; Aksamit & Haller 2022). We have extracted vortical momentum transport barriers without the need to advect arrays of trajectories and found such barriers across multiple spatial scales.

The generality of the presented algorithm provides a fertile ground for its use in the visualization of the important features in a number of diverse applications. A straightforward example is the barrier equation derived by Haller et al. (2020) for the transport of vorticity. Furthermore, the purely Lagrangian approach of Haller et al. (2020) could be combined with our algorithm to produce LCSs for unsteady simulations. This method, unlike the instantaneous approach used for the production of supplementary movie 2, will yield smoothly varying structures that are also experimentally verifiable. Another application of our algorithm could be the visualization of the magnetic field resulting from magnetohydrodynamic turbulence simulations (Biskamp 2003). Future research can also apply the algorithm to direction fields like those emerging from studies of 3-D LCSs (Oettinger & Haller 2016), to 3-D material barriers to diffusive transport (Haller, Karrasch & Kogelbauer 2018) or to uncover LCSs acting as barriers in diverse applications (Martínez et al. 2021). Such studies could extend the automated extraction of 2-D LCSs and transport barriers (Katsanoulis & Haller 2019; Katsanoulis et al. 2020) to 3-D flows, facilitating their integration to numerical simulation codes.

Alternative solutions to the same optimization problem for an approximate first integral are certainly viable to investigate in future work. These could involve different basis functions (such as Chebyshev polynomials, wavelets or radial basis functions) instead of Fourier expansions. Techniques to mitigate the Gibbs phenomenon tied to non-periodic domains (Gegenbauer polynomials) (Gottlieb & Shu 1997) are also feasible to consider. Similarly, the effect of using Chebyshev nodes, instead of the uniform grids we considered here, to reduce the impact Runge's phenomenon has near the boundaries could be examined.