A Computational Study of Residual KPP Front Speeds in Time-Periodic Cellular Flows in the Small Diffusion Limit
aa r X i v : . [ n li n . C D ] J un A Computational Study of Residual KPP FrontSpeeds in Time-Periodic Cellular Flows in theSmall Diffusion Limit
Penghe Zu, Long Chen, Jack Xin ∗ Abstract
The minimal speeds ( c ∗ ) of the Kolmogorov-Petrovsky-Piskunov (KPP)fronts at small diffusion ( ǫ ≪
1) in a class of time-periodic cellular flowswith chaotic streamlines is investigated in this paper. The variationalprinciple of c ∗ reduces the computation to that of a principal eigenvalueproblem on a periodic domain of a linear advection-diffusion operator withspace-time periodic coefficients and small diffusion. To solve the advec-tion dominated time-dependent eigenvalue problem efficiently over largetime, a combination of finite element and spectral methods, as well asthe associated fast solvers, are utilized to accelerate computation. In con-trast to the scaling c ∗ = O ( ǫ / ) in steady cellular flows, a new relation c ∗ = O (1) as ǫ ≪ Keywords:
KPP front speeds, time periodic cellular flows, chaotic stream-lines, sub-diffusion, residual front speeds.
PACS : 42.25.Dd, 02.50.Fz, 05.10.-a, 02.60.Cb ∗ Department of Mathematics, UC Irvine, Irvine, CA 92697, USA.Corresponding author: Penghe Zu. Phone number: 9493784420.Emails: (pzu,chenlong,jxin)@math.uci.edu.
Introduction
Front propagation in complex fluid flows arises in many scientific areas suchas combustion [22, 27], population growth of ecological communities (plankton)in the ocean [1], and reactive chemical front in liquids [21, 28]. An interestingproblem is the front speed enhancement in time dependent fluid flows withchaotic streamlines and random media [21, 28].In this paper, we shall consider a two dimensional time dependent cellularflow:(cos(2 πy ) + θ sin(2 πy ) cos( t ) , cos(2 πx ) + θ sin(2 πx ) cos( t )) , θ ∈ (0 , , (1.1)whose steady part (cos(2 πy ) , cos(2 πx )) is subject to time periodic perturba-tion that causes transition to Lagrangian chaos. Chaotic behavior of a similarflow field is qualitatively analyzed by formal dynamical system methods in [8].The enhanced residual diffusion in (1.1) is observed numerically in [6]. Theenhanced diffusion and propagation speeds in steady cellular flows with orderedstreamlines and their extensions have been extensively studied [2,3,9,12,14–16,20, 23, 25, 29, 30] among others. We shall take a statistical look at the chaoticstreamlines of (1.1) in terms of the scaling of the mean square displacementsfrom random initial data, and quantify the Lagrangian chaos of (1.1) as a sub-diffusion process.Consider then the advection-reaction-diffusion equation: ∂ t u = ǫ △ u + ~B ( x, t ) · ∇ u + 1 τ f ( u ) , x ∈ R , t > , (1.2)where f ( u ) = u (1 − u ) is the Kolmogorov-Petrovsky-Piskunov (KPP) nonlin-earity, ǫ is the molecular diffusion parameter, τ is the reaction rate and ~B ( x, t )is a space-time periodic, mean zero, and incompressible flow field such as (1.1).If the initial data for u is nonnegative and compactly supported, the large timebehavior of u is an outward propagating front, with speed c ∗ = c ∗ ( ~e ) in thedirection ~e = h , i . The variational principle of c ∗ is [18]: c ∗ = inf λ> µ ( λ ) λ , (1.3)where µ ( λ ) is the principal eigenvalue of the periodic-parabolic operator: L λ Φ = ǫ △ Φ + ( ~B + 2 λ~e ) · ∇ Φ + ( ǫλ + λ ~B · ~e + τ − f ′ (0))Φ − Φ t , (1.4)on the space-time periodic cell Ω. The principal eigenvalue µ ( λ ) can be com-puted by solving the evolution problem: w t = ǫ △ w + (2 ǫλ~e + ~B ( x, t )) · ∇ w + ( − C M + ǫλ + λ~e · ~B ( x, t ) + 1 τ f ′ (0)) ww ( x,
0) = 1 . (1.5)1ere the constant C M = ǫλ + λ || ~B || ∞ + || τ f ′ (0) || ∞ , so that w remains positiveand bounded by one. The µ ( λ ) is then given by: µ ( λ ) = C M + lim t →∞ t ln Z Ω w ( x, t ) dx (1.6)The number µ ( λ ) is also the principle Lyapunov exponent of the parabolicequation (1.5), and the formula (1.6) extends to the more general case when ~B is a stationary ergodic field [17]. The limit then holds almost surely and µ isdeterministic [17]. KPP fronts are examples of the so called “pulled fronts” [24]because their speed is determined by the behavior of the solution far beyondthe front interface, in the region where the solution is close to zero (the unstableequilibrium). The minimal speed of the planar wave solution of the linearizedequation near the unstable equilibrium gives (1.3), also known as the marginalstability criterion [24, 28].Equation (1.5) and formula (1.6) will be discretized for approximating µ atsufficiently large time for a range of λ values. The minimal point of µ ( λ ) /λ is searched by the golden section algorithm. At each new search of λ , thelarge time solution of (1.5) is computed. Because ǫ is small, (1.5) is advectiondominated. To this end, upwinding type finite element methods (EAFE [31])and spectral method (SM) are utilized. Spectral method has high accuracy butthe computation is slower because it takes small time step due to the restrictionof stability. On the other hand, EAFE can be solved faster using relative largertime step with implicit discretization. Without the compromise of accuracy, ourstrategy to obtain the minimal point is to narrow down the search interval byEAFE first, and solve more accurately in a smaller interval by SM.The rest of the paper is organized as follows. In section 2, we illustratethe difference of streamline geometry and dynamical properties of steady andunsteady cellular flows. We find that the chaotic streamlines of the time periodiccellular flow (1.1) can be quantified as a sub-diffusive process with distinct θ -dependent scaling exponents. The resulting motion appears ergodic inside theinvariant infinite channel domains. In section 3, we discuss numerical methodsfor advection-dominated problems such as EAFE, and semi-implicit SM, as wellas the stability of these methods and time step constraints. In section 4, weshow numerical results of KPP front speeds by a combination of these methodsand the existence of residual front speed in the sense that c ∗ ( ǫ ) = O (1) > ǫ ↓
0. In contrast, c ∗ = O ( ǫ / ) in steady cellular flows [20]. The streamlines ofthe steady cellular flows are ordered with closed orbits leading to a much slowerrate of enhancement for the transport. We also observe the presence of layer andcircular structures in the generalized eigenfunction w at large time, a reflectionof the advection-dominated transport in time periodic cellular flow (1.1). Thecorresponding energy spectrum of w shows decay towards high frequencies withan intermediate scaling range. In section 5, we give concluding remarks aboutour findings. 2 Properties of Cellular Flows
In this section, we illustrate the difference of streamline geometry and dynamicalproperties of steady and unsteady cellular flows. In particular, we quantify thechaotic behavior of the motion along streamlines of (1.1) using the empiricalmean square distance inside infinite invariant channels in the direction (1 , E [ | X ( t ) · (1 , | ], where X ∈ R denotes the particle trajectories in the flow. Weshow computationally that it scales with time as O ( t p ), p ∈ (0 , θ ∈ (0 , A phase portrait of the steady cellular flow B s = (cos(2 πy ) , cos(2 πx )) , (2.1)with a cell square is presented in Fig. 1. The phase portrait away from thesquare simply repeats. The streamlines are ordered. The closed orbits formelliptic part of the phase space. The saddles are located at half-integer points( n, m ) / n , m ∈ Z with connecting separatrices forming hyperbolic part of thephase space. Any particle trajectory is either a closed orbit or a separatrix,hence the motion is bounded inside the square. For the unsteady flow B u = (cos(2 πy ) + θ sin(2 πy ) cos( t ) , cos(2 πx ) + θ sin(2 πx ) cos( t )) , (2.2)the invariant manifold consists of lines y = x + n , n ∈ Z . The flow trajectoriesare restricted in the channels bounded by two neighboring lines y = x + n + 1and y = x + n . At θ >
0, the flow trajectory extends itself from one cell squareto another. The Lagrangian particle undergoes chaotic motion, see Fig. 2 foran illustration. The computation is done by a fourth order symplectic scheme.The local flow direction is colored red. The intensified red region in the middleindicates that the particle spends a lot of time wandering in and out of the cellsthere. Fig. 3 plots a projected trajectory in the direction (1 , / √ trajectories X ( t, ω ) over time interval [0 , ω denotes the random sam-ples. We calculate the empirical mean square distance ( E [ | X ( t, ω ) | ]) and theprojected mean square distance E [( X ( t, ω ) · (1 , ], and plot them as a func-tion of time on the logarithmic scale to recover scaling laws. For efficiency,the samples are obtained by a 500 node parallel implementation of the stan-dard 4th order Runge-Kutta method. Fig. 4 shows that the mean squaredistances E [ | X ( t, ω ) | ] ∼ O ( t p ), p ≈ . , .
673 at θ = 0 . , . π y), sin(2 π x)), with a unit cell square.x y Figure 1: Ordered and localized streamlines of steady cellular flow (2.1) with acell square containing four vertices. −0.5 0 0.5 1 1.5 20.20.40.60.811.21.41.61.822.2 The orbit of steady flow −1 0 1 2 3 400.511.522.533.544.5 The orbit of steady flow
Figure 2: Disordered and extended trajectories in time periodic cellular flow(2.2) over time [0 ,
20] (left) and [0 , y = x + 1 and y = x . The local flow direction is colored red.4
100 200 300 400 500 600 700 800 900 1000−8−6−4−20246 time t X \ c do t ( , ) ^ T ) /\ s q r t ( ) Figure 3: A projected trajectory X ( t ) · (1 , / √ θ = 1, showing stochastic character.5ence belonging to the sub-diffusive regime. Similar sub-diffusive behavior canbe observed in Fig. 5 where the projected mean square distance scales like E [( X ( t, ω ) · (1 , ] ∼ O ( t q ), q ≈ . , .
786 at θ = 0 . , . l og ( E | X ( t ) | ) log(E|X(t)| ) vs log(t) when theta = 0.1 RobustFit: 0.48347 0 1 2 3 4 5 6 7−3−2−10123 log(t) l og ( E | X ( t ) | ) log(E|X(t)| ) vs log(t) when theta = 0.4 RobustFit: 0.67307 Figure 4: Mean square distances E [ | X ( t, ω ) | ] vs. time on the logarithmic scalewith robust fit showing sub-diffusive scaling O ( t p ), p = 0 . , .
673 at θ = 0 . . l og ( E | X ( t ) ⋅ ( , ) T | ) log(E|X(t) ⋅ (1,1) T | ) vs log(t) when theta = 0.1 RobustFit: 0.6045 0 1 2 3 4 5 6 7−3−2−101234 log(t) l og ( E | X ( t ) ⋅ ( , ) T | ) log(E|X(t) ⋅ (1,1) T | ) vs log(t) when theta = 0.4 RobustFit: 0.78634 Figure 5: Projected mean square distances E [( X ( t, ω ) · (1 , ] vs. time onthe logarithmic scale with robust fit showing sub-diffusive scaling O ( t q ), q =0 . , .
786 at θ = 0 . . Numerical approximation of advection-dominated problems is a topic of inde-pendent interest. We shall use the Edge-Averaged Finite Element (EAFE)6ethod [31] due to the nice discrete maximum principle obeyed by EAFE andthe corresponding fast multigrid solvers. On the other hand, pseudo-spectralmethods are a class of highly accurate numerical methods for solving partialdifferential equations. In practice, pseudo-spectral method has excellent er-ror reduction properties with the so-called “exponential convergence” being thefastest possible as long as the solution is smooth. In this section, we brieflyreview these two methods used in our simulation.Recall that we are solving the following evolution problem over the unitsquare domain Ω = [0 , × [0 ,
1] with periodic boundary condition and a constantinitial condition: w t = ǫ △ w + (2 ǫλ~e + B ) · ∇ w − ( C M − ǫλ − λ~e · B − τ f ′ (0)) w,w (0 , y, t ) = w (1 , y, t ); w ( x, , t ) = w ( x, , t ) ,w ( x, y,
0) = 1 . (3.1) We first decompose the domain Ω into a triangulation T which is a set oftriangles such that ∪ τ ∈T τ = ¯Ω and ˚ τ i ∩ ˚ τ j = ∅ . (3.2)For the unit square, we simply use the uniform mesh obtained by setting lengthof each triangle h .We describe the EAFE method using a simple advection-diffusion equation (cid:26) −∇ · ( ∇ u + β ( x ) u ) = f ( x ) , x ∈ Ω u = 0 x ∈ ∂ Ω . (3.3)Associated with each T h , let V h ⊂ H (Ω) be the piecewise linear finite elementspace. The space H (Ω) is defined as the subspace of H (Ω) with zero trace on ∂ Ω.Given any edge E in the triangulation, we introduce a function φ E definedlocally on E (up to an arbitrary constant) by the relation: ∂φ E ∂t E = 1 | t E | ǫ − ( β · t E ) , (3.4)where E is the edge connecting two vertices q i and q j , and t E is the vector suchthat t E = q i − q j .The EAFE formulation of the problem (3.3) is: Find u h ∈ V h such that a h ( u h , v h ) = f ( v h ) for any v h ∈ V h (3.5)where a h ( u h , v h ) = X T ∈T h ( X E ⊂ T ω TE | t E | Z E e φ E J ( u ) · t E ds δ E v h ) and (cid:26) ω E = P E ⊂ T cot θ TE ≥ J ( u ) = ǫ ∇ u + βu θ TE is the angle between the edges that share the common edge E and δ E is related to the tangential derivative along E .The authors of [31] show that if the triangulation is a so-called Delaunaytriangulation, i.e., the summation of two angles of an interior edge is less thanor equal to π , the matrix of the linear system generated by EAFE would be anM-matrix with nonnegative row sum. The detail of the scheme can be foundin [31].When applied to equation (3.1), notice that ~B is divergence free, the opera-tor can be written in the form of (3.3) and the reaction term is always positive.We apply the mass lumping method [26] so that the discretization of the reac-tion term makes positive contribution to the diagonal. Therefore the M-matrixproperty and consequently the discrete maximum principle still holds. Remark 1
We have tested the Box method [5], EAFE method [31] and stream-line diffusion method [11] on the advection-dominated problem. The accuracyof EAFE scheme is better than Box method, but lower than the streamline dif-fusion method. We chose EAFE instead of streamline diffusion scheme becausethe matrix of linear system generated by EAFE is M-matrix. The M-matrix notonly preserves the discrete maximum principle (so that the numerical solutionstays between 0 and 1) but also benefits from the fast solvers.
To accelerate the speed of computation, Algebraic Multi-Grid (AMG) solveris used to solve the linear algebraic equation arising from each time step. Whenthe matrix is an M-matrix, the corresponding AMG solver is proven to be effi-cient. Specifically for the advection-diffusion equations, a multigrid method pre-serving the M-matrix property in coarse level is developed in [13]. Among manyAMG software packages, we use AGMG (aggregation-based algebraic multigridmethod) package [19]. The numerical test shows that AMG solver is 8 timesfaster than the direct solver.
To describe the pseudo-spectral method, we consider the problem: u t = ǫ ∆ u + ~B ( x, y, t ) · ∇ u + Cuu ( x, y,
0) = 1; u (0 , y, t ) = u (1 , y, t ) , (3.6) u ( x, , t ) = u ( x, , t ) . (3.7)The derivatives of solutions are ∇ u = F − ( i~k F ( u )) and △ u = F − (( i~k ) F ( u )).Here F is the discrete Fourier transformation and ~k represents the wave num-bers.The pseudo-spectral scheme with semi-explicit scheme is that diffusion termuses implicit scheme and advection term uses half step lagged explicit scheme:ˆ U n +1 − ˆ U n ∆ t = − ǫ~k ˆ U n +1 + F ( ~B · F − ( i~k ˆ U n )) + C ˆ U n (3.8)8ere the discretization on the spatial domain is N × N and ˆ U n is the frequencyof u at time step n . Pseudo-spectral method has higher order accuracy thanthe finite element methods. However, its time step size ∆ t depends on thesmall diffusion parameter ǫ . Therefore it is more costly to reach the large timesolution. The pseudo-spectral method is not so efficient when diffusion parameter ǫ goesto zero. Here we compare the upwinding finite element method and pseudo-spectral method. The analysis suggests a hybrid algorithm combining finiteelement method and pseudo-spectral method. The stability condition of EAFE method is equivalent to that of an upwindingscheme. To simplify analysis, the advection term is taken as one-dimensionaland the reaction term is ignored. The proto-type equation is u t = ǫ u xx + b u x (3.9)The upwinding semi-explicit scheme is ( U n +1 j − U nj ∆ t = ǫ ( U nj +1 + U nj − − U nj + U n +1 j +1 + U n +1 j − − U n +1 j ) / ∆ x + b U nj +1 − U nj ∆ x if b ≥ U n +1 j − U nj ∆ t = ǫ ( U nj +1 + U nj − − U nj + U n +1 j +1 + U n +1 j − − U n +1 j ) / ∆ x + b U nj − U nj − ∆ x if b < b + = max( b, , b − = min( b, , we carry out the von Neumann analysis to obtain the stability condition. Uponsubstitution U nj = ρ n e ijθ , (3.10) becomes: ρ − t = ǫ ( ρ + 1)(cos( θ ) − x + | b | (cos( θ ) − i sin( θ ))∆ x (3.11)Let µ = ǫ ∆ t ∆ x and CFL number C = | b | ∆ t ∆ x . The Peclet number P e = advectiondiffusion = | b | ∆ xǫ = C/µ.
Equation (3.11) simplifies to:(1 + µ (1 − cos θ )) ρ = (1 − µ (1 − cos θ )) + P e µ (cos θ − i sin θ )= (1 − µ (1 + P e )(1 − cos θ )) + i P e µ sin θ (3.12)The stable condition is: ρ ≤ ⇐⇒ µ ≤ − cos θ )(2 + P e )(1 +
P e ) (1 − cos θ ) + ( P e sin θ ) (3.13)9t suffices to impose: µ ≤ min θ ∈ [ − π,π ] 2(1 − cos θ )(2+ P e )(1+
P e ) (1 − cos θ ) +( P e sin θ ) = (2 + P e ) /P e → /P e while ǫ → ǫ ↓ t ≤ ∆ x | b | , (3.14)a CFL condition independent of ǫ . Consider pseudo-spectral method on the proto-type equation (3.9) again: wehave shown before: U n +1 − U n △ t = − ǫ~k U n +1 + b F − ( i~k F U n ) . (3.15)To prove the stability of pseudo-spectral method with Semi-Explicit scheme,we assume: u ( x, t ) = ∞ X k = −∞ ˆ u k ( t ) e ikx . (3.16)and we denote the projection operator:( P N ) u ( x, t ) = X − N
1, the time step is inthe order of O ( ǫ ) which is very restrictive. For each given λ , we compute the Lyapunov exponent of (1.4) by finding the longtime quasi-stationary solution u and applying formula (1.6). A quasi-stationarysolution with λ = 1 is shown in the left plot of Fig. 6. The layer structure of thenumerical solution moves along the y = x direction and oscillates between thetwo neighboring invariant manifolds (lines of unit slope). The right plot of Fig. 6shows the power spectrum of the auxiliary field (a genearalized eigenfunction) w (right) for the unsteady flow (2.2). We observe that the energy decreasestowards high frequency, with a linear decay (scaling range) in the intermediateregion of wave numbers.After calculating µ by the algorithm discussed in (1.6), the minimizationproblem (1.3) is solved by the golden-section method [4]. Golden-section methodis a method to find the minimal point of a particular function called unimodalfunction. The idea of Golden-section method is by successively narrowing downthe range of searching interval. A common definition of unimodal function is forsome value m , it is monotonically decreasing when x < m and monotonicallyincreasing when x ≥ m . In [18], the authors show that the function µ ( λ ) λ is indeeda unimodal function which is strictly convex with respect to λ , decreasing over11nterval [0 , λ ∗ ] and strictly increasing over interval [ λ ∗ , ∞ ], with λ ∗ the uniqueminimal point of µ ( λ ) λ . Notice that Newton’s method is not applicable since µ ′ ( λ ) is not known.For time periodic cellular flow (1.1), EAFE method and spectral method arecombined together to obtain the result in Fig. 7. The initial searching intervalis [0 , λ ∗ with length 2. The space size and time step is h = 1 / and dt = 1 / , respectively. Since the minimal diffusion parameter ǫ = O (10 − )and amplitude of advection term is O (10 ), the choice of dt = 1 / leads tonumerical method stable. The spectral method with space node N = 2 ineach direction and time step dt = 1 / is used to obtain the minimal pointwhen the length of searching interval is less than 2. The stopping criterion ofsearching algorithm is that the difference of two consecutive value of µ ( λ ) λ is lessthan 5 × − . The benefit of using two methods together is that we can findthe minimal point faster without the sacrifice of the accuracy.In Fig. 7, we plot c ∗ as a function of ǫ for three values of parameter θ =10 − , − ,
1. Recall that θ is the perturbation parameter which determines howclose the flow is to the steady cellular flow; see (1.1). From Fig. 7, we observethat for all three positive values of θ , c ∗ ( ǫ ) = O (1) as ǫ ↓
0. In contrast, it is wellknown that c ∗ ( ǫ ) = O ( ǫ / ) in steady cellular flows [20]. The presence of chaotictrajectories in time periodic cellular flows contributed to this phenomenon. Theyare much more mobile and far reaching than their counterparts of steady cellularflows. Fig. 8 shows the persistence of residual front speeds under grid refinementat θ = 1 and the close proximity of the data points between those at grid sizesof N = 2 and N = 2 validates the convergence of the results numerically. We have studied the KPP front speed asymptotics computationally in timeperiodic cellular flows with chaotic streamlines in the small molecular diffusivitylimit. The chaotic streamlines statistically resemble a sub-diffusion process andenhance the KPP front speed c ∗ significantly in the sense that c ∗ has a positivelimit as molecular diffusivity tends to zero. Such residual transport phenomenonis absent in steady cellular flows with ordered streamlines. To facilitate effectivecomputation in the advection dominated regime, we combined an upwindingfinite element methods with the spectral method for computing the principaleigenvalue of a time periodic parabolic operator with small diffusion and for asubsequent minimization. In future work, we plan to study residual KPP frontspeeds in three space dimensional chaotic flows. The work was partially supported by NSF grant DMS-1211179. We thankTyler McMillen for helpful conversations and visualization of chaotic dynamical12 c * theta=1theta=0.1theta=0.01 Figure 7: Residual KPP front speeds in the small diffusion limit ǫ ↓ * vs epsilon when theta=1 c * space node N = 2 space node N = 2 space node N = 2 Figure 8: Residual KPP front speeds at θ = 1 under grid refinement. References [1] E. Abraham, et al.,
Importance of stirring in the development of anironfertilized phytoplanton bloom , Nature 407 (2000), pp. 727–730.[2] M. Abel, M. Cencini, D. Vergni, A. Vulpiani,
Front speed enhance-ment in cellular flows , Chaos 12 (2002), pp. 481–488.[3] B. Audoly, H. Berestycki, Y. Pomeau,
R´eaction diffusion en´ecoulement stationnaire rapide , C.R. Acad. Sci. Paris, S´erie IIb, 328(2000), pp. 255–262.[4] R. Brent, “Algorithms for Minimization WithoutDerivatives”, Prentice Hall, 1973, DFMIN modulehttp://gams.nist.gov/serve.cgi/Module/NMS/DFMIN/5671/.[5] R. Bank, J. Burger, W. Fichtner, R. Smith,
Some upwinding tech-niques for finite element approximations of convection diffusion equa-tions , Numer. Math., 58 (1990), pp. 185–202.[6] L. Biferale, A. Cristini, M. Vergassola, A. Vulpiani,
Eddy diffusivitiesin scalar transport , Physics Fluids 7(11), 1995, pp. 2725–2734.147] A. Brooks, T. Hughes,
Streamline upwind/Petrov-Galerkin formu-lations for convection dominated flows with particular emphasis onthe incompressible Navier-Stokes equations.
Computer methods in ap-plied mechanics and engineering 32(1), 1982, pp. 199-259.[8] R. Camassa, S. Wiggins,
Chaotic advection in a Rayleigh-B´enardflow,
Physical Review A, 43(2), 1990, pp. 774–797.[9] S. Childress, A. Soward,
Scalar transport and alpha-effect for a familyof cat’s eye flows , J. Fluid Mech. 205 (1989), pp. 99–133.[10] J. Cooper, J. Butcher,
An Iterative Scheme for Implicit Runge-KuttaMethods,
IMA J. Numer. Anal. 3, pp. 127–140, 1983.[11] H. Elman, D. Silvester, A. Wathen, “Finite Elements and Fast Iter-ative Solvers with Applications in Incompressible Fluid Dynamics”,Oxford University Press, 2005.[12] A. Fannjiang, G. Papanicolaou,
Convection enhanced diffusion forperiodic flows,
SIAM J. Applied Math, 54(2), 1994, pp. 333-408.[13] H. Kim, J. Xu, and L. Zikatanov,
A multigrid method basedon graph matching for convection-diffusion equations,
Numer. Lin-ear Algebra Appl., (2003), pp. 181 – 195. More details are at:http://epubs.siam.org/doi/ref/10.1137/060659545[14] Y. Gorb, D. Nam, and A. Novikov,
Numerical simulations of diffusionin cellular flows at high Peclet numbers , DCDS-B, vol. 51(1), 2011,pp. 75–92.[15] S. Heinze,
Diffusion-advection in cellular flows with large Peclet num-bers,
Archive Rational Mech. Analysis, 168(4), 2003, pp. 329–342.[16] Y. Liu, J. Xin, Y. Yu,
A Numerical Study of Turbulent Flame Speedsof Curvature and Strain G-equations in Cellular Flows , Physica D,243(1), pp. 20-31, 2013.[17] J. Nolen, J. Xin,
Asymptotic Spreading of KPP Reactive Fronts inIncompressible Space-Time Random Flows , Ann Inst. H. Poincar´e,Analyse Non Lineaire, 26(3), 2009, pp. 815-839.[18] J. Nolen, M. Rudd, J. Xin,
Existence of KPP fronts in spatially-temporally periodic advection and variational principle for propaga-tion speeds , Dynamics of PDE, 2(1), pp. 1-24, 2005.[19] Y. Notay, AGMG software and documentation; seehttp://homepages.ulb.ac.be/ ∼ ynotay/AGMG[20] A. Novikov, L. Ryzhik, Boundary layers and KPP fronts in a cellularflow , Arch. Ration. Mech. Anal., 184(1), 2007, pp. 23 – 48.1521] M. Paoletti, T. Solomon,
Experimental studies of front propagationand mode-locking in an advection-reaction-diffusion system , Euro-phys. Lett. 69 (2005), pp. 819–825.[22] P. Ronney,
Some Open Issues in Premixed Turbulent Combustion ,Modeling in Combustion Science (J. D. Buckmaster and T. Takeno,Eds.), Lecture Notes In Physics, Vol. 449, Springer-Verlag, Berlin,1995, pp. 3-22.[23] L. Ryzhik, A. Zlatos,
KPP pulsating front speed-up by flows , Com-mun. Math. Sci. 5 (2007), pp. 575–593.[24] W. van Saarloos,
Front propagation into unstable states , Physics Re-ports 386(2), 2003, pp. 29-222.[25] L. Shen, A. Zhou, J. Xin,
Finite Element Computation of KPP FrontSpeeds in Cellular and Cat’s Eye Flows , Journal of Scientific Com-puting, 55(2), 2013, pp. 455–470.[26] P. Tong, T.H.H. Pian, L.L. Bucciarelli,
Mode shapes and frequenciesby finite element method using consistent and lumped masses , Com-put. Struct. 1 (1971), pp. 623-638.[27] F. Williams, “Combustion Theory”, Benjamin-Cummings, MenloPark, 1985.[28] J. Xin, “An Introduction to Fronts in Random Media”, Surveys andTutorials in Applied Math Sciences, Vol. 5, Springer, 2009.[29] J. Xin, Y. Yu,
Front Quenching in G-equation Model Induced byStraining of Cellular Flow , Arch. Rational Mechanics Analysis, toappear.[30] J. Xin, Y. Yu,
Asymptotic growth rates and strong bending of tur-bulent flame speeds of G-equation in steady two dimensional incom-pressible periodic flows , SIAM J. Math Analysis, to appear.[31] J. Xu, L. Zikatanov,