Lagrangian Transport Through Surfaces in Compressible Flows

07/03/2017 ∙ by Florian Hofherr, et al. ∙ Technische Universität München 0

A material-based, i.e., Lagrangian, methodology for exact integration of flux by volume-preserving flows through a surface has been developed recently in [Karrasch, SIAM J. Appl. Math., 76 (2016), pp. 1178-1190]. In the present paper, we first generalize this framework to general compressible flows, thereby solving the donating region problem in full generality. Second, we demonstrate the efficacy of this approach on a slightly idealized version of a classic two-dimensional mixing problem: transport in a cross-channel micromixer, as considered recently in [Balasuriya, SIAM J. Appl. Dyn. Syst., 16 (2017), pp. 1015-1044].



There are no comments yet.


page 14

Code Repositories


Flux Integration by Donating Regions in 2D

view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Advective transport of conserved quantities is a fundamental physical process. Its quantification is therefore of importance in a broad variety of applications. The classic approach proceeds as follows: first, a surface111We use the term surface as the short form of hypersurface, i.e., a codimension-one submanifold of space or space-time. In three dimensions, this refers to classic surfaces, in two dimensions to curves. of interest is specified, and second, the flux density is integrated over this surface. The flux density is composed of (i) the current concentration or density of the scalar quantity, and (ii) the normal component of the velocity, i.e., the component which is responsible for transport across the surface. This approach can be viewed as Eulerian, since it looks at transport purely from the spatial perspective of the surface, ignoring which particles contribute to the total transport.

This approach comes with drawbacks when complicated geometries need to be treated [20, 21], or when a restriction of the flux computation to certain sets of particles, material sets, is of interest. The latter is relevant in the determination of the physical relevance of coherent transport, i.e., transport by (Lagrangian) coherent structures [11, 13], see [8, 24] for two transport studies on the global ocean. As was abstractly argued in [12], and as we demonstrate in section 5.3.2, for flux integration restricted to material sets it is advantageous to convert the problem of flux integration over surfaces into a Lagrangian integration problem, because the implementation of the material restriction is straightforward. The conversion of the flux integration problem has been achieved rigorously for volume-preserving flows in arbitrary finite dimensions in [12] by the second author. For such problems, even though theoretically neat and computationally superior, this approach takes an unusual perspective and may therefore seem challenging. The aims of this paper are twofold: (i) to generalize the flux integral conversion formula to non-volume-preserving flows, and (ii) to show how to apply the Lagrangian transport methodology.

This paper is organized as follows. In the rest of this Introduction, we briefly recall some terminology. In section 2 we formulate the problem of flux computation and prove our main result, theorem 1, which solves the donating region problem posed by Zhang [21]; cf. also Zhang’s prior work [22, 23]. In section 3 we specify the general theory to two-dimensional flow problems and provide an algorithmic overview. Section 4 is devoted to the discussion of an instructive analytic example. In section 5 we discuss extensively a transport problem in a cross-channel micromixer, that has been inspired by Balasuriya’s recent work [2]. There, we demonstrate both the conceptual and numerical benefit of a Lagrangian perspective on transport. We conclude with section 6, and discuss further implementation details in appendix B.

1.1 Continuum-mechanical terminology

Continuum mechanics considers two frameworks for the description of motion and deformation of continua (bodies/fluids) in space, which are commonly referred to as Eulerian and Lagrangian. Eulerian coordinates are assigned to spatial points in a fixed frame of reference, Lagrangian coordinates label material points and are often taken as the Eulerian coordinates at some initial time, say, . The configuration of the body at time is given by the flow map , which assigns material points their location in space, and its motion is given by the one-parameter family of configurations/flow maps .

In fluid dynamics, there are two important characteristic curves associated with the flow [3]. A path line through is the time-curve of a fixed Lagrangian particle in Eulerian coordinates, i.e., . In our context, we refer to the time-curve of a fixed Eulerian location in Lagrangian coordinates, i.e., , as the streak line through . In other words, the streak line is a collection of material points that will occupy the Eulerian position at some time. More commonly, streak lines are viewed as the collection of material points (identified with the current Eulerian locations) that have passed the Eulerian position at some time in the past. It can be imagined as an instantaneous curve of Lagrangian markers, injected in the past at and passively advected by the flow, see [3].

2 Lagrangian transport by compressible flows

Suppose is conserved by the flow generated by a time-dependent velocity field , i.e., solves


where is the initial distribution, with no-flux (zero Neumann) boundary conditions. Then the flux of across a stationary surface over a time interval is given by the well-known Eulerian flux integral



is the unit normal vector of the surface

indicating the direction of positive flux. This formula generalizes naturally to smoothly varying surfaces in space-time

where is the velocity field and is the unit normal vector field of in space-time (or, extended state space).

In the case of incompressible velocity fields that generate volume-preserving flows, Karrasch [12] has shown that the above Eulerian flux integral admits a Lagrangian equivalent. Our first main result is to show that the assumption of volume-preservation by the flow may be dropped without altering the conclusion. Thus, we solve Zhang’s donating region problem [21, Def. 1.2] completely.

Theorem 1.

