An exact solution to dispersion of a passive scalar by a periodic shear flow
UUnder consideration for publication in J. Fluid Mech. An exact solution to dispersion of a passivescalar by a periodic shear flow
Miguel A. Jimenez-Urias † and Thomas W. N. Haine Earth and Planetary Sciences, Johns Hopkins University, Baltimore, MD 21218, USA(Received xx; revised xx; accepted xx)
We present an exact analytical solution to the problem of shear dispersion given ageneral initial condition. The solution is expressed as an infinite series expansion involvingMathieu functions and their eigenvalues. The eigenvalue system depends on the imaginaryparameter q = 2 ik Pe , with k the wavenumber that determines the tracer scale in theinitial condition and Pe the P´eclet number. Solutions are valid for all Pe , t >
0, and k > q = q EP(cid:96) called Exceptional Points (EPs), the first occurringat q EP ≈ . i . For values of q (cid:47) . i , all the eigenvalues are real, different andeigenfunctions decay with time, thus shear dispersion can be represented as a diffusiveprocess. For values of q (cid:39) . i , pairs of eigenvalues coalesce at EPs becoming complexconjugates, the eigenfunctions propagate and decay with time, and so shear dispersionis no longer a purely diffusive process.The limit q → k > k (cid:28) / Pe . The latter implies k →
0, strong separation of scales between the tracer and flow. The limit q → ∞ resultsfrom large P´eclet number for any k >
0, or from large k and non-vanishing Pe .We derive an exact closure that is continuous in wavenumber space. At small q , theclosure approaches a diffusion operator with an effective diffusivity proportional to U /κ ,for flow speed U and diffusivity κ . At large q , the closure approaches the sum of anadvection operator plus a half-derivative operator (differential operator of fractionalorder), the latter with coefficient proportional to √ κU .
1. Introduction
Shear dispersion describes the enhanced spreading of a conservative passive scalartracer due to the interaction of velocity gradients with weak diffusivities (Taylor 1953;Batchelor 1956; Aris 1956; Young & Jones 1991). Approximate models of shear dispersionsuggest that the problem can be split into two stages (Rhines & Young 1983). First, thereis a rapid stage during which stirring by the spatially-varying flow replaces the initialcondition by a generalized average along a streamline. Second, there is a slower stage during which diffusive mixing takes place with an effective diffusivity κ ∗ ∝ U /κ , where U is the characteristic flow speed and κ is the molecular diffusivity. The enhancementof diffusivity by shear dispersion, which is strongest for small κ , is referred to as G. I.Taylor’s effective diffusivity κ ∗ (Taylor 1953).Shear dispersion is often regarded as a special case of the method of homogenization(Haynes & Vanneste 2014; Hinch 1991 Ch. 7; Holmes 2012, Ch. 5). Homogenization ofperiodic flows exploits the separation of scales between the (small) scale of the flow andthe (large) scale of the tracer, thus it is restricted to timescales longer than the diffusivetimescale and so it fails to accurately represent the early-time behavior. Despite the lack † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] J a n M. A. Jimenez-Urias and T. W. N. Haine of a solution for all times t >
0, the effective diffusivity κ ∗ is generally regarded as validin the large P´eclet number limit (although claims have been made that it applies onlyfor small P´eclet number; see Young & Jones 1991 for a brief discussion). In the largeP´eclet number limit, approximate asymptotic solutions have been derived to extend theeffective diffusivity, which involve higher order corrections that attempt to model early-time behavior (Mercer & Roberts 1990; Young & Jones 1991; Haynes & Vanneste 2014).All these approximate solutions rely on homogenization theory, however, and thus fail todescribe the complete spatial evolution of the tracer for all t > et al. et al. et al. et al. et al. et al. t > Pe and wavenumbers k (except at combinations of k and Pe that result inspecific values of their product). The solution is obtained in accordance to Floquet theory,and it depends on a single parameter q , which is a function of the P´eclet number andthe wavenumbers of the initial condition. With the exact solution, we derive in section3 the exact closure to the advective term that is continuous across wavenumbers and wederive approximations to the exact closure based on its dependence on q . These closuresclarify the limits of validity of Taylor’s effective diffusivity and reveal a novel regime. Adiscussion of implications and future opportunities is in section 4 and conclusions are insection 5.
2. Plane Shear Dispersion
Problem Statement
Consider the evolution of an inactive passive tracer θ = θ ( x, y, t ) in a re-entrantchannel of width M . There is a steady, unidirectional, non-divergent horizontal velocity u = ( U ( y ) , U ( y ) = U cos (2 πy/M ), and molecular diffusion with diffusivitycoefficient κ . The advection diffusion equation determines the evolution of the tracer andreads ∂θ∂t + U cos (cid:18) πyM (cid:19) ∂θ∂x = κ ∇ θ, (2.1) θ ( x, y,
0) = θ cos (cid:18) π kxL (cid:19) Φ (cid:18) π yM (cid:19) , (2.2) hear Dispersion Figure 1.
Initial condition of a tracer θ ( x ∗ , y ∗ ,
0) in a channel that is periodic in x ∗ , and withslippery walls (shaded gray) at y ∗ = { , π } . All variables are non-dimensional. The thick blackline and arrows depict the shear flow defined by a single Fourier mode in y ∗ , with wavenumber (cid:96) = 1. The initial tracer is uniform in the y -direction, so that Φ ( y ∗ ) = 1. The reentrant channelspans 10 π < x ∗ < π , and so k (cid:62) . x ∗ for a dimensional length of L = 10. A larger choice of L increases the range of possible k -values. with periodic boundary condition in x and no-flux boundary conditions at y = { , M } .The initial condition (2.2) sets the scale of the tracer via the wavenumber k > x -length of the channel ( L ) along with the condition ofperiodicity in x , set the smallest possible value of k > L = 1 / min ( k )(see Fig. 1). Nonetheless, due to the re entrant boundary condition, the value of L isphysically irrelevant as it can be made a large as desired in order to permit any value of k > Φ in (2.2) defines the cross-channel ( y -) structure of the initial condition,with the added restriction for Φ to be square-integrable so that it can be expressedas a Fourier series. Thus, more complicated ( x, y )-dependent initial conditions can beconstructed as linear combinations of Fourier modes each with different values of k .We non-dimensionlize equation (2.1) and the initial condition according to the scaling y = M y ∗ π , x = Lx ∗ π , t = M t ∗ π κ , (2.3)resulting in the non-dimensional advection-diffusion equation (cid:15) ∂θ ∗ ∂t ∗ + cos( y ∗ ) ∂θ ∗ ∂x ∗ − Pe ∂ θ∂x ∗ − (cid:15) ∂ θ∂y ∗ = 0 , (2.4)with initial condition θ ( x ∗ , y ∗ ,
0) = θ cos( kx ∗ ) Φ ( y ∗ ). The non-dimensional domain isthen defined as ( − π (cid:54) kx ∗ (cid:54) π ) × (0 (cid:54) y ∗ (cid:54) π ), and t ∗ >
0. There are two non-dimensional parameters in (2.4): Pe = U L π κ , (cid:15) = t a t d = L M Pe , (2.5)where Pe is the P´eclet number, a measure of the strength of diffusive fluxes compared toadvective transport, and (cid:15) is the ratio between the advective and diffusive time scales, M. A. Jimenez-Urias and T. W. N. Haine given by t a = L/ (2 π U ) and t d = M / (4 π κ ), respectively. The ratio M/L is referred toas a scale-separation parameter and when considered small it is used to derive asymptoticsolutions to (2.4) (see Camassa et al. a priori assumption of scaleseparation and, without loss of generality we set L = M , which implies that the onlyphysically relevant length scale is given by the width of the varying flow. Then (cid:15) = Pe − and thus P´eclet number is the only (free) parameter that arises from the equation (2.4).The other independent parameter is the wavenumber of the initial condition k > k (cid:28) (cid:96) ≡ U ( y ) = cos( (cid:96)y ∗ )).Thus, Taylor’s effective diffusivity is valid in the limit Pe → ∞ , k (cid:28)
1, and t ∗ > Problem Solution
We seek x -periodic solutions † to the homogeneous problem (2.4) using the ansatz θ ( x, y, t ) = (cid:60) { A ( k, ω, y ) exp [ ikx − ωt ] } , (2.6)where the frequency ω = ω ( k ) is a function of k via the dispersion relation of (2.4). ‡ Ingeneral, ω may be complex, with the real part indicating decaying/growing solutions andthe imaginary part indicating propagation. This separable solution is motivated by thechoice of separable initial condition (2.2). Substituting (2.6) into (2.4) yields the equationfor amplitude A ( k, ω, y ) d Ady + (cid:2) ω − k − ik Pe cos( y ) (cid:3) A = 0 , (2.7)which is, by construction, independent of time. Re-scaling the domain as ˜ y = y/
2, anddefining the parameters a = 4 ω − k , q = 2 ik Pe , (2.8)converts (2.7) into the canonical Mathieu equation d Ad ˜ y + [ a − q cos(2˜ y )] A = 0 . (2.9)In accordace to Floquet theory, solutions to (2.9) that are π -periodic in ˜ y (2 π -periodic in y ) and satisfy the no-flux boundary conditions, are the even, cosine-elliptic Mathieufunctions ce n ( q, ˜ y ), n = 0 , , , · · · , also called Mathieu functions of the first kind(McLachlan 1947; Olver et al. n isthe eigenvalue a = a n ( q ) (also called the characteristic number). Both ce n and a n depend on the parameter q = 2 ik Pe (McLachlan 1947; Olver et al. q , which contains all therelevant physics. Because Mathieu functions converge to Fourier series as q vanishes,lim q → ce n ( q, ˜ y ) = cos(2 n ˜ y ), q can be considered an eccentricity parameter for eachcosine-elliptic function.From (2.8), the dispersion relation of (2.4) is ω n = a n ( q )4 + k , (2.10) † Note that the stars have been dropped, so henceforth all variables are non-dimensional. ‡ For example, in the 1D heat equation, which is the limit U → , k → ω = k (Deconinck et al. hear Dispersion l q EPl
Eigenvalue pairs0 1 . i a , a . i a , a . i a , a . i a , a Table 1.
Locations q EPl of the first four exceptional points (EPs) for Mathieu’s cosine-ellipticeigenfunctions (from Blanch & Clemm 1969). At EPs two successive eigenvalues coalesce(branch) in their real (imaginary) parts with increasing q . See Fig. 2. which sets the time dependence of each eigenfunction ce n . Specifically, from (2.6) thetime dependence of each ce n in A ( k, ω, y ) is given by the sum of k (1D diffusion)and the eigenvalue a n (see Fig. 2). For large enough k , ω n is dominated by the k term and the ce n mode decays exponentially with timescale τ ∼ / (cid:60){ ω n } ≈ /k .Otherwise, for small k , or large Pe (or both), ω n is dominated by the a n term and1D diffusion is negligible. Therefore, the functional dependence of a n with q (and hencewith k = q/ (2 i Pe )) is key to understanding the solution.The eigenvalues a n ( q ) exhibit a transition with increasing q when the parameter q is purely imaginary (as it is here). For q (cid:47) . i , the a n are all real and positive.Therefore, the ce n modes all decay exponentially with time. As q grows, successivepairs of eigenvalues merge into complex conjugate pairs (with positive real parts). Thesemergers happen at Exceptional Points (EPs), as seen in Fig. 2 (Mulholland & Goldstein1929; Hunter & Guerrieri 1981; Olver et al. et al. n occur at well known and fixed values of q (Table 1), and the q values of the EPs increase with n . The gravest eigenvalues a , a coalesce (in their realparts) and branch (in their imaginary parts) at q = q EP ≈ . i (Fig. 2, Table 1). Thenext pair of eigenvalues merges at q EP ≈ . i , and so on. Transition of a n past EPsmakes the eigenfrequency ω n for that mode complex. This indicates a different physicalbehavior: the ce n mode now propagates along the channel (and decays). The importanceof EPs has recently become evident in several areas of mathematical physics (Bender &Boettcher 1998; Bender 1999; Heiss 2004, 2012; Miri & Alu 2019). Shear dispersion isanother example.Mathieu functions ce n ( q, ˜ y ) can be defined by a Fourier series (McLachlan 1947;Arscott 2014). For ce n ( q, ˜ y ), this isce n ( q, ˜ y ) = ∞ (cid:88) r =0 A (2 n )2 r ( q ) cos(2 r ˜ y ) , n = 0 , , , · · · (2.11)For each n and q , the eigenvalue a n ( q ) and the sequence of Fourier coefficients A (2 n )2 r ( q )are calculated by standard procedures (see Appendix A and the github software reposi-tory for this paper).As q increases past each EP and the two eigenvalues coalesce, the associated eigen-functions also coalesce to become complex conjugates and for values q > q EPl , theseeigenfunctions satisfy a shifted-conjugate symmetry, e.g. ce ( q, ˜ y ) = ce ∗ ( q, ˜ y − π / ∗ denotes the complex conjugate (Fig. 3). Except at EPs, Mathieu eigenfunctionsfurnish a complete base, equiconvergent with Fourier series (see Olver et al. y -dependent function in the Hilbertspace of periodic (square integrable) functions. † † Exactly at the EP, the eigenvalues are repeated and the eigenfunctions no longer span
M. A. Jimenez-Urias and T. W. N. Haine
Figure 2.
Comparison of the two terms that define the eigenfrequency w n in (2.10) associatedwith the gravest four Mathieu eigenfunctions n = 0 , , ,
3, and their dependence on wavenumber k . The wavenumbers of the exceptional points (EPs) are defined by | q EPl | / Pe where the q EPl values are fixed (Table 1). The EPs shift towards low wavenumbers when Pe > Pe <
1. Colored lines represent Pe = 1, with black starred lines Pe = 10, and gray lines Pe = 100. The shaded region, delimited by the curve k , shows therange of wavenumbers for which diffusion dominates ω n , which increases the decay rate of theeigenmodes. Diffusion has no effect on the imaginary part of the eigenvalues. From the general definition (2.11), the exact solution to (2.4) is given by a linearcombination of even, cosine-elliptic Mathieu functions. This is θ ( x, ˜ y, t ) = (cid:60) (cid:40) ∞ (cid:88) n =0 α n ( q ) ce n ( q, ˜ y ) exp (cid:104) ikx − (cid:16) a n k (cid:17) t (cid:105)(cid:41) . (2.12)The coefficients α n are those requiring (2.12) to satisfy the initial condition ∞ (cid:88) n =0 α n ( q ) ce n ( q, ˜ y ) = Φ (˜ y ) . (2.13)Since Mathieu eigenfunctions ce n ( q, ˜ y ) are expressed as a Fourier cosine series, andassuming that Φ (˜ y ) has a Fourier cosine series, then the individual α n coefficients areconstructed from known Mathieu function identities (Olver et al. Φ = 1, we have the identity2 ∞ (cid:88) n =0 A (2 n )0 ( q ) ce n ( q, ˜ y ) = 1 , (2.14) the complete space. Nonetheless, the set of eigenfunctions can be made complete by defining ageneralized eigenfunction (Brimacombe et al. hear Dispersion Figure 3.
First four Mathieu eigenfunctions ce n ( q, y ) over a range of the parameter0 < q = 2 ik Pe < i , which extends past their EPs (Table 1). Colored lines arece n ( q = 900 i, y ) and dark gray lines are ce n ( q = 0 , y ). The ce n evaluated at q = 900 i (colored lines), can equally represent a solution when k = 10 and Pe = 45, or when k = 0 . Pe = 4500. Note that the abscissa is y = 2˜ y . where A (2 n )0 is the zeroth Fourier coefficient in the definition of each Mathieu functionce n in (2.11). Thus from (2.13) α n = 2 A (2 n )0 , and the exact solution to (2.4) with theinitial condition θ ( x, y,
0) = cos ( kx ) is θ ( x, ˜ y, t ) = (cid:60) (cid:40) ∞ (cid:88) n =0 A (2 n )0 ( q ) ce n ( q, ˜ y ) exp (cid:104) ikx − (cid:16) a n k (cid:17) t (cid:105)(cid:41) . (2.15)Similarly, for Φ (˜ y ) = cos(2 p ˜ y ), we have the following identity ∞ (cid:88) n =0 A (2 n )2 p ce n ( q, ˜ y ) = cos(2 p ˜ y ) , r (cid:54) = 0 , (2.16)where A (2 n )2 p is the p th Fourier coefficient in the definition of Mathieu functions in (2.11).Hence, α n = A (2 n )2 p and the exact solution to (2.4) with initial condition θ ( x, y,
0) =
M. A. Jimenez-Urias and T. W. N. Haine
Figure 4.
Time evolution of a passive tracer concentrated at x = 0 defined as a Gaussian in x ,and Φ (˜ y ) = 1 with Pe = 2000, in the presence of a periodic shear flow ( U (˜ y ) ,
0) = (cos(2˜ y ) , < k < .
25, possible by increasing the domain in x as − π (cid:54) x (cid:54) π . In this initial condition, every mode has a scale much larger than thatof the velocity field. cos( kx ) cos(2 p ˜ y ) is θ ( x, ˜ y, t ) = (cid:60) (cid:40) ∞ (cid:88) n =0 A (2 n )2 p ( q ) ce n ( q, ˜ y ) exp (cid:104) ikx − (cid:16) a n k (cid:17) t (cid:105)(cid:41) . (2.17)In general, let Φ (˜ y ) have a Fourier cosine series of the form Φ (˜ y ) = ∞ (cid:88) p =0 β p cos (2 p ˜ y ) . (2.18)Then, the general solution to the shear dispersion problem (2.1) with initial condition θ ( x, ˜ y,
0) = cos ( kx ) Φ (˜ y ) is θ ( x, ˜ y, t ) = (cid:60) (cid:40) ∞ (cid:88) n =0 ∞ (cid:88) p =0 (1 + δ p ) β p A (2 n )2 p ( q ) ce n ( q, ˜ y ) exp (cid:104) ikx − (cid:16) a n k (cid:17) t (cid:105)(cid:41) , (2.19)where δ p is the Kronecker delta † . This solution describes the initial condition consistingof a single Fourier mode in x . More general initial conditions can be defined as a linearcombination (in x ) of terms like (2.19) for different values of k . In other words, an initialcondition θ ( x, ˜ y,
0) = f ( x, ˜ y ) is handled by expressing f ( x, ˜ y ) as a double Fourier cosineseries in x and y and a triple sum that generalizes (2.19).2.3. Example Solution: a Gaussian initial condition
Although both the purely decaying ( q →
0) and propagating regimes ( q → ∞ ) admitsolutions associated with the large P´eclet number limit, only the propagating regime candescribe solutions associated with a wide range of wavenumbers in the initial condition.This is well illustrated by a Gaussian initial condition (Fig. 4). The initial conditionis given by a finite sum of cosine-modes in x , and the solution is a linear combinationof cosine solutions each of the form (2.15), with all modes satisfying k < .
25. Witha value of Pe = 2000, the range in the Mathieu parameter for each cosine-mode is20 i < q < i and so the gravest mode lies beyond the first two EPs (for Pe = 2000only the wavenumbers k < . q (cid:47) . i ). Similarly, for the longest mode k = 0 .
25 to lie before the first EP q EP ≈ . i , the P´eclet number must be Pe < . † The general solution implies α n ( q ) = (cid:80) ∞ p =0 (1 + δ p ) β p A (2 n )2 p ( q ) in (2.12) and (2.13) hear Dispersion q and large P´eclet number captures both the early and late timebehavior of the exact solution.
3. Closure for Average Shear Dispersion
The exact analytical solution to the advection diffusion equation, allows us to explorethe exact representation of the averaged equation and thus derive an exact closureequation that is uniform across scales. From the exact closure equation, approximationsare derived with particular applicability and utility.3.1.
Closure Statement
Consider the family of weighted averages defined by (cid:104) f (˜ y ) (cid:105) m = 1 π (cid:90) π f (˜ y ) cos(2 m ˜ y ) d ˜ y, m = 0 , , , · · · , (3.1)where cos(2 m ˜ y ) represents the weight/filter. Applying these operators to an initialcondition whose ˜ y -dependence is a Fourier cosine series, as in (2.18), gives (cid:104) Φ (cid:105) m = β m (cid:18) δ m (cid:19) (3.2)Thus, the average (3.1) selects the p = m cosine mode in the general initial condition.The average of Mathieu functions (2.11) yields (cid:104) ce n ( q, ˜ y ) (cid:105) m = A (2 n )2 m (cid:18) δ m (cid:19) (3.3)Similarly, the average (cid:104) ce n (cid:105) m selects the single r = m cosine mode in the infinite sumthat defines each eigenfunction from (2.11). Note that both in (3.2) and (3.3) the averageswith m (cid:54) = 0 have a factor of 1 / (cid:104) θ (cid:105) m = (cid:60) (cid:40) ∞ (cid:88) n =0 ∞ (cid:88) p =0 (1 + δ p ) β p A (2 n )2 p A (2 n )2 m ((1 + δ m ) /
2) exp (cid:104) ikx − (cid:16) a n k (cid:17) t (cid:105)(cid:41) . (3.4)The initial condition (cid:104) θ ( x, y, (cid:105) m = (cid:104) cos( kx ) Φ ( y ) (cid:105) m = cos( kx ) β m ((1 + δ m ) /
2) isguaranteed for all choices of m as long as Φ is defined as in (2.18) due to the orthogonalityof the A (2 n )2 p coefficients that defined Mathieu functions (see Appendix A and (A 3)).Applying the average (3.1) to the advection-diffusion equation (2.4) yields ∂ (cid:104) θ (cid:105) m ∂t + Pe ∂ ( (cid:104) θ cos(2˜ y ) (cid:105) m ) ∂x = ∂ (cid:104) θ (cid:105) m ∂x − m (cid:104) θ (cid:105) m (3.5)for all m = 0 , , , · · · . Interest lies in the closure of the advection term, namely, thesecond term on the left hand side of (3.5). We aim to write it exactly as a differentialoperator acting on the mean tracer field and to approximate it, so as to render a simplercompact form.0 M. A. Jimenez-Urias and T. W. N. Haine
Closure Solution
The expression of interest is ∂ ( (cid:104) θ cos(2˜ y ) (cid:105) m ) ∂x ≡ Λ m [ (cid:104) θ (cid:105) m ] , (3.6)where Λ m is the closure operator acting on the mean tracer (cid:104) θ (cid:105) m . The left hand sideof (3.6) is given exactly by (cid:60) (cid:40) ik exp ( ikx ) ∞ (cid:88) n =0 ∞ (cid:88) p =0 (1 + δ p ) A (2 n )2 p (cid:104) ce n cos (2˜ y ) (cid:105) m exp (cid:104) − (cid:16) a n k (cid:17) t (cid:105)(cid:41) (3.7)Since the exact averaged solution (3.4) can be written as the infinite sum (cid:104) θ (cid:105) m = (cid:60) (cid:8)(cid:80) ∞ n =0 (cid:104) θ (2 n ) (cid:105) m (cid:9) , then (3.7) implies that the closure equation (3.6) is written foreach (cid:104) θ (2 n ) (cid:105) m as follows ik (cid:104) ce n ( q, ˜ y ) cos(2˜ y ) (cid:105) m = Λ (2 n )2 m A (2 n )2 m (cid:18) δ m (cid:19) (3.8)Thus the averaged, unclosed advection-diffusion equation (3.5) can be written as anautonomous, closed system of equations for each n -mode, which in Fourier space readslike ∂ (cid:104) θ (2 n ) (cid:105) m ∂t = − (cid:16) k + 4 m + Pe Λ (2 n )2 m (cid:17) (cid:104) θ (2 n ) (cid:105) m . (3.9)The n -modes are subject to the initial condition (cid:60) (cid:40) ∞ (cid:88) n =0 (cid:104) θ (2 n ) (cid:105) m ( x, t = 0) (cid:41) = β m (cid:18) δ m (cid:19) cos( kx ) . (3.10)We thus proceed by first computing the left hand side of (3.8), which contains all the ˜ y dependence in θ . This is2 (cid:104) cos(2˜ y )ce n ( q, ˜ y ) (cid:105) m = 2 ∞ (cid:88) r =0 A (2 n )2 r (cid:104) cos(2˜ y ) cos(2 r ˜ y ) (cid:105) m , = ∞ (cid:88) r =0 A (2 n )2 r (cid:104) cos (2 [ r −
1] ˜ y ) + cos (2 [ r + 1] ˜ y ) (cid:105) m , (3.11)where the Mathieu function definition (2.11) is used. For m = 0 in the weighted average(3.1) (arithmetic cross-channel mean), the only term in the infinite sum that survives isthe r = 1 Fourier coefficient, for each n = 0 , , , · · · . Thus, when m = 0, (cid:104) cos(2˜ y )ce n ( q, ˜ y ) (cid:105) = A (2 n )2 . (3.12)For m >
0, the (cid:104)·(cid:105) m in (3.11) involves the two cosine terms multiplying the weightcos(2 m ˜ y ). They are rewritten as2 cos (2 [ r ±
1] ˜ y ) cos(2 m ˜ y ) = cos (2 [ r + m ±
1] ˜ y ) + cos (2 [ r − m ±
1] ˜ y ) , (3.13)again using trigonometric identities. The non-vanishing terms in the average are thosefor which the arguments of either cosine on the right hand side of (3.13) vanish. Considerfirst m = 1 after which we generalize to all m (cid:62)
2. When m = 1, only r = 0 contributesto (3.13) with the plus sign, whereas both r = 0 and r = 2 contribute with the minus hear Dispersion m = 1 (cid:104) cos(2˜ y )ce n ( q, ˜ y ) (cid:105) = 2 A (2 n )0 + A (2 n )4 . (3.14)When m (cid:62)
2, only r = m − r = m + 1contributes with the minus sign. Therefore, when m (cid:62) (cid:104) cos(2˜ y ) ce n ( q, ˜ y ) (cid:105) m = A (2 n )2( m − + A (2 n )2( m +1) . (3.15)With these, we proceed to calculate the exact closure operator in (3.8) for each choiceof m and n as follows ikA (2 n )2 = 2 Λ (2 n )0 A (2 n )0 m = 0 , (3.16) ik (cid:104) A (2 n )0 + A (2 n )4 (cid:105) = 2 Λ (2 n )2 A (2 n )2 m = 1 , (3.17) ik (cid:104) A (2 n )2( m − + A (2 n )2( m +1) (cid:105) = 2 Λ (2 n )2 m A (2 n )2 m m (cid:62) . (3.18)Or, written as a recurrence for mΛ (2 n )0 A (2 n )0 − (cid:18) ik (cid:19) A (2 n )2 = 0 m = 0 ,Λ (2 n )2 A (2 n )2 − (cid:18) ik (cid:19) (cid:104) A (2 n )0 + A (2 n )4 (cid:105) = 0 m = 1 ,Λ (2 n )2 m A (2 n )2 m − (cid:18) ik (cid:19) (cid:104) A (2 n )2( m − + A (2 n )2( m +1) (cid:105) = 0 m (cid:62) . (3.19)Remarkably, this recurrence relation matches term by term the recurrence associatedwith Mathieu’s eigenvalue system (see Appendix A, (A 1)). Hence, for each choice ofmode m in the averaging filter (3.1), we arrive at a family of exact closures for eacheigenfunction ce n . Each closure operator is proportional to the eigenvalue a n of theassociated eigenfunction: Λ (2 n )2 m = a n − m Pe (3.20)The closure above applies for all m = 0 , , , · · · .It is instructive to write the closed equation (3.9), say for m = 0, in physical space andin terms of dimensional variables. From (2.3) and (2.5) the closed dimensional equationis ∂ (cid:104) θ (2 n ) (cid:105) ∂t = κ ∂ (cid:104) θ (2 n ) (cid:105) ∂x − π κ Pe M Λ (2 n )0 (cid:104) (cid:104) θ (2 n ) (cid:105) (cid:105) , (3.21)where Λ (2 n )0 is an operator acting on the averaged mode (cid:104) θ (2 n ) (cid:105) .3.3. Asymptotic Closure for Small and Large Canonical Parameters
In (3.20) we found that the exact closure for each eigenfunction is proportional to theeigenvalue a n ( q ), which is a function of the canonical parameter q = 2 ik Pe . To expressthe closures analytically, and to better understand the closures across the two regimes ofthe solution, consider the asymptotic limits of small and large q .For q → a n and therefore the closures Λ (2 n )2 m can be approximatedasymptotically as an even power series in q (see Appendix B). The first term in theexpansion of the closure operator for the first four eigenfunctions for the cross-channel2 M. A. Jimenez-Urias and T. W. N. Haine arithmetic mean m = 0 reads Λ (0)0 ∼ k Pe + O ( k Pe ) ,Λ (2)0 ∼ Pe − k Pe + O ( k Pe ) ,Λ (4)0 ∼ Pe − k Pe + O ( k Pe ) ,Λ (6)0 ∼ Pe − k Pe + O ( k Pe ) . (3.22)The presence of the k Pe terms suggests a second derivative that enhances the along-stream diffusion in the original advection-diffusion equation (2.4). This follows by con-necting Fourier space and physical space with ik ←→ ∂ x . Thus, in the q → Λ (0)0 , written dimensionally following (3.21), enhances the molecular diffusion as κ (cid:18) Pe (cid:19) ∂ (cid:104) θ (0) (cid:105) ∂x . (3.23)Using (2.5) this term is κ Pe M U π κ , (3.24)which is Taylor’s effective diffusivity (Taylor 1953). This approximation is valid in thelimit | q | = 2 k Pe →
0, which is attained either by large Pe under the scale-separationcondition 2 k (cid:28) / Pe , or by small P´eclet number for finite k >
0. This latter limit ischaracterized by strongly diffusive, weakly advective flows.The closure in (3.23) is for the gravest, slowest-decaying mode ( n = 0). A similaroperator applies to the other modes n >
0, with a term proportional to Pe . Higher orderterms in the expansion (see Appendix B, equation (B 1)) provide further corrections todescribe the behavior of the exact closures (3.20) in the small q limit.In the asymptotic limit q → ∞ , the eigenfrequencies ω n are complex, with theimaginary part given by (cid:61){ a n } / ∂ (cid:61){ ω , } ∂k ∼ ± (cid:114) Pe k ∓ Pe , (3.25) (cid:61){ ω , } k ∼ ± (cid:114) Pe k ∓ Pe , (3.26)They differ only by a constant pre-factor of two, independent of k and Pe , signalingdispersive wave behavior. However, as k → ∞ , the propagating modes become non-dispersive, in the kinematic sense. Such behavior is typical of all eigenmodes in the q → ∞ regime.For the first four eigenmodes, q → ∞ , and m = 0, the closures are Λ (0 , ∼ (cid:114) k Pe − Pe ± i (cid:34) (cid:114) k Pe − k (cid:35) ,Λ (4 , ∼ (cid:114) k Pe − Pe ± i (cid:34) (cid:114) k Pe − k (cid:35) , (3.27)where we have paired the closures of coalesced eigenvalues for simplicity of expression. Inphysical space, the closure Λ (0)0 (with plus sign in imaginary component), can be written hear Dispersion − U ∂ (cid:104) θ (0) (cid:105) ∂x − π √ κU M ∂ / (cid:104) θ (0) (cid:105) ∂x / . (3.28)where we have used (1 + i ) / √ √ i and ( ik ) α ˆ f ( k ) ←→ ∂ α f /∂x α when writing thederivative of fractional order with α = 1 / et al. ± ik ) term in (3.27). The second term in (3.28) is referred to in the literature asa semi-derivative or half-derivative operator (see Oldham & Spanier 1974 Ch. 3). Itcaptures both the real part of the closure associated with the decay rate of the mode,and the imaginary part associated with the (kinematically) dispersive behavior of themodes (e.g. see 3.25).3.4. Approximation to Closure Operator for Arbitrary Wavenumber and P´eclet number
The asymptotic approaches above fail to provide a closure that is continuous acrossthe EPs in wavenumber space, as they only apply for values of q far from EPs. To derivea continuous closure operator, we approximate the eigenvalues as coalescing pairs anddivide the parameter dependence into three regions based on the behavior of eigenmodes:I) A region up to the EP that captures the merging of the eigenvalues, q (cid:54) q EPl . II) Anintermediate matching region for q EPl < q (cid:54) q EPl . And III) a region that matchesthe asymptotic behavior of eigenvalues q > q EPl . The parameter dependence of theeigenvalues in regions I and II is captured by a quadratic approximation to the coalescingpair based on the equation describing a half ellipse † , shown in Figure 5.To capture the coalescing eigenvalue pair that occurs in region I, two pieces ofinformation are needed: 1) The location of each eigenvalue at q = 0, which is given exactlyby a n = (2 n ) , and 2) the location of the EPs on the q -axis (i.e. q = q EP(cid:96) ), which can allbe easily calculated independent of any (tunable) physical parameter (Blanch & Clemm1969; Hunter & Guerrieri 1981). The quadratic approximation to each pair is therefore Γ l = γ l ± ∆ l | q EPl | (cid:113) | q EPl | − k Pe , (3.29)where Γ l ≈ { a l , a l +2 } approximates the eigenvalue pair, q EPl is the location of the EPon the q -axis (Table 1), and ∆ l = 8 l + 2 and γ l = 16 l + 8 l + 2 are known constants thatmatch the q -dependence of the eigenvalue pair (see Fig. 5 for a geometrical interpretationof the terms in the quadratic approximation).Approximation captures the coalescing of eigenvalue pairs in region I and providesa continuous transition across (the exact) EPs in wavenumber space from purely real(decaying) to purely imaginary (propagating) in region II. For values of q > q EP(cid:96) thereal part of the approximation if a constant ( γ (cid:96) ), in disagreement with the behavior ofthe eigenvalues (see Fig. 2). Furthermore, the approximation (3.29) underestimates theamplitude of imaginary part of the coalesced eigenvalues at large q (also beyond the EP),as it approximates the modes as (kinematically) non-dispersive waves. Thus, in regionII we add a term based on the asymptotic behaviour of the eigenvalues at large q , andregion III matches the asymptotic behavior of the coalesced eigenvalue pair at large q.With the above corrections, the approximation to the exact closure operator (3.20)that remains continuous in wavenumber space across EPs, captures the small and large † see also Ziener et al. M. A. Jimenez-Urias and T. W. N. Haine
Figure 5.
Quadratic approximation to the coalescing eigenvalue pair using a half ellipse.With the known value of a n = (2 n ) at q = 0 for all n = 0 , , , · · · , we have ∆ l = (4 l + 2) − (4 l ) = 8 l + 2 and γ l = (4 l ) + ∆ l = 16 l + 8 l + 2. These parameters areindependent of Pe and k . For q > q EPl , the square root becomes purely imaginary and, as q → ∞ , it approximates the linear asymptotic behavior of the coalesced eigenvalues, namely ∼ iq at large q (see Appendix B). q -limits, and is valid for all t > Λ l ) = γ l Pe ± ∆ l | q EPl | Pe (cid:113) | q EPl | − k Pe if 2 k Pe (cid:54) | q EPl | , β l (cid:114) k Pe ± α (cid:96) ∆ l | q EPl | Pe (cid:113) | q EPl | − k Pe if | q EPl | < k Pe (cid:54) | q EPl | , β l (cid:114) k Pe ± i (cid:32) β l (cid:114) k Pe − η l k (cid:33) if 2 k Pe > | q EPl | . (3.30)In the above approximation, region II is delimited by the range q EP(cid:96) < q < q EP(cid:96) ,although a different upper limit could have been chosen. All of α (cid:96) , β (cid:96) and η (cid:96) are constantsthat match the asymptotic limits derived in (3.27) in a way that assures the closures arecontinuous across regions, although a different definition of upper limit in region II wouldaffect these values. Fig. 6 shows the closure approximation to the first two eigenvaluepairs l = 0, l = 1.The closure associated with the square-root term in regions I and II in (3.30) thatcaptures the coalescing of eigenvalues retains the limits of large P´eclet number. This is,applying L’Hopital’s rule lim P e →∞ (cid:113) | q EPl | − k Pe Pe ≈ − k Pe . (3.31)The limit is proportional to the effective diffusivity derived by G. I. Taylor as demon-strated in (3.22), (3.23) and (3.24). Furthermore, the argument in the square-root in(3.31) must remain positive in order to represent a purely diffusive process. This implies4 k Pe (cid:47) .
468 (as q EP ≈ . i ), which is a statement of strong separation of scalesin the large P´eclet number limit, in line with homogenization theory (Majda & Kramer1999; Haynes & Vanneste 2014).The operator in (3.30) that captures the eigenvalue dependence on q = 2 ik Pe acrossEPs, i.e. across regions I and II, is a fractional order operator in physical space. It can hear Dispersion Figure 6.
Quadratic approximation across EPs to eigenvalue pairs { a , a } (red and blue linesin (c) and (d)), and { a , a } (green and purple lines in (a) and (b)) in the three regions asdescribed in (3.30). The parameters associated with the quadratic approximation (3.29) are ∆ = 2, ∆ = 10 and γ = 2 and γ = 26. Depicted in the figures are the closure equations inregions II and III. be written 1 Pe (cid:114) | q EPl |I + 4 Pe ∂ ∂x (3.32)where I is an identity operator (when discretized, it is an identity matrix). The operatorsthat capture the asymptotic limit for large q -values, i.e. those is region III in (3.30) aregiven by the square-root derivatives of the form (3.28).
4. Discussion
The results obtained in this study are based on the exact solution to the problem oftracer dispersion by a periodic shear flow for all t >
0. Our solution method has no a priori assumption of separation of scales, but the periodic shear flow is time-independent andso the velocity field does not evolve according to the Navier Stokes equation, a commonpractice in studies of the advection-diffusion equation (Majda & Kramer 1999). As aresult, the closures derived are kinematic, and not dynamic. We neglect any externalforcing and reaction terms, but they are readily incorporated. Our method of solutionrelies on an ansatz (2.6) that reduces the original advection-diffusion equation into aperiodic, second order ordinary differential equation for the amplitude. When the velocityfield is defined by a single cosine mode, as in this study, the amplitude equation reducesto the Mathieu equation (2.9), which depends on a single parameter q . The general casespecifies the velocity field by a sum of (infinitely) many modes. The resulting amplitudeequation is Hill’s equation, from which Mathieu’s equation is the simplest case, and bothare subject to Floquet theory (Olver et al. M. A. Jimenez-Urias and T. W. N. Haine complexity of the system. Nonetheless, Hill’s equation generalizes the problem studiedhere and thus represents a promising avenue of research. Solutions experience a transitionin wavenumber space across Exceptional Points, which are associated with the mergerof eigenvalue pairs. Exactly at at the EPs, one eigenfunction is no longer independent,and so the eigenfunctions no longer form a complete set in the space of square-integrablefunctions. Nonetheless, a generalized eigenfunction can be constructed that completes theset, which solves the shear dispersion problem for all q , and thus all possible combinationsof k and Pe (see Brimacombe et al. n determines the solution. As a result, Φ in the initial condition must have a series expansion in terms of ce n , i.e. a cosine Fourierseries. In the case of Dirichlet boundary conditions, the eigenfunctions are the odd, π -periodic sine-elliptic Mathieu functions se n +2 , with eigenvalues denoted by b n +2 . Eventhough ce n and se n +2 do not share the same eigenvalues ( a n (cid:54) = b n +2 ), the asymptoticlimits at small and large q of the two have the same dependence with q , with only slightdifferences in the coefficients (see Olver et al. q = 2 ik Pe is insensitive to the choice of boundary conditions, the eigenvalue system (inAppendix A) remains non-Hermitian, and EPs still occur, albeit at different q values.The P´eclet number dependence of our derived closures elucidates tracer transport byjets. Smith (2005) postulated a predictive theory associated with tracer transport bycoherent jets in β -plane turbulence, in which the across-jet transport is described bymixing-length theory and the along-jet transport is controlled by shear dispersion, witha background diffusivity determined by the across-jet mixing. Such theory, valid in thelimits of homogenization ( Pe large but | q | (cid:28) q ), the tracer experiences linear advectionthat is independent of the diffusivity, as well as a non-local advection and non-localdiffusion, both with a coefficient that is proportional to √ κU . Thus an increase tothe background diffusivity associated with the across-jet mixing, increases the coefficientassociated with the along-jet non-local diffusion. An attractive problem to consider is howflow instabilities or background turbulence associated with coherent jets can introducespatial inhomogeneities to the background diffusivity, and whether a closed equationof the form (3.28) can be written for a new spatially varying κ associated with thebackground turbulent flow. Lastly, it remains to be seen how the derived closures forlarge Pe values, which can be considered Eulerian, fair against diagnostic Lagrangianmethods that provide instantaneous eddy diffusion coefficients even in the presence ofchaotic advection (Nakamura 1996).
5. Conclusion
In this study we have provided, to our knowledge, the first analytical solution to theproblem of shear dispersion for all t >
0. We thereby calculate exact closure operators(3.20) for a wide range of averaging kernels, and derive analytical approximations to theexact closures (3.30), which are expressed as the sum of two differential operators offractional order. The approximations capture the behavior of the solution for large andsmall P´eclet number, and are continuous in wavenumber space. These results establishthe limits of validity of G. I. Taylor’s effective diffusivity κ ∗ ∝ U /κ , where U is thecharacteristic flow speed and κ is the molecular diffusivity. It is valid in the small P´eclet hear Dispersion k >
0. Itis also valid in the large P´eclet number limit, provided that the wavenumber of theinitial condition satisfies 2 k (cid:28) / Pe , where Pe is the P´eclet number. This requirementguarantees a wide scale separation between the flow and the tracer field, and is consistentwith homogenization theory. For wavenumbers that do not satisfy this requirement, weestablish a new large P´eclet number limit that is associated with propagating behavior inthe kinematic solution. The closure is then composed of a linear advection operator plusa differential operator of fractional order (a semi-derivative) that compactly represents(non-local) diffusion and advection, with a coefficient proportional to √ κU . We findthat closures in agreement with homogenization theory (at large P´eclet number) cannotcapture the early time evolution of the averaged tracer, but it is accurately described bythe new large P´eclet number regime. Acknowledgments . The authors acknowledge helpful conversations with C. Mene-veau and A. Mani at an early stage of the research. This project was funded by the JohnsHopkins University Institute for Data Intensive Engineering and Science Seed FundingInitiative and NSF grant 1835640. The code used to calculate Mathieu functions, eigenval-ues, and solutions is available at github.com/Mikejmnez/Adv Diff Supplemental Materials . Appendix A. Calculation of Mathieu eigenfunctions andcharacteristic values
Following Chaos-Cador & Ley-Koo (2002); Olver et al. (2010); Ziener et al. (2012),Mathieu functions are calculated from an infinite, tridiagonal eigenvalue system thatarises after substituting the definition (2.11) into Mathieu’s equation, and grouping bycosine terms. We get for each cosine termcos(0 y ) : 0 = a n A (2 n )0 − qA (2 n )2 , cos(2 y ) : 0 = ( − a n ) A (2 n )2 − q (cid:104) A (2 n )0 + A (2 n )4 (cid:105) , cos(4 y ) : 0 = ( −
16 + a n ) A (2 n )4 − q (cid:104) A (2 n )2 + A n (cid:105) , cos(6 y ) : 0 = ( −
36 + a n ) A (2 n )6 − q (cid:104) A (2 n )4 + A (2 n )8 (cid:105) , cos(2 ry ) : 0 = (cid:0) − r + a n (cid:1) A (2 n )2 r − q (cid:104) A (2 n )2 r +2 + A (2 n )2 r − (cid:105) . (A 1)Therefore, √ q · · ·√ q q q q
00 0 q q . . . . . . . . . q r q . . . √ A (2 n )0 A (2 n )2 A (2 n )4 A (2 n )6 ... A (2 n )2 r ... = a n √ A (2 n )0 A (2 n )2 A (2 n )4 A (2 n )6 ... A (2 n )2 r ... , (A 2)8 M. A. Jimenez-Urias and T. W. N. Haine where q = 2 ik Pe . The elements of the eigenvector are the Fourier coefficients thatdetermine each Mathieu function ce n ( q, ˜ y ) as defined by (2.11). They satisfy the or-thonormality relationship (Ziener et al. ∞ (cid:88) n =0 A (2 n )2 p A (2 n )2 r = δ pr − δ p δ r . (A 3)This formula is used to test convergence of the tridiagonal eigenvalue system (A 2), as itmust be truncated to numerically calculate the eigenvalues and eigenfunctions. See (2.14)and (2.16) for other Mathieu function identities. Appendix B. Asymptotic Approximations to Mathieu’s Eigenvalues
Expansions of the eigenvalues a n ( q ) as a function of parameter q = 2 ik Pe give a = 2( k Pe ) + 78 ( k Pe ) − k Pe ) + O (( k Pe ) ) , (B 1) a = 4 −
54 ( k Pe ) − k Pe ) + 10024011244160 ( k Pe ) + O (( k Pe ) ) , (B 2) a = 16 −
215 ( k Pe ) + 43354000 ( k Pe ) − k Pe ) + O (( k Pe ) ) , (B 3) a = 36 −
235 ( k Pe ) + 1872744000 ( k Pe ) + 67436171452124800000 ( k Pe ) + O (( k Pe ) ) , (B 4)which are valid for k Pe (cid:28) / et al. q , the limitlim q → a n ( q ) = (2 n ) ensures convergence of Mathieu functions into Fourier functions.For large q , the first four eigenvalues asymptote to a , a ∼ √ k Pe − ± i (cid:16) √ k Pe − k Pe (cid:17) , (B 5) a , a ∼ √ k Pe − ± i (cid:16) √ k Pe − k Pe (cid:17) (B 6)(Hunter & Guerrieri 1981; Ziener et al. REFERENCESAris, Rutherford
Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences (1200), 67–77.
Arscott, Felix Medland
Periodic differential equations: an introduction to Mathieu,Lam´e, and allied functions , , vol. 66. Elsevier.
Batchelor, Go K
Journal of Fluid Mechanics (2), 177–190. Bender, Carl M
Physics Reports (1-3), 27–40.
Bender, Carl M & Boettcher, Stefan
Physical Review Letters (24), 5243. Blanch, G & Clemm, DS
Mathematics of Computation (105), 97–108. Brimacombe, Chris, Corless, Robert M & Zamir, Mair arXiv preprintarXiv:2008.01812 . Camassa, Roberto, Lin, Zhi & McLaughlin, Richard M
Communications in Mathematical Sciences (2),601–626. hear Dispersion Chaos-Cador, Lorea & Ley-Koo, E
Revista mexicana de f´ısica (1), 67–75. Chen, Wen
Chaos: An Interdisciplinary Journal of Nonlinear Science (2), 023126. Deconinck, Bernard, Trogdon, Thomas & Vasan, Vishal siam REVIEW (1), 159–186. Gent, Peter R & Mcwilliams, James C
Journal of Physical Oceanography (1), 150–155. Gent, Peter R, Willebrand, Jurgen, McDougall, Trevor J & McWilliams, James C
Journalof Physical Oceanography (4), 463–474. Grooms, Ian, Smith, K Shafer & Majda, Andrew J
Dynamics of atmospheres and oceans , 95–107. Haynes, PH & Vanneste, J arXiv preprint arXiv:1401.6665 . Heiss, WD
Journal of Physics A:Mathematical and General (6), 2455. Heiss, WD
Journal of Physics A: Mathematical andTheoretical (44), 444016. Hinch, EJ
Perturbation Methods . Cambridge University Press and references therein.
Holmes, Mark H
Introduction to perturbation methods , , vol. 20. Springer Science &Business Media.
Hunter, C & Guerrieri, B
Studies in Applied Mathematics (2), 113–141. Lischke, Anna, Pang, Guofei, Gulian, Mamikon, Song, Fangying, Glusa, Christian,Zheng, Xiaoning, Mao, Zhiping, Cai, Wei, Meerschaert, Mark M, Ainsworth,Mark & others
Journal of Computational Physics , 109009.
Majda, Andrew J & Kramer, Peter R
Physics reports (4-5), 237–574.
Mani, Ali & Park, Danah
Marshall, John & Speer, Kevin
Nature Geoscience (3), 171–180. McLachlan, N.W.
Theory and Application of Mathieu Functions . Oxford UniversityPress.
Mercer, GN & Roberts, AJ
SIAM Journal on Applied Mathematics (6),1547–1565. Miri, Mohammad-Ali & Alu, Andrea
Science (6422).
Mulholland, HP & Goldstein, S
The London, Edinburgh, and Dublin PhilosophicalMagazine and Journal of Science (53), 834–840. Nakamura, Noboru
Journal of the atmospheric sciences (11), 1524–1537. Oldham, Keith & Spanier, Jerome
The fractional calculus theory and applications ofdifferentiation and integration to arbitrary order . Elsevier.
Olver, Frank W J, Lozier, Daniel W, Boisvert, Ronald F & Clark, Charles W
NIST handbook of mathematical functions hardback and CD-ROM . Cambridge Universitypress.
Park, Danah, Shirian, Yasaman & Mani, Ali
Bulletin of the American Physical Society . Podlubny, Igor
Fractional differential equations: an introduction to fractional derivatives,fractional differential equations, to methods of their solution and some of theirapplications . Elsevier. M. A. Jimenez-Urias and T. W. N. Haine
Rhines, Peter B & Young, William R
Journal of Fluid Mechanics , 133–145.
Samko, Stefan G, Kilbas, Anatoly A, Marichev, Oleg I & others
Fractionalintegrals and derivatives , , vol. 1. Gordon and Breach Science Publishers, YverdonYverdon-les-Bains, Switzerland.
Smith, K Shafer
Journal of Fluid Mechanics , 133–142.
Taylor, Geoffrey Ingram
Proceedings of the Royal Society of London. Series A. Mathematical and PhysicalSciences (1137), 186–203.
Tseng, Chien-Cheng, Pei, Soo-Chang & Hsia, Shih-Chang
Signal Processing (1),151–159. Young, WR a & Jones, Scott
Physics of Fluids A: Fluid Dynamics (5), 1087–1101. Ziener, CH, R¨uckl, M, Kampf, T, Bauer, WR & Schlemmer, Heinz-Peter
Journal of Computational and AppliedMathematics236