Developing new approaches for the path tubes method

Applied Mathematical Modelling 35 (2011) 285–302 Contents lists available at ScienceDirect Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm Developing new approaches for the path tubes method Nélio Henderson a,*, Marcelo Sampaio b, Luciana Pena c a b c Thermodynamics and Optimization Group (TOG), Instituto Politécnico, Universidade do Estado do Rio de Janeiro, 28601-970 Nova Friburgo, RJ, Brazil Eletronuclear – Eletrobrás Termonuclear S.A, 23948-00 Angra dos Reis, RJ, Brazil Laboratório de Ciências Matemáticas, Universidade Estadual do Norte Fluminense, 28013-602 Campos dos Goytacazes, RJ, Brazil a r t i c l e i n f o Article history: Received 29 April 2008 Received in revised form 23 May 2010 Accepted 4 June 2010 Available online 10 June 2010 Keywords: Advection equations Path tubes method Semi-Lagrangian algorithm a b s t r a c t The path tubes (PT) method for advection equations is a semi-Lagrangian algorithm recently introduced in the literature. PT method is simple, physically intuitive, backward-in-time and mass-conserving algorithm, whose formulation is based on the fundaments of the mechanics of continuous media, in the light of Lagrangian methodology. In present work, we describe new approaches for PT methodology, including an absolutely accurate semi-Lagrangian scheme for 1D problems and a new and attractive discretization for multidimensional problems. In addition, we show that PT method is unconditionally stable. Ó 2010 Elsevier Inc. All rights reserved. 1. Introduction The advection equation can be written as follows: @C þ r  ðC v Þ ¼ 0: @t ð1Þ Here, C = C(x, t) is a scalar function that represents the concentration of a given substance, v = v(x, t) is a specified velocity field, x is the spatial variable and t is the time. Since the advection equations govern the evolution of a variety of physical phenomena, Eq. (1) is of great practical importance in the science and engineering. The velocity field can be an intricate function of the space and time, and the solution cannot be analytically calculated. Therefore, in situations of practical interest, numerical approximations to the advection equation are necessary for solving real problems. However, it is well recognized that some methods can produce numerical diffusion, oscillatory behavior and physically unrealistic concentrations, and/or are subject to the Courant–Friedrichs–Lewy (CFL) condition, which imposes restrictions in the time step, see [1–3]. In order to overcome these problems, in a recent article we introduce the path tubes (PT) method, a simple, physically intuitive, backward-in-time, semi-Lagrangian and mass-conserving algorithm, [4]. Since the 1950s, a variety of articles have appeared in the literature, which address aspects of semi-Lagrangian algorithms in different fields of science and technology. Hence, the present method has a special relation with works of several other authors, among them we can mention, for example, Wiin-Nielsen [5], Trulio [6], Hirt et al. [7], Robert [8], Staniforth and Côté [9], Smolarkiewicz and Pudykiewicz [10], Williamson and Olson [11], Gravel and Staniforth [12] and Chilakapati [13]. In present work, we describe and analyze new approaches for PT method. In fact: (i) for the first time, we develop a 1D integral formulation for PT method which uses the Leibnitz formula, (ii) from the three-point Gauss quadrature rule, we * Corresponding author. Fax: +55 22 25288536. E-mail address: [email protected] (N. Henderson). 0307-904X/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2010.06.004 286 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 present an absolutely accurate semi-Lagrangian scheme for 1D advection equations, (iii) we show that this 1D scheme is unconditionally stable, and (iv) using the change of variable formula and the trapezoidal rule, we deduct a new and attractive multidimensional discretization. Moreover, a variety of numerical results with classical (and extremely difficult) 1D and 2D examples are shown to illustrate the robustness of the new formulations. 2. Basic principles of continuum mechanics The path tubes method is based on the fundaments of the mechanics of continuous media in the light of Lagrangian methodology. In view of this, we introduce some physical concepts of interest related with the basic kinematics [14,15]. In continuum mechanics, a material body is characterized by the possibility of occupying different regions of the Euclidean space at different times. Given t 2 R, a three-dimensional body, for example, is idealized as a connected set in R3 , called the configuration of the body at time t. Thus, in different instants t and ~t, the same body X is represented by two possibly different configurations Xt and X~t . During the motion of the body, the location in space is given by a vector x ¼ xðf; tÞ; ð2Þ where f denotes the position of an arbitrary material point f of X in the configuration related with the instant t, see Fig. 1. The spacetime is identified with the Cartesian product R3  R of the Euclidean space R3 and the time axis R. The curve in spacetime along which an arbitrary material particle f travels is referred to as the path line for the particle f. Thus, the path line of f is determined by following set If ¼ fðx; tÞ 2 R3  R; x ¼ xðf; tÞ 2 Xt ; t 2 Rg: ð3Þ The parametric equations of a particle path are the solutions of the system of differential equations dx ¼ v: dt ð4Þ The required additional condition to solve Eq. (4) may be obtained by choosing a reference configuration, a configuration that the material particle assumed at some particular time ~t, ~: xðf; ~tÞ ¼ x ð5Þ [ ð6Þ A path tube R is a region in spacetime formed by union of path lines relative to all particles f of the given material body X, i.e., R¼ If : f2X 3. PT formulation in 1D 3.1. 1D Lagrangian formulation In 1D, the position of f will be denoted by Fig. 1. Positions of an arbitrary material point f. N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 x ¼ xðf; tÞ; 287 ð7Þ and the path line of f is determined as following If ¼ fðx; tÞ 2 R  R; x ¼ xðf; tÞ 2 Xt ; t 2 Rg; ð8Þ where the material body takes the particular form Xt = [x1 (t),x2(t)]. In this case, the parametric equations of a particle path are the solutions of the differential equation dx ¼ v; dt where ð9Þ v = v(x, t) is a specified velocity field in 1D. The required additional condition to solve Eq. (9) is given by xð~tÞ ¼ ~x: ð10Þ The 1D advection equation can be written as follows, @C @ þ ðC v Þ ¼ 0: @t @x ð11Þ Integrating Eq. (11) along the path tube R¼ [ f2X fðx; ~tÞ 2 R  R; x ¼ xðf; ~tÞ 2 X~t ¼ ½x1 ð~tÞ; x2 ð~tÞ; t 6 ~t 6 t þ Dtg; ð12Þ which is defined in the interval [t, t + Dt],we obtain Z tþDt t x2 ðtÞ Z x1 ðtÞ   @C @ðv CÞ dxdt ¼ 0; þ @t @x ð13Þ Eq. (13) can be written as Z tþDt t "Z x2 ðtÞ x1 ðtÞ # @C dx þ ½v ðx2 ðtÞ; tÞCðx2 ðtÞ; tÞ  v ðx1 ðtÞ; tÞCðx1 ðtÞ; tÞ dt ¼ 0: @t ð14Þ On the other hand, by Leibnitz formula (see [16]) we have Z x2 ðtÞ x1 ðtÞ   Z x2 ðtÞ @C dx2 ðtÞ dx1 ðtÞ d Cðx; tÞdx: dx þ Cðx2 ðtÞ; tÞ  Cðx1 ðtÞ; tÞ ¼ @t dt dt dt x1 ðtÞ ð15Þ Using Eq. (9), the relation in Eq. (15) becomes Z x2 ðtÞ x1 ðtÞ @C d dx þ ½v ðx2 ðtÞ; tÞCðx2 ðtÞ; tÞ  v ðx1 ðtÞ; tÞCðx1 ðtÞ; tÞ ¼ @t dt Z x2 ðtÞ Cðx; tÞdx: x1 ðtÞ ð16Þ From Eqs. (14) and (16), we obtain Z t tþDt " d dt Z x2 ðtÞ x1 ðtÞ # Cðx; tÞdx dt ¼ 0: ð17Þ Thus, by the fundamental theorem of calculus, Eq. (17) can be written as Z x2 ðtþDtÞ x1 ðtþDtÞ Cðx; t þ DtÞdx ¼ Z x2 ðtÞ Cðx; tÞdx: x1 ðtÞ ð18Þ This 1D Lagrangian formulation (Eq. (18)) shows that the mass is conserved along the path tube. 3.2. 1D discretization scheme In order to obtain a discretization scheme, in any instant, we use uniform cell centered grids, as shown in Fig. 2. The interval Ij = [xj1/2,xj+1/2] represents an arbitrary grid cell, whose length is Dx = xj+1/2  xj1/2, for all j. In a backtracking step, we solve ordinary differential equations, as indicated in Eq. (9), subject to final condition xðt þ DtÞ ¼ ~ x. In agreement with the illustration in Fg. 3, using backtracking of the path lines, we obtain the interval [x, x+], the image of the cell Ij = [xj1/2, xj+1/2] in the instant t. Thus, it follows that the points x and x+ denote the solutions obtained by solving the ordinary differential equation dx/dt = v, with final conditions ~ x ¼ xi1=2 and ~ x ¼ xiþ1=2 , respectively. Hence, we can write the 1D Lagrangian formulation (Eq. (18)) for each path tube obtained by discretization, as shown in Fig. 3. Thus, for all j, we have 288 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 Fig. 2. The 1D grid. xjþ1=2 Z xj1=2 Cðx; t þ DtÞdx ¼ Z xþ Cðx; tÞdx: ð19Þ x In each instant t + Dt, consider ðnþ1Þ Cj  1 Dx Z xjþ1=2 xj1=2 Cðx; t þ DtÞdx; ð20Þ the average value of the concentration on the jth cell Ij = [xj1/2,xj+1/2]. From Eqs. (19) and (20), we obtain the explicit scheme ðnþ1Þ Cj ¼ 1 Dx Z xþ Cðx; tÞdx: ð21Þ x This 1D scheme (Eq. (21)) is a numerical adaptation of the Lagrangian formulation summarized initially in Eq. (18), being (locally) a mass-conserving algorithm. Rx In order to calculate C in the time level (n + 1), the schemes in Eq. (21) demands the calculation of integrals xþ Cðx; tÞdx. In the present work, these integrals are approximated by the three-point Gauss quadrature rule, see, for example [17,18]. Thus, we obtain Z xþ x Cðx; tÞdx ffi ðxþ  x Þ ½5ðC a þ C b Þ þ 8C c : 18 ð22Þ In Eq. (22), we have Ca  C(xa, t),Cb  C(xb, t) and Cc  C(xc, t); where Fig. 3. A path tube in R  R. N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 pffiffiffiffiffiffiffiffi xa ¼ ðxþ þ x Þ=2  ½ðxþ  x Þ=2 3=5; pffiffiffiffiffiffiffiffi xb ¼ ðxþ þ x Þ=2  ½ðxþ þ x Þ=2 3=5; 289 ð23Þ ð24Þ xc ¼ ðxþ þ x Þ=2: ð25Þ ðnÞ In each instant t, the only known values are C j , for all j, i.e., the values of C in the center of the cells. Using these concentrations, Ca, Cb and Cc are calculated in the following way: if, in instant t, Iq = [xq1/2,xq+1/2] denotes the cell that contains xa, then Ca is approximated by cubic Lagrange interpolation of the values of the concentration in Iq and in the neighbour cells, Iq1, Iq+1 and Iq+2 (or Iq2). The values Cb and Cc are calculated in similar form. In order to avoid oscillatory behavior, our computational experience shows that this procedure requires a small adjustment. The result of the interpolation cannot exceed the concentration values in the interpolation cells. Thus, when necessary, we replaced the calculated value by the maximum value of the concentration in the interpolation cells. A similar fact happens with the minimum of C in the interpolation cells. Thus, we used the following limiter: if xa 2 Iq = [xq1/2,xq+1/2], then C(xa,t) = Ca is approximated by relations where 8 if P3 ðxa Þ < m > < m; C a ffi P3 ðxa Þ; if m 6 P3 ðxa Þ 6 M > : M; if P3 ðxa Þ > M P3 ðxÞ ¼ qþ2 X ðnÞ Ci i¼q1 ð26Þ qþ2 Y x  xk x i  xk k–i ð27Þ k¼q1 is the interpolation polynomial, and the values m and M are given by ðnÞ ðnÞ ðnÞ m ¼ Min fC q1 ; C ðnÞ q ; C qþ1 ; C qþ2 g; ð28Þ M ¼ Max ð29Þ ðnÞ ðnÞ ðnÞ fC q1 ; C qðnÞ ; C qþ1 ; C qþ2 g: In boundary-initial-value problems, if Iq is a boundary cell, the boundary condition value is included in this modified cubic Lagrange interpolation. This naturally incorporates the boundary conditions, when they exist. If m 6 P3(xa) 6 M, then an error E(xa)  Ca  P3 (xa) occur in the approximation Ca ffi P3(xa). This interpolation error can be estimated as follows. Proposition 1. Assume that each of the partial derivatives @ kC(x,t)/@xk,1 6 k 6 5, exists and is continuous in [xq1,xq+2]. Then for all point xa 2 Iq = [xq1/2, xq+1/2] there exists a number na 2 [xq1,xq+2] such that jEðxa Þj 6   3ðDxÞ4 @ 5 Cðna ; tÞ :  64  @x5  ð30Þ Proof. Let x(x) = (xq1  x)(xq  x)(xq+1  x)(xq+2  x). It is well-known that, in these conditions, see [18], there exists na 2 [xq1,xq+2] satisfying Eðxa Þ  C a  P3 ðxa Þ ¼ xðxa Þ @ 5 Cðna ; tÞ 120 @x5 ð31Þ : Since jxðxa Þj 6 45 ðDxÞ4 ; 8 ð32Þ then, from Eqs. (31) and (32), follows the estimate in Eq. (30). As we can see, Eq. (30) shows that jE(xa)j ? 0, whenever Dx ? 0. h 4. Stability analysis In the present work, the stability analysis for PT method is performed in 1D, in agreement with the scheme shown in Eq. (21). To proceed with the stability analysis, we consider a sufficiently fine grid, such that the scheme in Eq. (21) can be approximated by ðnþ1Þ Cj ffi 1 Cð^x; tÞ Dx Z xþ x dx; ð33Þ 290 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 where ^ x is a point in [x, x+] obtained by backtracking of the center of Ij = [xj1/2, xj+1/2]. We also suppose that the flow is incompressible. In this case, the cell Ij = [xj1/2, xj+1/2] is not deformed by backtracking. Thus, xþ Z x dx ¼ Z xjþ1=2 xj1=2 dx ¼ Dx: ð34Þ From Eqs. (33) and (34), we can write ðnþ1Þ Cj ¼ Cð^x; tÞ: ð35Þ Since the flow is incompressible, then @v/@x = 0, i.e., v = constant. In the present wok, we will assume v > 0. Therefore, by backtracking, we have ^x Z xj dx ¼ v Z t dt ð36Þ tþDt or ^x ¼ xj  v Dt: ð37Þ From Eq. (37), we can affirm that there exists an index p P 0 such that xjp1 6 ^x 6 xjp ; ð38Þ where xjp is the center of the cell Ijp = [xjp1/2, xjp+1/2]. In other words ^x ¼ ½j  pDx  aDx; for some a 2 ½0; 1: ð39Þ Thus, we obtain the scheme ðnþ1Þ Cj ¼ C ðnÞ ð½j  pDx  aDxÞ: ð40Þ Since we can have a – 0 and a – 1, then we need an interpolation process. Considering the cubic Lagrange interpolation, we can write ðnþ1Þ Cj ¼ C ðnÞ ð½j  pDx  aDxÞ ¼ jpþ2 X ðnÞ Ck k¼jp1 jpþ2 Y s–k s¼jp1 x  xs : xk  xs ð41Þ After some algebraic effort, we can show that C ðnÞ ð½j  pDx  aDxÞ ¼ jpþ2 X k¼jp1 ðnÞ Ck jpþ2 Y s–k s¼jp1 x  xs ðnÞ ðnÞ ðnÞ ðnÞ ¼ g1 C jp1 þ g2 C jp þ g3 C jpþ1 þ g4 C jpþ2 ; xk  xs ð42Þ i.e., ðnþ1Þ Cj ðnÞ ðnÞ ðnÞ ðnÞ ¼ g1 C jp1 þ g2 C jp þ g3 C jpþ1 þ g4 C jpþ2 ; ð43Þ where g1 ¼ aða þ 1Þða þ 2Þ ; 6 ð1  aÞða þ 1Þða þ 2Þ ; g2 ¼ 2 aða  1Þða þ 2Þ ; g3 ¼ 2 að1  aÞða þ 1Þ g4 ¼ : 6 ð44Þ ð45Þ ð46Þ ð47Þ To study the stability of the scheme in Eq. (43), we consider von Neumann’s method, see [2], for example. Thus, we need investigate a propagation of the term ðnÞ C k ¼ nn eikbDx ; pffiffiffiffiffiffiffi where i ¼ 1; b is the wave number and n is the amplification factor. We can write Eq. (48) in the form ðnÞ C k ¼ nn eik/ : where / = bDx is the phase angle. ð48Þ ð49Þ N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 291 ðnÞ Substitution of C k , as defined in Eq. (49), onto both sides of Eq. (43) shows that n ¼ eip/ ðg1 ei/ þ g2 þ g3 ei/ þ g4 ei2/ Þ: ð50Þ Eq. (50) can be written in the equivalent form, n ¼ eip/ ðA þ iBÞ; ð51Þ A ¼ ðg1 þ g3 Þ cos / þ g4 cos 2/ þ g2 ; B ¼ ðg3  g1 Þ sin / þ g4 sin 2/: ð52Þ ð53Þ where It is not difficult to show the identity   2 A2 þ B2 ¼ g21 þ g22 þ g23 þ g24 þ 2g1 g3 cos2 /  2g1 g3 sin / þ 2½g3 g4 þ g2 ðg1 þ g3 Þ cos / þ 2g2 g4 cos 2/ þ 2g1 g4 cos 3/: ð54Þ Thus,   A2 þ B2 6 g21 þ g22 þ g23 þ g24 þ 2ðg1 g2 þ g1 g3 þ g1 g4 þ g2 g3 þ g2 g4 þ g3 g4 Þ ¼ ðg1 þ g2 þ g3 þ g4 Þ2 : ð55Þ Therefore, we have qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A2 þ B2 6 jðg1 þ g2 þ g3 þ g4 Þj: ð56Þ g1 þ g2 þ g3 þ g4 ¼ 1  a; a 2 ½0; 1: ð57Þ jnj ¼ But, from Eqs. (44)–(47), we obtain Thus, Eqs. (56) and (57) show that PT method satisfies the necessary and sufficient condition for stability, jnj 6 1; for all / and a 2 ½0; 1: ð58Þ We have proved the following result. Proposition 2. For incompressible flows, PT method is unconditionally stable. 5. PT formulation in multidimensions 5.1. The multidimensional Lagrangian formulation The following proposition is the key result used in the determination of the Lagrangian formulation for PT method in multidimensions [4]. Proposition 3 (Reynolds transport theorem). Let Xt be a material volume representing the configuration of a body at some instant t, which moves with velocity v = v(x,t). If W = W (x,t) is any continuously differentiable function of time and position, then Z  Xt  Z @W d Wdx; þ r  ðWv Þ dx ¼ dt Xt @t ð59Þ where d/dt denotes the Lagrangian derivative (or material derivative). Proof. See [14]. Integrating the advection equation (Eq. (1)) along the multidimensional path tube R¼ [ f2X fðx; ~tÞ 2 Rn  R; where n ¼ 2; 3; x ¼ xðf; ~tÞ 2 X~t  Rn and t 6 ~t 6 t þ Dtg; ð60Þ we obtain Z tþDt t Z  Xt  @C þ r  ðC v Þ dxdt ¼ 0: @t ð61Þ Eq. (61) and the Reynolds transport theorem indicate that Z t tþDt  d dt Z Xt Cdx dt ¼ 0: ð62Þ 292 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 By definition of material derivative, see [14], we have d dt Z Xt Cdx ¼ lim Dt!0 "R R Cdx  XtþDt Dt Xt # Cdx ð63Þ : Then, by the fundamental theorem of calculus, we may express Eq. (62) as 0¼ Z tþDt t  d dt Z Xt Cdx dt ¼ Z XtþDt Cdx  Z Cdx: Xt ð64Þ Thus, we obtain the multidimensional Lagrangian formulation Z XtþDt Cðx; t þ DtÞdx ¼ Z Cðx; tÞdx: ð65Þ Xt Eq. (65) shows that mass is conserved also along the multidimensional path tubes. 5.2. A new multidimensional discretization Let jXtþDt j ¼ C ðnþ1Þ  R dx be the measure of the set Xt+Dt, and define XtþDt 1 jXtþDt j Z XtþDt Cðx; t þ DtÞdx ð66Þ as the average value of the concentration on Xt+Dt. From Eqs. (65) and (66), we obtain the multidimensional explicit formula for PT method: C ðnþ1Þ ¼ 1 jXtþDt j Z Cðx; tÞdx: ð67Þ Xt In present section, we introduce a new discretization for Eq. (67). This new scheme for the path tubes method will be developed using the change of variable formula for multiple integrals, a multidimensional integration rule, and an interpolation technique in multidimensions. We consider an uniform grid with blocks centered, where (on two-dimensional spaces, for example) an arbitrary block is described by Xi,j = [xi1/2,xi+1/2]  [yj1/2,yj+1/2], with grid spacing Dx = xi+1/2  xi1/2 and Dy = yj+1/2  yj1/2, and area jXi,jj = DxDy, for all i and j, see Fig. 4. Thus, according to the scheme in Eq. (67), in the level of time n + 1 the average value of the concentration C on the block R Xi,j is C ðnþ1Þ ¼ jX1 j Dt Cðx; tÞdx, i.e., i;j i;j ðnþ1Þ C i;j ¼ 1 D x Dy Z Cðx; y; tÞdxdy; ð68Þ Dt where Dt represent the image set of Xi,j (at some time t) obtained by backtracking of the path lines that start on the time t + Dt. More specifically, we can say that Dt is the image set of the block Xi,j under ut, where ut : Xi;j  R2 ! Dt , given by ut ðx; yÞ ¼ utð1Þ ðx; yÞ; uð2Þ t ðx; yÞ , is a bijective mapping of Xi,j onto Dt, see [4]. This mapping (which determines the backtrack- Fig. 4. The two-dimensional grid. N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 293 ing of the path lines) is obtained by solving the system of differential equations shown in Eq. (4), subject to final conditions as ~. Thus, using the change of variable formula, see [16], we can rewrite the integral in Eq. (68) as xðt þ DtÞ ¼ ~ x and yðt þ DtÞ ¼ y an integral over the block Xi,j, i.e., ðnþ1Þ C i;j ¼ 1 DxDy Z yiþ1=2 yi1=2 Z xiþ1=2 Cðx; y; tÞjJ t ðx; yÞjdxdy; xi1=2 ð69Þ where  ð1Þ  @ u ðx;yÞ  t  @x J t ðx; yÞ ¼  ð2Þ  @ ut ðx;yÞ  @x ð1Þ  @ ut ðx;yÞ   @y  ð2Þ  @ ut ðx;yÞ   @y ð70Þ is the Jacobian determinant of /t. As in Section 3.2, the integral in Eq. (69) can also be approximated by the three-point Gauss quadrature rule. Since, for multidimensional problems, this method requires many integration points, in present work we use the trapezoidal rule, see [18]. Thus, Eq. (69) becomes ðnþ1Þ C i;j ¼ 1 jJ t ðxiþ1=2 ; yjþ1=2 ÞjC ðnÞ ðut ðxiþ1=2 ; yjþ1=2 ÞÞ þ jJ t ðxi1=2 ; yjþ1=2 ÞjC ðnÞ ðut ðxi1=2 ; yjþ1=2 ÞÞ 4 þjJ t ðxiþ1=2 ; yj1=2 ÞjC ðnÞ ðut ðxiþ1=2 ; yj1=2 ÞÞ þ jJ t ðxi1=2 ; yj1=2 ÞjC ðnÞ ðut ðxi1=2 ; yj1=2 ÞÞ : ð71Þ In Eq. (71), the four values of C(n) are obtained by backtracking of the path lines that start on the four vertexes of the block Xi,j = [xi1/2,xi+1/2]  [yj1/2,yj+1/2]. Fig. 5. A point obtained by backtracking is not necessarily the center of a block. 294 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 Here, in Eq. (71), the partial derivatives that constitute the determinants are approximate by finite differences, as shown in the examples below: @ ut ðxiþ1=2 ; yjþ1=2 Þ /t ðxiþ1=2 þ hx ; yjþ1=2 Þ  ut ðxiþ1=2 ; yjþ1=2 Þ ’ ; @x hx @ ut ðxi1=2 ; yj1=2 Þ ut ðxi1=2 ; yj1=2 þ hy Þ  ut ðxi1=2 ; yj1=2 Þ : ’ @y hy ð72Þ ð73Þ From Eqs. (72) and (73) we can see that each of the partial derivative above requires a new backtracking step. To avoid these additional efforts, and to keep the mapping ut constrained to the vertexes of each block, we consider hx = Dx and hy = Dy. Thus, the examples in Eqs. (72) and (73) become, respectively: @ ut ðxiþ1=2 ; yjþ1=2 Þ ut ðxiþ3=2 ; yjþ1=2 Þ  ut ðxiþ1=2 ; yjþ1=2 Þ ’ ; @x Dx @ ut ðxi1=2 ; yj1=2 Þ ut ðxi1=2 ; yjþ1=2 Þ  ut ðxi1=2 ; yj1=2 Þ : ’ Dy @y ð74Þ ð75Þ Now, we have a new numerical semi-Lagrangian explicit scheme. However, since a point obtained by backtracking is not necessarily the center of a block, see Fig. 5, we observe that this scheme is not still complete. To complete it, we consider a multidimensional numerical interpolation. Here, we use a two-dimensional interpolation which has an entirely similar form to 1D modified cubic Lagrangian interpolation given in Eq. (26). Thus, if (xa, yb) 2 [xq1/2, xq+1/2]  [yr1/2, yr+1/2] and (xa, yb) – (xq, yr), then C(xa, yb, t) = Ca,b is approximated by the following interpolation: C a;b where 8 if P3 ðxa ; yb Þ < mx;y > < mx;y ; ffi P3 ðxa ; yb Þ; if mx;y 6 P3 ðxa ; yb Þ 6 M x;y > : Mx;y ; if P3 ðxa ; yb Þ > Mx;y P3 ðx; yÞ ¼ rþ2 X j¼r1 0 B B @ qþ2 X i¼q1 C ðnÞ ðxi ; yj Þ mx;y ¼ Min fC ðnÞ ðxi ; yj Þ; ðnÞ Mx;y ¼ Max fC ðxi ; yj Þ; qþ2 Y k–i k¼q1 ð76Þ 1 r þ2 Y x  xk C y  ys C ; A xi  x k y j  ys s–j ð77Þ s¼r1 8i ¼ q  1; . . . ; q þ 2 and j ¼ r  1; . . . ; r þ 2g; ð78Þ 8i ¼ q  1; . . . ; q þ 2 and j ¼ r  1; . . . ; r þ 2g: ð79Þ 6. Numerical examples in 1D In this section, we present examples involving the 1D transport of concentration-hills by incompressible (Example 1–3) and compressible (Example 4) flows. We consider the differential equation @C @ðC v Þ þ ¼ 0; @t @x ð80Þ with appropriate conditions. In all the examples, we use a grid with 400 points and the path lines are evaluated analytically. 6.1. Example 1: Advancing front (with discontinuity) In this example, we consider a boundary-initial-value problem: Eq. (80) with velocity Cðx; 0Þ ¼ 0; x 2 ½0; 12800 v = 0.5, initial data ð81Þ and subject to boundary conditions Cð0; tÞ ¼ 1; 0 < t 6 tf ; Cð12800; tÞ ¼ 0; 0 < t 6 tf : ð82Þ ð83Þ In Fig. 6, the solutions computed by PT method are shown for several values of Dt. We see that all the numerical solutions to the above problem at tf = 9600 place the jump in approximated the correct place, x = 4800. An analysis of Fig. 6 indicates that our methodology is capable to describe the discontinuity, mainly if Dt is sufficiently large. We can observe that the more accurate numerical solution is reached just when Dt = tf = 9600, where the Courant number (NCo) is 300. All the time steps shown in Fig. 6 are highly prohibitive for methods subject to CFL condition, given by N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 NCo  jv jDt 6 1: Dx 295 ð84Þ 6.2. Example 2: Transport of irregular signal (with discontinuity) Here, we consider a Cauchy problem: Eq. (80) with velocity v = 1 and initial data      2px 2px 1 þ 0:4 sin ; Cðx; 0Þ ¼ 2 þ rðxÞ 1 þ 0:3 sin 9 10 8x 2 ½8; 39; ð85Þ where rðxÞ ¼ 1; 8 6 x 6 28 1; 28 < x 6 39 : Fig. 6. Numerical solution of Example 1 at tf = 9600 for different values of Dt. ð86Þ 296 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 In Fig. 7, the solutions computed by PT method are compared with the exact solutions, for several values of Dt. The numerical solutions are shown as curves with circles. We see that the numerical solutions at tf = 5 are in agreement with the exact solutions, for several values of Dt. In this example, we observe that NCo = 1.29, if Dt = 0.1, and NCo = 64.52, if Dt = tf = 5. Again, the more accurate numerical solution is reached exactly in the maximum value Dt = tf. 6.3. Example 3: Transport of a triangular pulse In this test, we consider Eq. (80) with velocity Cðx; 0Þ ¼ 1  jx  x0 j; jx  x0 j < l0 0; otherwise v = 0.5, x 2 [0,12800] and initial data ; ð87Þ where x0 = 2000 and l0 = 264. Exact and numerical solutions are exhibited in Fig. 8 at tf = 9600, for different values of Dt. As in Example 1, this example uses very large time steps. When Dt = tf = 9600 (and the Courant number is 300) the numerical solution presents better agreement with the exact solution, even in the sharp corners. Fig. 7. Exact (–) and numerical ( ) solutions of Example 2 at tf = 5 for different values of Dt. N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 297 6.4. Example 4: Compressible flow In this case, we simulate the transport of a concentration-hill (with corners) by compressible flow. We consider Eq. (80) with v = x  0.5. The computation was performed in the domain 0 6 x 6 1. As initial data we take the hill of the form Cðx; 0Þ ¼ 0:25½1 þ cosðpkxÞ; ð88Þ k ¼ minðjx  0:5j; 0:2Þ=0:2: ð89Þ where As shown in Fig. 9, the solutions computed by PT method, at tf = 0.25, are in agreement with the exact solutions (even in the corners), for several values of Dt. In this example, we observe that (NCo)max = 0.25, if Dt = 0.00125, and (NCo)max = 50, if Dt = tf = 0.25, where ðNCo Þmax  Dt max jv ðxÞj: Dx 06x61 Fig. 8. Exact (–) and numerical ( ) solutions of Example 3 at tf = 9600 for different values of Dt. ð90Þ 298 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 7. Numerical examples in two dimensions 7.1. Example 5: Rotation of a circular cone Solid body rotation is a classical test for advection algorithms. Generally, it is considered a domain B = [0,Lx]  [0,Ly] and a function f : B  R2 ! R, whose the graph defines the position of a three-dimensional body at some initial instant. The body is rotated (counterclockwise) with angular velocity x, around a fixed axis that passes through the point (x0, y0) 2 B and has direction normal to the xy-plane. Hence, the values that represent the velocities in the x and y directions are u ¼ xðy  y0 Þ; ð91Þ v ¼ xðx  x0 Þ: ð92Þ The partial differential equation that describe this movement is the following advection equation Fig. 9. Exact (–) and numerical ( ) solutions of Example 4 at tf = 0.25 for different values of Dt. 299 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 @C @ðuCÞ @ðv CÞ þ þ ¼ 0; @t @x @y ð93Þ with initial condition given by Cðx; y; 0Þ ¼ f ðx; yÞ: ð94Þ We consider the rotation of a circular cone, where B = [0,2 ]  [0,2], x0 = y0 = 1,x = 1 and 8 pffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi < 3 r  ðx  x Þ2 þ ðy  y Þ2 ; if ðx  x Þ2 þ ðy  y Þ2 6 r 2 c c c c ; f ðx; yÞ ¼ : 0; otherwise ð95Þ with xc = 1, yc = 1.25 and r = 0.70. In this example, the path lines are evaluated analytically. The plots in Fig. 10 show the numerical and exact solutions after two evolutions, with Dt = 4p and Dx = Dy = 0.005. Here, in addition, the accuracy of PT method is measured using the numerical error e2, given by e2 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uP 2 u numer:  C theor: u i;j C i;j i;j ¼u : t 2 P theor: C i;j i;j ð96Þ Fig. 10. Exact and numerical solutions of Example 5. Fig. 11. The error e2 after twenty complete evolutions of the circular cone, with three different values of Dt. 300 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 Fig. 11 shows the error e2 after twenty complete evolutions of the circular cone for three values of Dt, where Nx and Ny are the block numbers in the directions x and y, respectively. We note that, for a given grid, the error increases as the time step decreases, and vice versa. Besides, e2 ? 0, whenever Dx,Dy ? 0. 7.2. Example 6: Deformation of a cylinder Following Leveque [19], we consider the swirling deformation flow: 2 uðx; y; tÞ ¼ sin ðpxÞ sinð2pyÞgðtÞ; v ðx; y; tÞ ¼  sin2 ðpyÞ sinð2pxÞgðtÞ: Fig. 12. Five stages of the swirling deformation flow. ð97Þ ð98Þ N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 301 This flow is defined on the domain B = [0,1]  [0,1]. The function g(t) = cos (pt/tf) is used to introduce time dependence on the interval 0 6 t 6 tf. This flow field reverses direction at the instant that t = tf/2 (instant of larger deformation) in such a way that the initial condition is recovered at final time tf, i.e., C(x,y,tf) = C(x,y,0). Here we use tf = 2. The initial condition for the differential equation described in Eq. (93) is (essentially) a cylinder, given by Cðx; y; 0Þ ¼ ( 1; if ðx  0:3Þ2 þ ðy  0:3Þ2 6 ð0:2Þ2 0; otherwise ð99Þ We consider Dx = Dy = 0.02 and Dt = 0.4. Note that this spatial grid is relatively coarse. In fact, in Example 7.1 we use Dx = Dy = 0.005. In the backtracking step, given v = (u,v) we solve the system of differential equations dx/dt = v, subject to final conditions ~, using an appropriate backward algorithm. Such algorithm combines n steps of the classical xðt þ DtÞ ¼ ~ x and yðt þ DtÞ ¼ y Runge–Kutta method of order 4 with m  n steps of an implicit method, also of order 4, see [20]. This backward algorithm is described below, where we consider m = 80 and n = 3. Since Dt = 0.4, in this case we have h  D t/m = 0.005. Backward Algorithm ~ ¼ yðt þ DtÞ; m P 3; and 3 6 n 6 m: Data:t; Dt; ~x ¼ xðt þ DtÞ; y Step 1. Set h = Dt/m and tm = t + Dt. ~, and x(m) = (x(m),y(m)). x; yðmÞ ¼ y Step 2. Set xðmÞ ¼ ~ Step 3. For j = m util m  (n  1), set:   ðjÞ ðjÞ ðjÞ k1 ¼ v xðjÞ ; yðjÞ ; tj  k11 ; k12 ;  h ðjÞ h ðjÞ h ðjÞ ðjÞ ðjÞ  k21 ; k22 ; k2 ¼ v xðjÞ  k11 ; yðjÞ  k12 ; t j  2 2 2  h ðjÞ h ðjÞ h ðjÞ ðjÞ ðjÞ  k31 ; k32 ; k3 ¼ v xðjÞ  k21 ; yðjÞ  k22 ; t j  2 2 2 ðjÞ ðjÞ ðjÞ ðjÞ ðjÞ k4 ¼ v ðxðjÞ  hk31 ; yðjÞ  hk32 ; tj  hÞ  k41 ; k42 ;   h ðjÞ ðjÞ ðjÞ ðjÞ xðj1Þ ¼ xðjÞ  ðk1 þ 2k2 þ 2k3 þ k4 Þ  xðj1Þ ; yðj1Þ ; 6 t j1 ¼ t j  h: Step 4. For j = m  n util 1, To obtain x(j1) and y(j1), solve the following nonlinear system: ( 1 xðj1Þ ¼ 25 ½xðjÞ  36xðjþ1Þ þ 16xðjþ2Þ  3xðjþ3Þ  12huðxðj1Þ ; yðj1Þ ; t ðj1Þ Þ 1 yðj1Þ ¼ 25 ½yðjÞ  36yðjþ1Þ þ 16yðjþ2Þ  3yðjþ3Þ  12hv ðxðj1Þ ; yðj1Þ ; t ðj1Þ Þ Set xðj1Þ ¼ ðxðj1Þ ; yðj1Þ Þ Set t j1 ¼ tj  h: Five steps of this flow obtained by the present methodology are shown in Fig. 12. The numerical solutions show that the flow indeed reverses directions. Besides, at final time t = 2, the path tubes method is capable to recover (approximately) the initial condition indicated in Eq. (99). 8. Conclusion We have proposed new approaches for the path tubes (PT) method, a semi-Lagrangian metodology for advection equations that was presented in a recent article by Henderson et al. [4]. Summarizing, in this work, the following objectives have been achieved: (a) it was developed a 1D integral formulation for PT method, (b) this formulation led to an accurate discretization scheme for one-dimentional advection equations, which uses the three-point Gauss quadrature rule and a modified cubic Lagrangian interpolation with flux limiters, (c) using the change of variable formula and the trapezoidal rule, it was deduced a new and attractive multidimensional discretization scheme, which works also with a modified multidimensional cubic Lagrangian interpolation with flux limiters, (d) considering an 1D scheme, it was shown that PT method is unconditionally stable, (e) using both 1D and 2D problems, the new schemes were tested on a classical numerical benchmark with difficult examples. We would like to emphasize that PT method is an explicit scheme, which is not subject to CFL condition. Thus, different from many conventional schemes, PT method can work with time steps that are considered as of great magnitude. Besides, PT scheme is a (locally) conservative formulation, whose behavior of the generated solution is non-oscillatory and nondiffusive. 302 N. Henderson et al. / Applied Mathematical Modelling 35 (2011) 285–302 In spite of none three-dimensional problem to have been solved here, it should be observed that the application of the multidimensional scheme shown in Section 5.2 to 3D problems is an immediate consequence. Such applications will be considered in our future works. Acknowledgements The research by N. Henderson has been carried out in the framework of project PROCIENCIA-UERJ financed by FAPERJ. Partial support for this work was also provided to M. Sampaio by FAPERJ. N. Henderson gratefully acknowledges the research grant provided by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico, Ministry of Science & Technology, Brazil). References [1] R.J. Leveque, Finite volume methods for hyperbolic problems, Cambridge Texts in Applied Mathematics, Cambridge University Press, 2002. [2] J.W. Thomas, Numerical partial differential equations: finite difference methods, Texts in Applied Mathematics, vol. 22, Springer, New York, 1995. [3] J.W. Thomas, Numerical partial differential equations: conservation laws and elliptic equations, Texts in Applied Mathematics, vol. 33, Springer, New York, 1999. [4] N. Henderson, M. Sampaio, L. Pena, Path tubes method: a semi-Lagrangian approach for linear advection equations, Chem. Eng. Sci. 64 (2009) 3138. [5] A. Wiin-Nielsen, On the applications of trajectory methods in numerical forecasting, Tellus 11 (1959) 180. [6] J.G. Trulio, Report AFWL-TR-66-19, Air Force Weapons Laboratory, Kirtland Air Force, USA, 1965. [7] C.W. Hirt, J.L. Cook, T.D. Butler, A Lagrangian method for calculation the dynamics of an incompressible fluid with free surface, J. Comput. Phys. 5 (1970) 103. [8] A. Robert, A semi-Lagrangian semi-implicit numerical integration scheme for the primitive meteorological equations, J. Meteorolog. Soc. Jpn. 60 (1982) 319. [9] A. Staniforth, J. Côté, Semi-Lagrangian integration schemes for atmospheric models – a review, Monthly Weather Rev. 119 (1991) 2206. [10] P.K. Smolarkiewicz, J.A. Pudykiewicz, A class of semi-Lagrangian approximations for fluids, J. Atmos. Sci. 49 (1992) 2082. [11] D.L. Williamson, J.G. Olson, A comparison of semi-Lagrangian and Eulerian polar climate simulation, Monthly Weather Rev. 126 (1998) 991. [12] S. Gravel, A. Staniforth, A mass conservative semi-Lagrangian scheme for the shallow wave equation, Monthly Weather Rev. 122 (1994) 243. [13] A. Chilakapati, A characteristic-conservative model for Darcian advection, Adv. Water Resour. 22 (1999) 597. [14] M.E. Gurtin, An Introduction to Continuun Mechanics, Academic Press, New York, 1981. [15] C. Truesdell, R.A. Toupin, Classical field theories of mechanics, Handbuch der Physik III/1, Springer, Berlin, 1960. [16] R.G. Bartle, The Elements of Real Analysis, secon ed., Wiley, New York, 1976. [17] M. Abramowitz, I.M. Stegun, Handbook of Mathematical Functions, Dover, New York, 1970. [18] J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, Springer-Verlag, New York, 1980. [19] R.J. Leveque, High-resolution conservative algorithms for advection in incompressible flow, SIAM J. Numer. Anal. 33 (1996) 627. [20] W.E. Boyce, R.C. DiPrima, Elementary Differential Equations and Boundary Value Problems, eighth ed., John Wiley, New York, 2005.