For a given regular, time-dependent vector field , let be a conserved quantity satisfying the Eulerian scalar conservation law eq. 1. Let be a compact, connected, embedded codimension-one surface in (configuration space) , and be a compact time interval. Then there exists a decomposition , indexed by , of Lagrangian particles, identified by their spatial locations at time , that covers a set of full measure, such that


where is the unit normal vector field to characterizing the direction of positive flux.

Remark 1.

Equation 3 corrects the corresponding [12, Eq. (1)] by the factor . In [12] orientation issues have been ignored, such that the given equation there is correct only up to sign; however, it is fully correct for two spatial dimensions, i.e., .

Remark 2.

As in [12], we may generalize theorem 1 to smoothly moving surfaces , i.e., is everywhere transversal to time fibers . We omit the straightforward formulation and refer to [12, Problem 2] for the volume-preserving case.

Remark 3.

The decomposition is constructed at the initial time . Since it is material, it induces an equivalent decomposition at any other time instance , in particular at the final time instance.

We follow the strategy of [12], but provide a concise self-contained proof.



denote the map assigning to any crossing event on the extended surface the corresponding crossing particle , identified by the initial location of the latter. Here, is the flow map induced by taking particles from their spatial location at time to their spatial location at time . With the flow map at hand, it is well-known (see, e.g., [5, p. 11]) that the scalar density at a later time instance can be represented in terms of the initial distribution by


where denotes the derivative of at . In other words, the scalar density is anti-proportional (relative to the initial density) to the volume distortion by the flow. We show in appendix A that


Combining all together, we have


denotes the (topological) degree of [16] at a regular point , and the last equality holds by virtue of a direct consequence [10, Sect. 3.1.5, Thm. 6] of the area formula ([9, Thm. 3.2.3] and [14, Thm. 5.3.7]); cf. also [12].

Strictly speaking, the degree function may not be defined on the whole of : excluded are (i) the set of critical values of , i.e., in our context points whose trajectory has at least one non-transversal crossing of [12], and (ii) the image of the boundary of under . Both sets have measure zero: the first by Sard’s theorem, the latter as an -dimensional subset in [15, Ch. 6]. In other words, the contribution of these particles to transport is negligible. With the notation as in [12], , , we may decompose the last integral into

3 Application to two-dimensional flows

As indicated in [12], the use of the area formula is to account for multiple crossings of each particle from the donating region as well as the orientation of each crossing. In two-dimensional flows, bookkeeping of section222To avoid confusion, we will refer to as section (and not surface) in the following, since we will be concerned with two-dimensional applications. In this case, is a curve. crossings and their orientation can be performed easily. It turns out that the net number of crossings of a given particle , i.e., , coincides with its winding number relative to the closed curve [7, Sect. 6.6]; cf. [12, 22, 23]. Generally, does not need to be the topological boundary of the donating region , but we keep the notation by abuse of notation, and refer to it as the bounding curve in the sequel. The bounding curve is composed of the section, its back-advected image and the two streak lines emanating from the two endpoints of ; see fig. 1 for a visualization of the involved objects. The bounding curve may have a non-trivial topology due to self-intersections [12].

Figure 1: Visualization of the donating region and its bounding curve , corresponding to the example from section 4. The bounding curve is constructed as the concatenation of section, streak line, back-advected section and the other streak line. Positive flux is defined to be downward along (red). The induced orientation on (indicated by arrows) is such that points to the right. The coloring of the loops corresponds to their winding numbers, each connected piece corresponds to a simple loop (not all are labelled to avoid visual clutter). Enclosed white regions belong to , i.e., those particles have a vanishing number of net crossings and do not contribute to the total transport. The shown trajectory (dashed) originates from and has two (positive) crossings of .

Recall that the winding number of a point not on 333There exist different ways of defining a winding number for points lying on the curve (one of which is employed in [23]), but since the curve has measure/area zero this does not affect the integral equality, and there is no loss of generality by omitting this discussion. relative to is defined as the number of counter-clockwise turns that a vector connecting to does, as it follows the bounding curve in its positive direction/orientation. The bounding curve is equipped with the pullback orientation induced by [15, Prop. 15.15], i.e., it inherits the orientation of its pre-image ; for the latter we refer the reader to appendix A. In short, the orientation of is such that points to the right along in the plane of particles; cf. fig. 1.

As is well-known, the winding number is constant on the interior of simple loops of the donating region: any two particles that can be connected by a continuous curve without crossing the bounding curve admit the same winding number, which is why we may introduce the winding number of a simple loop relative to . Notably, the winding number is fully determined by the bounding curve alone.

With these considerations, we may re-write eq. 3 as


As has been mentioned in [12] and as we demonstrate later, eq. 6 is advantageous when transport by a material set of interest is to be computed. Since both the definition of and the Lagrangian flux integral are given in Lagrangian coordinates, the restriction of the integral to the set is direct and explicit. We summarize the above considerations for two-dimensional flow problems in algorithm 1. Some implementation details are given in appendix B.

