Numerical simulations of diffusion in cellular flows at high Péclet numbers

We study numerically the solutions of the steady advection-diffusion problem in bounded domains with prescribed boundary conditions when the Péclet number Pe is large. We approximate the solution at high, but finite Péclet numbers by the solution to a certain asymptotic problem in the limit Pe → ∞. The asymptotic problem is a system of coupled 1-dimensional heat equations on the graph of streamline-separatrices of the cellular flow, that was developed in [24]. This asymptotic model is implemented numerically using a finite volume scheme with exponential grids. We conclude that the asymptotic model provides for a good approximation of the solutions of the steady advection-diffusion problem at large Péclet numbers, and even when Pe is not too large.


1.
Introduction. In the simply connected, bounded domain Ω ⊂ R 2 consider the steady advection-diffusion problem where φ ε is a temperature. We assume here that the time-independent flow v is incompressible: ∇· v = 0, and it does not penetrate through the boundary of Ω: v · n = 0 at ∂Ω. The small parameter ε = Pe −1 ≪ 1 is the inverse of the Péclet number, the relative strength of the advection, compared to the diffusion. Since the flow v is assumed to be incompressible, there is a stream function Ψ(x, y) such that v = ∇ ⊥ Ψ = ∂Ψ ∂y , − ∂Ψ ∂x .
(2) solutions of the asymptotic model of [21] with numerical simulations, obtained in the MATLAB PDE Toolbox TM . We observe that the asymptotic model approximates well the solution of (1)-(3) even when the Péclet number is about 10 2 .
The paper is organized as follows. In the next section, the summary of the asymptotic model [21] is presented. The asymptotic system of heat equations is solved numerically using a finite volume method on exponential grids. We discuss its numerical implementation in Section 3. We then perform numerical tests in order to determine accuracy of the asymptotic model. The numerical results are presented in Section 4. 2. Asymptotic model. At high Péclet numbers the numerical solution to (1) becomes expensive if a discretization method does not take into account diffusive boundary layers. These boundary layers arise where the stream function Ψ(x, y) for the incompressible flow v (see (2)) vanishes: Ψ(x, y) = 0 2 . We propose here to use a detailed analysis of boundary layers for (1) developed in [21] as ε = 1/Pe → 0. Let us summarize some results obtained in [21].
As ε tends to zero, the solution to (1) converges to a solution of a coupled onedimensional heat equations on a graph. The edges of the graph, that arise in this diffusion on the graph model, are associated with the separatrices between different flow cells, and the vertices of this graph are associated with saddle-points of the stream function Ψ(x, y). See e.g. the two-cell case on Figure 1.
The diffusion on the graph model is independent of Péclet number, and the relative H 1 error between the solution to (1) and the solution to the diffusion on the graph model was estimated a priori as O(ε 1/2 ) (Theorem 6.1 in [21]). Let us first describe the simplest case, when we have just two vertices.
2.1. The two-cell case. We describe the asymptotic problem first on the simplest example of a domain Ω that consists of two cells C 1 and C 2 depicted in the left figure of Figure 1. We denote by e j0 = ∂Ω ∩ ∂C j , j = 1, 2, the part of the boundary of Ω along the cell C j and by e 12 the common edge of the two cells. We also introduce the boundary layer coordinates h and θ 12 , θ j0 , j = 1, 2. The coordinate h = Ψ/ √ ε, and it represents the (rescaled) distance to the level-set {(x, y)|Ψ(x, y) = 0}. The coordinate θ 12 represents parametrization along the edge e 12 = {h = 0} ∩ {0 ≤ θ 12 ≤ l 12 }, while the coordinates θ j0 parameterize along the boundaries e j0 = {h = 0} ∩ {l 12 ≤ θ j0 ≤ l j0 }. We first solve the heat equation "along e 12 ": with a prescribed initial data f 0 12 and A natural way to impose this decay to constants is to prescribe homogeneous Neumann boundary conditions: Then we solve two half-space problems "along the outer boundaries e j0 " with the prescribed Dirichlet data that comes from (3): and with (7) h = ±∞, and with the Dirichlet data f j0 (0, θ j0 ) = T 0 (θ j0 ) at h = 0. The initial data for (8) and (9) comes from (5): Finally we glue together the functions f 10 (h, l 10 ), h ≤ 0 and f 20 (h, l 20 ), h ≥ 0: The asymptotic problem is to construct a periodic solution of the above, that is, find a function f 0 12 (h) so that f 0 12 (h) = f g 12 (h), h ∈ R. This problem is described schematically in the right figure of Figure 1. There exists a unique function f 0 12 ∈ L 2 (R) such that f 0 12 = f g 12 (Proposition 5.1 in [21]). An alternative approach to the proof of existence of a periodic solution of (5)- (11), that is somewhat less transparent in the two-cell case but is easier to generalize to the case of multiple cells is as follows. Let us define the operator L 12 : L 2 (R) → L 2 (R) by L 12 : f 0 12 → f 12 (l 12 ), that is, the solution operator of (5). We also let L 10 and L 20 be solution operators for (8) and (9), respectively with homogeneous boundary data T 0 = 0. We introduce an operator L = L 12 ⊗ L 10 ⊗ L 20 defined on We also define a redistribution operator R on the same space Above the operators R ± restrict a function defined on R to the positive and negative semi-axes, respectively, while the gluing operator G glues together two functions defined on those axes: as in (11). Then the existence of a unique function f 0 12 ∈ L 2 (R) such that f 0 where the function g(h) is obtained by solving (5)-(11) with f 0 12 = 0 and inhomogeneous boundary conditions.