Input: Section , velocity-field , time interval , boundary of the conditioning set of particles , distribution of at

  1. Calculate the bounding curve

    1. Calculate the streak lines through the endpoints of

    2. Calculate the back-advected section

    3. Assemble by concatenation of (back-advected) section and streak lines

  2. Divide bounding curve to obtain simple loops

    1. Find intersections of the bounding polygon

    2. Use intersections to assemble simple loops

  3. Calculate the winding numbers

    1. Find an interior point of each loop

    2. calculate the winding number of this point with respect to the bounding polygon

  4. Evaluate the r.h.s. of eq. 6, potentially restricted to

    1. Intersect each loop with

    2. Integrate the initial distribution of over the intersections, or calculate the area of the intersections (if )

    3. Sum up the integrals weighted by the corresponding winding numbers

Output: Accumulated flux due to material originating from

Algorithm 1 Lagrangian transport approach for two-dimensional flow problems

4 An analytic example

We want to test algorithm 1 on a simple example that still shows some nontrivial features like compressibility and higher winding numbers with mixed signs. To obtain higher winding numbers, we need particles that cross a given section multiple times in the same direction. It is therefore natural to consider a divergent rotational linear flow, similar to but simpler than the example in [12] and inspired by the tests performed in [20].

To this end, consider a superposition of rotation and strain , i.e.,


with the fluid density


That is we release a Gaussian fluid density , centered at the origin, at time . It is easily checked that solves the conservation law eq. 1.

We consider the transport problem for the section over the time interval , where positive flux is downward. The geometry of the donating region as computed by algorithm 1 is shown in fig. 1. Close to the origin, we see a region of winding number -1, which corresponds to particles which cross the section three times upward (i.e., negative crossings) and twice downward, with a net number of crossings equal to -1. The enclosed white regions correspond to particles that cross as often upward as downward, and therefore their transport contribution annihilates. The region with winding number 2 corresponds to particles that cross twice more often in the positive than in the negative direction. This is achieved by reflux from below that is bypassing the section ; see the shown trajectory (dashed).

Distance threshold Relative error Computation time [s]
0.05 3.18e-03 0.39
0.01 1.44e-04 0.56
0.005 4.87e-06 0.92
0.001 2.75e-07 4.16
Table 1: Analysis of accuracy and computation time (measured on one core of an Intel Core i5 2.3 GHz processor) for the Lagrangian transport computation, depending on a distance threshold used in the adaptive construction of streak lines and the back-advected section.

The relative error to the analytic result and the computation times for our implementation of algorithm 1 are given in table 1, depending on a threshold for the maximal distance of consecutive points along the streak lines and the back-advected section. This threshold is used in their adaptive calculation: the smaller the threshold, the finer resolved are streak lines and back-advected section; see appendix B for more details. A thorough numerical analysis of a similar algorithm has been performed by Zhang in [20].

5 Transport in a cross-channel micromixer

5.1 Set-up

We consider the following two-dimensional setup shown in fig. 2, which is strongly inspired by [2]; cf. also [19]. Two incompressible fluids enter a straight channel on the left side, fluid 1 in the lower region and fluid 2 in the upper region of the channel. In the absence of an external disturbance the two fluids remain separated by their common interface indicated by the horizontal line (green). Physically, the only possible mixing is due to diffusion along the interface, which is, however, not as efficient as often desired. To improve mixing for microfluidic applications, cross-channels are introduced, that are disturbing the horizontal flow in the vertical direction by periodic sucking and pumping. The resulting apparatus is called a cross-channel micromixer [19]. Consistently with [2], we neglect diffusive transport in the following, and discuss the validity of this assumption in section 6.

Figure 2: Setup of a cross-channel micromixer with 5 cross channels. The two fluids enter from the left and flow to the right. At the end of the channel, the unperturbed background velocity profile is shown. In the upper half, the full instantaneous velocity field at some instance in time is indicated. For further details, refer to the text.

We assume that the velocity field is modeled as the superposition of a steady horizontal background field and a time-dependent, -modulated vertical cross-channel velocity agitation, i.e.,

For the flow in the main channel, we assume a parabolic flow-pattern with flow speed at the interface and no-slip boundary conditions at the walls. The corresponding velocity profile is shown on the right of fig. 2. For the given geometric parameters, the steady background velocity field takes the form

where and are the widths of fluid 1 and fluid 2, respectively, see fig. 2. The assumed profile for the cross-channel flows is likewise parabolic and reads as


Here, is the maximal velocity magnitude at the center of the -th channel, is half of the width of the channel, is the position of the center of the channel, is the frequency of the vertical disturbance and represents the phase shift. The flow in the main channel is supposed to be dominant, i.e., we have and . An exemplary part of the superposed velocity field is shown in the upper half of the mixer in fig. 2 for some instant of time.

Due to the (assumed) absence of diffusive mixing, we will quantify mixing of the two fluids under the unsteady velocity agitation by the amount of fluid 1 that leaves cross-channel micromixer via the upper part of the channel (exit 2) which would be occupied by fluid 2 without the velocity agitation (and vice versa). Consequently, quantifying mixing is equivalent to quantifying transport, i.e., accumulated flux, of fluid 1 across the section , see fig. 2.

5.2 The Eulerian transport approach

5.2.1 Choice of section and other issues

In light of Eq. (2), it may seem that the question of quantifying transport of fluid 1 through exit 2 can be answered directly by applying the flux integral formula to some section like , that connects a point on the original fluid interface behind the cross-channel region with the upper channel wall; see fig. 2. This way, however, we would measure the out-flux of both fluid 1 and fluid 2 through exit 2, which is a useless assessment of mixing in the micromixer.

There are two possible attempts to fix this: (i) choose a different section to distinguish the two fluids according to their origin, and/or (ii) multiply the fluid density

by a characteristic function

corresponding to the region occupied by fluid 1 in space over time.

Regarding the first option, a natural candidate would be the initial interface itself, bounded by a point left of the cross-channel section and a point right of it, with normal direction ; that is, the portion of the black dashed line between and in fig. 3. The rationale behind such a choice is that any material carried by fluid 1 which leaves the micromixer through the upper exit 2 has to cross this section. Generally, this is a minimal requirement for a meaningful choice of section for this problem: essentially any section connecting to or to the upper wall has to be crossed by fluid 1 on the way to exit 2 and is therefore admissible. The issue here is, however, that with a direct application of the Eulerian flux integral to admissible sections we would again measure any flux of both fluid 1 and fluid 2. Also, it does not help to restrict attention to the positive flux only, i.e., replacing by , because this way we would neglect reflux of fluid 1 but would take into account reflux of fluid 2 that has invaded the lower region before.

In contrast, the second option would in principle do the job with any choice of section just discussed, but is numerically extremely challenging. While it is simple to decide at time which spatial region is occupied by fluid 1 and which one by fluid 2, this knowledge is expensive to obtain for later times as it corresponds to solving numerically the conservation law (1). Alternatively, one may proceed reversely: choose an admissible section and discretize to apply some numerical quadrature scheme to solve eq. 2. For each point , compute its past flow image at the initial time to see whether the particle occupying at time originated from fluid 1 or fluid 2. Correspondingly, its instantaneous flux value is, respectively, included and excluded from the numerical quadrature scheme. This kind of numerical quadrature turns out to be expensive even for short time intervals , see section 5.3.2.

5.2.2 The streak line approach

Recently, Balasuriya [2] proposed a more involved construction of a time-dependent section for a flow problem that is very related to the one described in section 5.1. He considers two anchor points and on the original fluid interface, which are located before and after the cross-channel region, respectively; see fig. 3.

Figure 3: Construction of the nominal fluid interface at a given time , following [2]: the upstream streak line through (red) and the vertical gate (cyan) at . The original interface for the undisturbed case is indicated by the horizontal thin dashed line. Particles located above the downstream streak line belong to fluid 2. The shown situation corresponds to an instantaneous flux of fluid 2 into the exit 1 region.

At any time instance , the so-called upstream streak line of (red) is considered, i.e., a collection of Eulerian positions of material points which have passed through before time . This corresponds to the classic notion of streak lines eluded to in section 1.1. The streak line is continued up to the point when it hits a line perpendicular to the original fluid interface and attached to . The nominal fluid interface is then constructed by a concatenation of upstream streak line and a straight, vertical gate (blue) connecting the hitting point with ; see [2]. Finally, one builds a continuously varying section in space-time by taking the union of all nominal fluid interfaces . This section can be considered as half-Eulerian and half-Lagrangian, since it is build partially from a material curve, the streak line, and partially from a spatially defined curve, the gate.

By construction, the streak line segment of the nominal fluid interface behaves like a material curve and, therefore, admits no flux. Therefore, the upstream streak line forms the fluid interface throughout, and depending on whether it hits the line above or below , some of fluid 1 is leaving through exit 2 or fluid 2 through exit 1. In the first case, the gate has an orientation such that flux to the right is measured positively, in the second case flux is measured negatively. By restricting to one sign of flux only, one may compute transport, say, of fluid 1 through exit 2.

Since the streak line segment of the nominal interface is material, any flux across is across the gate. The aim of [2] is to expand the (instantaneous) flux across the gate, i.e.,

in powers of the perturbation strength parameter (implicit in ) and to determine the leading order term444By construction, the nominal interface has only tangential variation over time, which is why the (normal component of the) relative velocity equals the original velocity .. This is a feasible approach provided that the streak line has a unique intersection with the line, which can be satisfied by choosing sufficiently small. In cases when the streak line has multiple intersections with a horizontal line (say, an S-shaped pattern), the gate would consist of several pieces, and the perturbative approach breaks down.

5.3 The Lagrangian transport approach

In this section, we demonstrate the efficacy of the Lagrangian transport formalism on the cross-channel micromixer with 5 cross-channels described in section 5.1. The parameters are given in table 2. For the experiments, we assume that a steady flow—with the cross-channels turned off—is established in the main channel before the experiment starts. Thus, fluid 1 occupies the lower, fluid 2 the upper region, and both are separated by their mutual interface at . Possible turbulence due to the sharp edges at the cross-channels is neglected. At time the cross-channels are switched on and the mixing process begins. For convenience, the initial distribution of the conserved quantity is set to throughout this section.