2.2.
The general multiple-cell case. Equation (14) has a straightforward generalization to the case of more than two cells. Suppose the domain Ω consists of a finite number of cells. The asymptotic model is described in terms of an oriented graph constructed using the stream function Ψ as shown on Figures 2. The vertices of this graph are associated with the saddle points of Ψ. The edges e ij of the graph are associated with the separatrices of the the stream function. The direction of an edge is determined by the direction of the velocity field on the corresponding separatrix. The length of an edge is determined by the length of the separatrix in the boundary layer coordinate θ associated with Ψ. The boundary edges are those that are associated with the separatrices at the boundary of the domain. The cells C i are quadrangles bounded by minimal cycles of the graph. The interior edges (drawn as solid arrows on the right figure of Figure 2) are indexed so that a common edge of two cells C i and C j is denoted by e ij . The boundary edges (drawn as dotted arrows on the right figure of Figure 2) are indexed so that the outer part of a boundary cell C i is denoted by e i0 . The boundary value problem is: Given the values of the temperature T 0 on the boundary edges e i0 , determine the values of the temperature f ij on all the edges. Note that the value of f ij may vary along each edge.
Given the values of f ik on all the edges, find the solutions f i of the Childress' problem [3] for each cell C i : where the index k takes four values of the adjacent cells, l i = l ik1 + · · · + l ik4 is the length in θ of the four edges e ik1 ,. . . e ik4 , bounding C i and f ik (θ) = f ik1 (θ), . . . , f ik (θ) = f ik4 (θ) are the values of the temperature on respective edges.

YULIYA GORB, DUKJIN NAM AND ALEXEI NOVIKOV
• [iii] When any two cells C i and C j share a common edge, the normal derivatives from the left and from the right match point-wise on this edge:  (1) is the square (0, π) × (0, π) in R 2 . We choose Ψ(x, y) = sin(k 1 πx) sin(k 2 πy) with integer k 1 and k 2 . The domain Ω will contain either one (k 1 = 1, k 2 = 1), two (k 1 = 2, k 2 = 1), or four flow cells (k 1 = 2, k 2 = 2). The second boundary layer coordinate is defined as before: h = Ψ/ √ ε. In the case of two cells the only interior separatrix is the line x = π 2 (as on the left plot of Figure 5). For the four-cell problem the interior separatrices are the lines x = π 2 and y = π 2 (as on the right plot of Figure 5). With a slight abuse of notation we use f for the periodic solutions of (15) as well as for the numerical approximations of (15) in the rest of this paper.
For the spatial discretization we use a finite volume method [7], and we assume that the solution is piecewise constant. The finite difference backward Euler method is used for the temporal discretization.
3.1. Boundary layer coordinates. It was shown in e.g. [3,8,21] that the solution of (1)-(4) can be well-approximated by the stretched asymptotic boundary layer solution f ε (x, y) = f (h(x, y), θ(x, y)), where f (h, θ) is the unique solution of the Childress' problem (15). In this section we briefly discuss the construction of the local coordinates (h, θ). We concentrate on the 2-cell case (see the left plot in Fig.  5), because the one-cell case is discussed in [8, p.349], and the four-cell case is similar to the one-cell case.
Thus we use a standard exponential grid [25,11] to obtain a numerical discretization of (15). The grid is defined as follows. Let {0 = y 0 , y 1 , · · · , y N = 1} be the uniformly distributed points in the interval [0, 1]. Then for an appropriate positive constant C which is a stretching parameter, the positive part of the exponential grid points are defined by t j = −C ln y j , for each j = 1, · · · , N.
The negative part of the exponential grid points is obtained by reflecting the points {t j } N i=1 about the origin. All the grid points can be obtained by merging those two parts 3.3. Discretization of the heat equations. As described in Section 2.1, the Childress' problem consists of solving the heat equations along the common edges of the two cells (e.g. see equation (5)) and along the exterior boundaries (e.g. see equations (8), (9)).
Consider an interior edge. We impose the Neumann boundary conditions (7) at the endpoints of the interval as Integrating the equation over each subinterval where ∆x = x i+1 − x i . Since f is piecewise constant we denote its value on the subinterval (x i , x i+1 ) by f i+ 1 2 . Then (18) is Using the backward Euler method with the step size ∆θ, we obtain the difference equation at the (n + 1)th time step: for i = 2, · · · , 2N − 3. For i = 1, 2N − 2 we impose Neumann boundary conditions at end points to obtain: We similarly discretize the half-space problems (8) and (9). In the latter case we used only a half of the grid points: (8) and (9), respectively. Since the functions in finite volume method are piecewise defined the redistribution operator R (see (12)), and the gluing operator G (see (13)) can be easily implemented. Finally, we need to impose the corresponding boundary conditions on the outer edges, i.e. at h = 0.