Parameter Variable Value(s)
Velocity at the interface of the channel 1
Height occupied by the lower fluid 1
Height occupied by the upper fluid 0.7
Positions of the centers of the cross-channels
Vertex-velocities of the cross-channels
Half-widths of cross-channels
Phase shift between cross-channels
Sucking/pumping frequency 4
Table 2: Parameters for the cross-channel mixer. The strength of the disturbance is chosen separately for each experiment.

5.3.1 Visual evaluation of donating region

First, we want to evaluate graphically how well the bounding curve captures the particles that are crossing the chosen section. To this end, we choose the section as the straight line along the initial interface from at the start of the first cross channel to at the end of the last cross channel. We set , which is considered as huge in the perturbative setting of [2].

Specifically, we compute the path lines for a grid of fluid 1 type particles over the time interval . If the end position of a particle has a positive -component, the particle has crossed and, therefore, contributes to the transport within the time interval. If it has a negative -component, it does not contribute to transport. In fig. 4, the particles are plotted in their initial positions and colored depending on whether they end up above (yellow) or below (gray) the initial interface. Independently, the first step of algorithm 1 is performed with 200 discretization points for and each. The resulting bounding curve of the donating regions is plotted in the same figure as the particles.

Figure 4: (a) Particles in their initial positions colored depending on whether they end up above (yellow) or below (gray) the initial interface after the time interval . The bounding curve is composed of section (black), streak lines (red) and back-advected section (dashed). (b) Close-up of fig. (a)a.

The results of the two independent experiments match perfectly, see fig. (a)a. All crossing sample particles (yellow) are enclosed correctly by . Note that the loops of the bounding curve that are located above the initial interface contain all the fluid-2 type particles that end up in the lower half. By increasing the particle resolution we again find perfect matching, see fig. (b)b.

An animation of this simulation experiment is provided in the Supplementary Material SM1.

5.3.2 Comparison of Eulerian and Lagrangian transport approaches

Next, we compare computational results and effort for evaluating the Lagrangian and the Eulerian versions of the flux integral. To this end, we set and aim at computing the transport of fluid 1 across the horizontal section described in section 5.3.1 over the time interval .

For the Eulerian flux integral an indicator function is used for the distinction of the two fluids, as detailed in section 5.2.1. Since , the term simply yields the -component of . The characteristic function is evaluated by back-advection of points of some grid discretization of back to the initial time , and subsequently checking if the -component of the result is less than 0. In this case, equals 1, otherwise it is 0. Finally, a numerical quadrature based on the trapezoidal rule was applied to compute the integral non-adaptively, for increasing numbers of integration points.

For the Lagrangian approach, and were discretized initially by 100 equidistant points each, and then refined adaptively, to yield 168 and 201 points, respectively; see appendix B for details on the adaptive integration scheme. The conditioned set of particles restricting the transport calculation to fluid 1 is specified by a sufficiently large rectangular polygon enclosing the lower half of the channel.

No. of integration points Relative error Integration time [s]
1k 1.4744e-01 0.5
5k 7.1727e-02 0.8
10k 4.9097e-02 1.2
25k 2.5313e-02 2.2
50k 1.7816e-02 3.4
75k 1.4843e-02 4.5
100k 1.1801e-02 5.5
250k 7.3630e-03 10.9
500k 4.9478e-03 20.6
750k 3.7236e-03 30.4
1,000k 3.3845e-03 37.9
2,000k 2.1410e-03 76.9
3,000k 1.7187e-03 102.4
5,000k 1.2459e-03 171.1
10,000k 7.3614e-04 364.5
30,000k 2.6525e-04 1036.8
Table 3: Analysis of accuracy and numerical effort for the conditional Eulerian transport computation, depending on the number of integration points.

The results are shown in table 3. The numerical quadrature of the Eulerian integral converges to the result of the Lagrangian calculation, the relative error decreases from approx. 15% (i.e., one correct digit) for 103=1k integration points to 0.026% (4 correct digits) for 3107 integration points. Thus, many sampling points are required to produce reasonably accurate results. The computation times555The timings are measured on an Intel Core i5 (dual core) with 2.3 GHz. increase from about 1 second for 103 up to more than 1,000 seconds for 3107 integration points. For comparison, the Lagrangian computation takes less than 2 seconds. In fact, in our specific example in the perturbed situation, transport can occur only in the range of the cross-channels. Therefore, it would be sufficient to apply numerical quadrature only to those parts of the section. This would roughly halve the computational effort; of course, such information is not available, in general.

5.3.3 Different choices of sections

Finally, we compare the transport calculation for different choices of sections to demonstrate that any section constructed according to the criteria discussed in section 5.2.1 measure essentially the same transport; see fig. 5 for the considered sections. The point is at the beginning of the first channel and is at the end of the last channel and therefore all sections are admissible. For all sections, the accumulated flux of fluid 1 is determined for different time intervals using the Lagrangian approach. The restriction to fluid 1 is achieved as in section 5.3.2.

Figure 5: Different choices of sections. The point is at the -position of the beginning of the first cross-channel, is at the -position of the ending of the last cross-channel.

We choose several values for from to in -steps. This implies that the time intervals become large which—for large —leads to possibly poor approximations of the donating region. To ensure proper computation, we set .

The resulting conditional accumulated fluxes are shown in fig. (a)a. As expected, the lines follow the same trend and an offset for the different sections is noticeable, which is, however, not constant over time. For each , the transport across section is maximal and the transport across is minimal. This is due to the fact that fluid 1 has to cross first anyway before getting to . Therefore, the difference between the two curves at is exactly the amount of fluid that has crossed but not yet by the time . The other offsets are explained analogously.

To get a finer intuition, consider the evolution of conditional fluxes across the sections, shown in fig. (b)b. As before, we are conditioning on flux of fluid 1 only. These curves are obtained from finite-differencing the conditional accumulated flux curves shown in fig. (a)a, and are therefore a short-term flux approximation to the actual instantaneous conditional flux densities. As expected, section is showing the strongest reflux, i.e., negative flux, of fluid 1 downward, since it is everywhere normal to the velocity agitation in the cross-channels. Up to small numerical inaccuracies around , when the boundaries of the donating region become extremely complicated, the flux across is always non-negative. This is trivial when taking an Eulerian viewpoint, since the velocity field is—at any time instance and at any point on —pointing rightwards (or vanishing at the wall), irrespective of whether fluid 1 or fluid 2 is crossing. From the Lagrangian viewpoint, recall that our computations are based on the area of induced donating regions.

We provide animations of the evolution of the donating regions over increasing time intervals for (SM2) and (SM3) for the parameters used in this section in the Supplementary Material.

Figure 6: Transport of fluid 1 across the sections from fig. 5. (a) Accumulated flux over the time intervals . (b) Instantaneous flux at , obtained from finite-differencing of the accumulated flux function in (a).

To see that the different sections indeed capture essentially the same transport, consider fig. 7. For all sections, the constituents of are shown for : The back-advected section (cyan), the two streak lines (brown and pink) and the section itself (in the color corresponding to fig. 5). The restricting region of origin is shown by the dotted rectangular area. In the Lagrangian approach, the transport of fluid 1 to exit 2 is calculated by integrating the initial distribution of over the simple loops of that have a non-empty intersection with the dotted area.

Figure 7: Components of for the sections shown in fig. 5 and for integration time : the sections (colored as in fig. 5), the two streak lines (pink and brown) and the back-advected section (cyan). The region of origin is shown by the dotted area. Note that for and one of the streak lines degenerates to a point, for in the respective upper endpoint of those sections.

Now consider point . As it is the endpoint of all four sections, its corresponding streak line (pink) is exactly the same. Thus, this part of the bounding curve yields the same contribution to the transport for all cases. The difference in the accumulated flux arises from the varying back-advected surfaces (cyan). The contribution to the transport by this part of the bounding polygon is the highest for and the lowest for which matches the results shown in fig. 6. So we can quantify the differences between the curves by calculating the contribution of the cyan part of the bounding polygon to the transport. As the shape of the back-advected section is not constant for different , this contribution is neither, which explains the non-constant differences between the conditional transport curves.

The above discussion applies verbatim to the time-dependent nominal interface introduced by Balasuriya in [2] and described in section 5.2.2. Here, at any time instance both endpoints are stationary and coincide with and , i.e., the ones of and . For any flow interval, the streak line parts of the boundary of the donating regions generated by , or are the same, and the only (time-dependent) differences are due to the final sections and their back-advected images (for and ) on the one hand, and the initial and final nominal interfaces on the other hand. Up to these explainable differences due to the exact choice of section, the same transport is measured by the Lagrangian transport approach and the Eulerian flux across the (possibly split) gate from the streak line approach [2].

6 Concluding remarks

We have extended the Lagrangian transport approach from the volume-preserving flow case [12] to the general case. Also, we demonstrated both the efficacy and accuracy of the Lagrangian approach to transport across surfaces in a cross-channel micromixer problem. As argued theoretically in [12], this approach is advantageous when the interest in the computation of transport (integrated flux) is restricted to a specific material set. Roughly speaking, instead of determining which points on the surface are occupied by material points of interest and subsequently computing the conditional Eulerian transport integral as usual, we transform the Eulerian transport integral into a Lagrangian one, whose restriction to the material set of interest is straightforward.

Consistently with [2], we have neglected diffusive transport in the micromixer problem. This assumption, however, appears to be rather unrealistic. As discussed in detail in [17] (and the references therein), diffusion is present, and may be even dominating advection as measured via the Péclet number. Nevertheless, it is often not efficient enough for the desired mixing purpose. This is why mixing by diffusion is further enhanced by the use of complicated channel geometries or cross-channel velocity agitations, or, generally speaking, by chaotic advection, [1]. The inclusion of the effect of diffusion turns the conservation law eq. 1 into an advection-diffusion equation, which is no longer of conservation/transport type. Still, a (generalized) conservation law can be formulated, to which our framework can be applied again, and work is in progress in this direction.