3.4.
Discretization of the boundary conditions on the outer edges. As mentioned earlier, the homogeneous Neumann conditions and the Dirichlet conditions have been imposed on the outer edges. For the Neumann boundary conditions we use for i = N − 1. For the Dirichlet boundary conditions we use where for i = N − 1, and f n+1 N are the known values which are Dirichlet data.
4. Numerical results. In this section, we present some numerical results, and compare them with simulations obtained in the MATLAB PDE Toolbox TM . We start with the reduced Childress' problem (15) for the one-cell case, Ψ(x, y) = sin x sin y. Figure 3 illustrates the agreement between the initial condition (at θ = 0), and the solution after one period (at θ = 8). The which shows that we can find the periodic solution with our numerical scheme. In Figure 4 we plot this periodic solution. The left plot represents the solution of the asymptotic Childress' problem (15) in the boundary layer coordinates h > 0, 0 ≤ θ ≤ 8. The right plot shows the same solution with ε = 10 −3 in the Cartesian coordinates.

4.1.
The case of the mixed Dirichlet-Neumann boundary conditions. The domain for the two-cell case with Ψ(x, y) = sin 2x sin y is shown on the left part of Figure 5. In the boundary layer coordinate system, as ε → 0, the problem for the two-cell case becomes the following system of equations over the rectangular region h ∈ R, 0 < θ < 10. First, for 0 < θ < 4: with the homogeneous Neumann boundary conditions at ±∞ and an unknown initial data f 0 (h) which will be determined below in (27).   For 4 < θ < 5 we solve two problems for h > 0 and h < 0, respectively: f + (0, θ) = π, 5 < θ < 9, ∂f + ∂h (0, θ) = 0, 9 < θ < 10, and with the initial data where f solves (24), and the homogeneous Neumann boundary conditions at ±∞. Hence due to the periodicity we have The stretching parameter C (16) for the exponential grid was chosen to be 5. In Figure 6 we present the asymptotic solution for the two-cell problem (24)- (26). The left plot shows that we indeed obtain a periodic solution. The right plot shows a comparison with the solution of (1) with ε = 10 −3 .
In the four-cell case the structure of the separatrices is as in Figure 5. In the right plot of Figure 7 we show the 2-D solution of (1) in the Cartesian coordinates. In the left plot of Figure 7 we compare this 2-D solution with our asymptotic solution of the Childress' problem in the boundary layer variables and at fixed θ = 0.   For these boundary conditions, the system of the heat equations (24)-(26) is changed slightly, and it is as follows. For 0 < θ < 4 the equation is as before (24). For 4 < θ < 5 we consider splitting for positive and negative half spaces: ∂f + ∂θ − ∂ 2 f + ∂h 2 = 0, h > 0, 4 < θ < 10, f + (0, θ) = π − arcsin( √ 5 − θ), 4 < θ < 5, f + (0, θ) = π, 5 < θ < 9, f + (0, θ) = π − arcsin( √ θ − 9), 9 < θ < 10, and with homogeneous Neumann boundary conditions at h = ±∞. Figures 8 -11 show the solution of the system of the heat equations (24), (28), (29). In particular, the curve in the left plot in the Figure 8 is the periodic initial condition which is obtained by solving the Childress' problem for the two-cell case. Finally, in Figure 12 we compare our solution (solid curve from the left plot of Figure 8) with the corresponding solution of (1) with ε = 10 −3 .