For another instructive comparison of the classic Eulerian with the Lagrangian perspective, consider the animations from the Supplementary Material, SM1 and SM3. In both cases, the exact same process as an evolution in time is shown: from the Eulerian perspective (SM1), in which particles are moving through space, colored according to whether they end up in the upper half at the end of the observation time interval; and from the Lagrangian perspective (SM3), where particles are colored according to whether they have crossed the horizontal section by the current time. Effectively, the fact that they have crossed the section is represented through coloring, not by drawing them in their current spatial position.

The Matlab implementation of algorithm 1, which is available on github666, is robust. In our many numerical experiments, only extremely thin filaments and near-tangent line segments appeared to be challenging for the winding number computation, which are generally known challenges in (time-dependent) polygon computations. Future work on the computational front will be concerned with implementations for three-dimensional flows.

Appendix A Proof of Equation 5

We prove

To this end, we need to first specify orientations; for an introduction to orientations, see, e.g., [15, Ch. 15].

For definiteness, we equip space and time with their respective standard orientations [15, Example 15.2]. Next, we equip with the product orientation, cf. [15, Prop. 15.7]: let be spatially (positively) oriented, then is oriented in space-time. For example, in two spatial dimensions time points upward as in the space-time figures in [12].

We equip the extended surface with an orientation induced by the unit normal vector field as follows [15, Prop. 15.21]. An ordered basis of the tangent space of is oriented if and only if is oriented in space-time. Finally, we equip the boundary of with the Stokes orientation [15, p. 386] by means of an outward-pointing vector field : an ordered basis of the boundary tangent space is oriented if and only if is oriented in , which holds if and only if is oriented in space-time.

For instance, when , at the “lower” boundary , the negative time vector is an outward-pointing vector, and the orientation along is given by such that is positively oriented in three-dimensional space-time. A simple right-hand check reveals that then points to the right of , or, equivalently, that is a positively oriented basis in space . Indeed, if is oriented in space-time, then so is , and by the above definition of the space-time orientation, it follows that is oriented in space.777Since acts on like the identity map, the induced orientation on is the one for which points to the right along .

Next, fix a regular point . Then acts diffeomorphically between an open neighborhood of in and its image in . We may choose local coordinates on such that in we have a positively oriented orthonormal basis, where , , and the (unit) positive time direction. The following parallelepipeds


span (positive) unit volume parallelepipeds, whose volume we also denote by the wedge product for notational simplicity.

On the one hand, we have


where the second equality holds true including the sign because of eq. 10. The determinant is an -determinant, where all vectors can be represented, for example, in the -basis. On the other hand, we first observe that

due to the multiplicativity of the determinant. The last determinant is again an -determinant, where the image vectors of can be represented in the -basis. Thus, we may compute as the change of volume under the composed action of on the oriented normalized parallelepiped (recall eq. 11):


where we have used (i) that acts like the identity map on tangent vectors , , of , and (ii) an identity for the streak vector field

which was proven in [12, Sec. 3]. To finish the proof, it remains to switch the sign of and to shift it to the first position in eq. 13, which requires reversions in orientation, hence the -coefficient.

Appendix B Implementation details

Step 1: Construction of

We use Matlab’s ode45 adaptive fourth-order Runge-Kutta ODE solver to compute streak lines directly from their definition and for the back-advected section. Additionally, we implemented an adaptive point insertion scheme to ensure that (i) points along the streak lines and the back-advected section are not too far apart and that (ii) consecutive line segments do not form an angle too far away from the straight one. At the same time, point insertion is prevented when consecutive points are already sufficiently close. This is to prevent self-overlapping zigzag-patterns in regions of high numerical sensitivity.

Step 2: Detection of intersections and of simple loops

Currently, self-intersections of the bounding polygon are computed by checking pairwise line segments for crossings. The intersection points are inserted into the polygon, and simple loops are extracted by going “as left as possible” at intersections in a counter-clockwise orientation, until one arrives at the starting point. A further decrease in computational effort could be achieved through the sweepline methodology [6].

Step 3: Construction of interior points and computation of winding numbers

Once simple loops have been computed, the next aim is to determine each individual winding number . To this end, we need to construct a representative interior point and compute its winding number classically by counting the number of turns of a line segment connecting to consecutive points along . For robustness, we construct many candidates for interior points, slightly on the left of each line segment of in a counter-clockwise parametrization. If all those winding numbers coincide, we take it as . Problems with this construction may occur in thin filaments, where the little left-offset may already yield an exterior point. In such a case, we construct the skeleton of the simple polygon and use its branching points as interior point candidates. As before, if their computed winding numbers all coincide we take it as .

Step 4: Integration of the Lagrangian flux formula

For uniform initial material densities and volume-preserving flows, the integration of over each simple loop reduces essentially to the computation of its respective area, for which Matlab admits the fast built-in function polyarea.

For general non-uniform initial material densities a numerical quadrature of has to be performed on each . There are at least two fundamental approaches to this. First, one may triangulate/mesh and then apply Gauss quadrature of some chosen degree on each triangle. There exist Matlab functions such as quadpts from the iFEM package [4], which give the location of the quadrature points in a triangle (in barycentric coordinates) and their corresponding weights. Second, there are theory and Matlab packages available that perform Gauss quadrature directly on polygons, PolyGauss [18]. For convex sets, the quadrature points are guaranteed to lie inside the domain of integration; for non-convex sets this may no longer be true, though, see [18].


We would like to thank Christian Ludwig for useful discussions and for support for his Matlab suite VISUALCOMPLEX. D.K. also acknowledges stimulating discussions with Jörg Schumacher.


  • [1] H. Aref. The development of chaotic advection. Phys. Fluids, 14(4):1315–1325, 2002. doi:10.1063/1.1458932.
  • [2] Sanjeeva Balasuriya. Transport between Two Fluids across Their Mutual Flow Interface: The Streakline Approach. SIAM Journal on Applied Dynamical Systems, 16(2):1015–1044, 2017. doi:10.1137/16M1089253.
  • [3] G. K. Batchelor. An Introduction to Fluid Dynamics. Cambridge Mathematical Library. Cambridge University Press, Cambridge, 2000. doi:10.1017/CBO9780511800955.
  • [4] L. Chen. iFEM: an integrated finite element methods package in MATLAB. 2008.
  • [5] A. J. Chorin and J. E. Marsden. A Mathematical Introduction to Fluid Mechanics, volume 4 of Texts in Applied Mathematics. Springer, 1993. doi:10.1007/978-1-4612-0883-9.
  • [6] M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars. Computational Geometry: Algorithms and Applications. Springer, 3 edition, 2008. doi:10.1007/978-3-540-77974-2.
  • [7] K. Deimling. Nonlinear Functional Analysis. Springer, 1985. doi:10.1007/978-3-662-00547-7.
  • [8] C. Dong, J. C. McWilliams, Y. Liu, and D. Chen. Global heat and salt transports by eddy movement. Nature Communications, 5(3294):1–6, 2014. doi:10.1038/ncomms4294.
  • [9] H. Federer. Geometric Measure Theory. Classics in Mathematics. Springer, 1996. doi:10.1007/978-3-642-62010-2.
  • [10] M. Giaquinta, G. Modica, and J. Souček. Cartesian Currents in the Calculus of Variations I: Cartesian Currents, volume 37 of Ergebnisse der Mathematik und ihrer Grenzgebiete. Springer, 1998.
  • [11] G. Haller. Langrangian Coherent Structures. Annu. Rev. Fluid Mech., 47(1):137–161, 2015. doi:10.1146/annurev-fluid-010313-141322.
  • [12] D. Karrasch. Lagrangian transport through surfaces in volume-preserving flows. SIAM J. Appl. Math., 76(3):1178–1190, 2016. doi:10.1137/15M1051348.
  • [13] D. Karrasch and J. Keller. A geometric heat-flow theory of Lagrangian coherent structures. 2016. arXiv:1608.05598.
  • [14] S. G. Krantz and H. Parks. Geometric Integration Theory. Cornerstones. Birkhäuser, 2008. doi:10.1007/978-0-8176-4679-0.
  • [15] J. M. Lee. Introduction to Smooth Manifolds, volume 218 of Graduate Texts in Mathematics. Springer, 2nd edition, 2012. doi:10.1007/978-1-4419-9982-5.
  • [16] J. W. Milnor. Topology from the Differentiable Viewpoint. The University Press of Virginia, Charlottesville, VA, 1965.
  • [17] N.-T. Nguyen and Z. Wu. Micromixers—a review. Journal of Micromechanics and Microengineering, 15(2):R1–R16, 2004. doi:10.1088/0960-1317/15/2/r01.
  • [18] A. Sommariva and M. Vianello. Product Gauss cubature over polygons based on Green’s integration formula. BIT Numer. Math., 47(2):441–453, 2007. doi:10.1007/s10543-007-0131-2.
  • [19] P. Tabeling, M. Chabert, A. Dodge, C. Jullien, and F. Okkels. Chaotic mixing in cross-channel micromixers. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 362(1818):987–1000, 2004. doi:10.1098/rsta.2003.1358.
  • [20] Q. Zhang. Highly accurate Lagrangian flux calculation via algebraic quadratures on spline-approximated donating regions. Comput. Methods Appl. Mech. Engrg., 264:191–204, 2013. doi:10.1016/j.cma.2013.05.024.
  • [21] Q. Zhang. On a Family of Unsplit Advection Algorithms for Volume-of-Fluid Methods. SIAM J. Numer. Anal., 51(5):2822–2850, 2013. doi:10.1137/120897882.
  • [22] Q. Zhang. On Donating Regions: Lagrangian Flux through a Fixed Curve. SIAM Rev., 55(3):443–461, 2013. doi:10.1137/100796406.
  • [23] Q. Zhang.

    On generalized donating regions: Classifying Lagrangian fluxing particles through a fixed curve in the plane.

    J. Math. Anal. Appl., 424(2):861–877, 2015. doi:10.1016/j.jmaa.2014.11.043.
  • [24] Z. Zhang, W. Wang, and B. Qiu. Oceanic mass transport by mesoscale eddies. Science, 345(6194):322–324, 2014. doi:10.1126/science.1252418.