Numerical convergence.
In this section we present the numerical convergence of the asymptotic problem on the graph given by (15). The convergence tests are presented in Tables 1-3 for each of the three cases discussed in Subsections 4.1, 4.2. In these tables, N denotes the number of grid points of the exponential grid in the spatial domain [−30, 30], and T is chosen so that 1/T is the "time-step" size in the circulation variable θ to solve the heat equations. The solution for the case of N = 800 and T = 400 is used as a reference to calculate the errors of solutions over coarser grids and with bigger "time-step" sizes. As it is expected from our first-order discretization, we observe the O(1/N, 1/T ) numerical convergence. Higher order methods can be used to solve (15), but we  (15) approximates (1) for large Péclet numbers. We emphasize that unlike other numerical schemes for advection-diffusion problems (such as e.g. upwind methods which are not suited for problems with high Péclet number) our numerical scheme was developed for the asymptotic problem (15), which is independent of the Péclet number [21]. Thus it provides an effective numerical tool to solve the problem at high Péclet numbers. We now turn to assessing the accuracy of the numerical scheme. In the next section we discuss the convergence of our method in terms of ε = Pe −1 ≪ 1.   Figure 14. Comparison of the asymptotic solution at θ = 2 or (x, y) = (π/2, π/2) with the solution of (1) with various ε (left), a magnified region (right). PDE Toolbox TM for various ε = 10 −1 , 10 −2 , 10 −3 . For smaller ε the MATLAB's solution becomes unreliable.
The comparison results are shown for the two-cell case (24)- (27) at two different values of θ: for θ = 0 in Figure 13, and for θ = π/2 in Figure 14. In order to put all four functions on the same graph we rescaled the solutions of (1) obtained with MATLAB TM . As it is evident from these plots the asymptotic model is reliable for ε as large as 10 −2 . For ε = 10 −1 shows less fidelity, but the approximate solution is still fairly close to the solution of (1).
We used a finite volume method with exponential grid to discretize the asymptotic system of the heat equations (15). It turns out that piecewise constant functions in our finite volume method allow us to handle complicated restriction and gluing transformations easily. By imposing homogeneous Neumann conditions at end points of the restricted domains (17), we can capture the constant states (6) automatically. Figure 15. Approximate solutions by the "water-pipe network". Spectral methods, e.g. Galerkin approach [25,11] are used for the problems on unbounded domains because spectral approach has lower computational cost and higher accuracy. It is, however, difficult to implement for our problem, because complicated gluing conditions, described in Section 3, need to be explicitly specified for Galerkin basis functions, i.e. Laguerre and Hermite polynomials. We are currently implementing a spectral method with Laguerre polynomials for advection-diffusion at high Péclet numbers.
It was shown in [21] that it was also possible to approximate the solutions of (1) by the "water-pipe network". Namely, restrict the domain to the region of width K √ ε near the separatrices for some positive fixed number K, and denote it by Ω ε K = Ω ∩ {|Ψ(x)| ≤ K √ ε}. The boundary of Ω ε K consists of ∂Ω and the level curves l K = {x ∈ Ω : Ψ(x) = K √ ε}. Then we approximate the solution of (1) by the solution of ε∆φ ε K − u · ∇φ ε K = 0, x ∈ Ω ε K , with the boundary conditions We performed numerical tests for the "water-pipe network" using MATLAB PDE Toolbox TM . We focused on the one-cell and four-cell cases in these tests. The square (0, 1) × (0, 1) and (−1, 1) × (−1, 1) in R 2 are considered as the domains for one-cell and four-cell cases, respectively. The stream function Ψ(x) by Ψ(x, y) = sin(πx) sin(πy)/ √ 2, and the small molecular diffusivity by ε = 10 −2 and ε = 10 −3 for one-cell and four-cell case, respectively. The positive number K was chosen by 5 in both cases. The results of numerical simulations are depicted in Figure 15. For such relatively large ε the "water-pipe network" approximates well the solution of (1). It turns out, however, that this method is not accurate enough for smaller ε due to difficulty in mesh generation in the narrow pipes.