Finite element methods of an operator splitting applied to population ...

March 2, 2011

16:10

Finite element methods of an operator splitting applied to population balance equations

Naveed Ahmed Institut f¨ ur Analysis und Numerik, Otto-von-Guericke-Universit¨ at Magdeburg, Postfach 4120, D-39016 Magdeburg, Germany. [email protected] Gunar Matthies Universit¨ at Kassel, Fachbereich 10 Mathematik und Naturwissenschaften, Institut f¨ ur Mathematik, Heinrich-Plett-Straße 40, 34132 Kassel, Germany. [email protected] Lutz Tobiska Institut f¨ ur Analysis und Numerik, Otto-von-Guericke-Universit¨ at Magdeburg, Postfach 4120, D-39016 Magdeburg, Germany. [email protected]

In population balance equations, the distribution of the entities depends not only on space and time but also on their own properties referred to as internal coordinates. The operator splitting method is used to transform the whole time-dependent problem into two unsteady subproblems of a smaller complexity. The first subproblem is a timedependent convection-diffusion problem while the second one is a transient transport problem with pure advection. We use the backward Euler method to discretize the subproblems in time. Since the first problem is convection-dominated, the local projection method is applied as stabilization in space. The transport problem in the one-dimensional internal coordinate is solved by a discontinuous Galerkin method. The unconditional stability of the method will be presented. Optimal error estimates are given. Numerical results confirm the theoretical predictions. Keywords: Operator splitting; discontinuous Galerkin; stabilized finite elements; population balance equations. AMS Subject Classification: 65M12, 65M15, 65M60

1. Introduction In this paper, we advocate the operator splitting method to approximate the solutions of population balance equations (PBE). This type of problems arises e.g. from models in the simulation of industrial crystallization process.20 In PBE, the distribution of entities depends not only on space and time but also on their own properties referred to as internal coordinates and the source term usually involves an integral operator. For efficient methods to handle integral operators we refer to 1

March 2, 2011

2

16:10

N. Ahmed, G. Matthies, L. Tobiska

Ref. 29. In this work we consider the source as a known function but still we have a high dimensional system of equations which is challenging from the computational point of view. In order to overcome the difficulty of higher dimensions, the operator splitting method is used to decompose the original problem into two unsteady subproblems of smaller complexity. The first subproblem is a time-dependent convection-diffusion problem while the second one is a transport problem with pure advection. Operator splitting methods are widely used for time integration of unsteady problems. The basic theory of operator splitting for one-dimensional problems can be found in Ref. 41, 45. The concept of operator splitting for time-dependent problems is to decompose the spatial operator into a sum of two or more operators. For example in Ref. 33, the decomposition of convection-diffusion-reaction problem into pure convection and diffusion-reaction problems was studied. For more details about operator splitting methods for linear and non-linear convection-diffusion problems, see Ref. 22, 23, 24, 25, 30. A detailed analysis of an alternating direction implicit (or operator-splitting) scheme is demonstrated in Ref. 26 for the Fokker-Planck equation. The basic idea in Ref. 26 is to split the high dimensional problem into two low dimensional problems corresponding to the configuration and the physical spaces. The solution of the convection-diffusion type problem in configuration space is obtained by a Galerkin spectral method at each quadrature point corresponding to the physical domain. Furthermore, a type of L2 projection is used to update the right-hand side vector at the second stage where the solution of advection equation in physical space is approximated by a finite element method. The main advantage of such splitting is that each subproblem can be discretized and stabilized separately by the best suitable method independently of the other subproblem(s). For example in Ref. 13 the Streamline-Upwind Petrov-Galerkin method (SUPG) has been combined with the standard Galerkin method. The main disadvantage of SUPG scheme, in particular for unsteady problems, is that several terms which include the time derivative, the source term, and second order derivatives have to be added into the stabilizing term in order to ensure the consistency of the resulting method. There are several other stabilization techniques as alternative to SUPG. We mention the local projection stabilization (LPS),3,4,34 the continuous interior penalty method (CIP),5,6,7 the subgrid scale modeling (SGS),15,32 and the orthogonal subscales method (OSS).10,11 The LPS method was originally proposed for the Stokes problem in Ref. 2 as a two-level approach and extended to transport problems in Ref. 3. An analysis of the local projection stabilization applied to Oseen problems can be found in Ref. 4, 34 and for scalar convection-diffusion problems with mixed boundary condition in Ref. 35. The LPS in space combined with a discontinuous Galerkin (dG) method in time for transient convection-diffusion-reaction equations was studied in Ref. 1. A comparison of one- and two-level approaches of local projection stabilization for linear advection-diffusion-reaction problem is presented in

March 2, 2011

16:10

Operator splitting method for population balance equation

3

Ref. 27. For more details about local projection stabilization we refer to Ref. 43 where an overview on recent development for this class of stabilization method applied to scalar convection-diffusion, Stokes and linearized Navier-Stokes problems is given. In this article, we will concentrate on the one-level local projection stabilization technique. The second subproblem in our splitting method is a transport problem with pure advection, so one suitable choice is to approximate it by the discontinuous Galerkin (dG) method. The dG method was first introduced for the neutron transport problem in Ref. 39 and then analyzed in Ref. 31. The theoretical analysis of the dG method for scalar hyperbolic equations can be found in Ref. 21 and for the space-time dG method in Ref. 12. For an introduction to discontinuous Galerkin method we refer to Ref. 9. The application of dG method makes the mass matrix corresponding to the internal coordinate diagonal which leads to the feasibility of parallel implementation without any projection steps between the two sub-steps during the computing process. The aim of the present paper is to combine the local projection stabilization method in space with discontinuous Galerkin method in internal coordinate. We will give the stability and convergence estimates for fully discrete two-step scheme based on an operator splitting method. The format of the paper is as follows: Section 2 introduces the model problem under consideration and defines basic notations. In Section 3, the operator splitting technique is applied to decompose the problem into two simpler ones. We shall formulate the backward Euler discretization and derive the weak form of the two subproblems. Further, we derive the unconditional stability of the two-step method. We then discretize the subproblems in space and internal coordinate using local projection stabilization and discontinuous Galerkin methods, respectively, in Section 4. We show the unconditional stability of the fully discrete two-step method. Section 5 presents the error analysis of the fully discrete scheme. Some implementation issues of the method are discussed in Section 6. Finally, we present in Section 7 some computational results supporting our theoretical results.

2. Model problem Let Ωx be a domain in Rd (d = 2 or 3) with boundary ∂Ωx , Ω` = [`min , `max ] ⊂ R and T > 0. The state of individual particle in population balance equation may consists of external coordinate x, referring to its position in the physical space, and internal coordinate `, representing the properties of particles, such as size, temperature, volume etc. A population balance for a solid process such as crystallization with one internal coordinate can be described by the following partial differential equation:

March 2, 2011

4

16:10

N. Ahmed, G. Matthies, L. Tobiska

Find z : (0, T ) × Ω` × Ωx → R such that  ∂(Gz) ∂z   + − ε∆x z + b(x) · ∇x z = f   ∂`   ∂t z(0, ·) = z0   z|`min = zmin     z=0

in (0, T ] × Ω` × Ωx , in Ω` × Ωx ,

(2.1)

on (0, T ] × Ωx , on (0, T ] × Ω` × ∂Ωx ,

where the diffusion coefficient ε > 0 is a given constant, ∆x and ∇x represent the Laplacian and gradient with respect to x, respectively, b is a given velocity field satisfying ∇x · b = 0, and f is a source function. Here G > 0 represents the growth rate of the particles that depends on ` but is independent of x and t, we also assume that ∂` G ≥ 0, see Ref. 36, 37. Furthermore, let us consider the data of the problem G, b, f , z0 and zmin to be sufficiently smooth functions. Let us introduce some standard notations. Let H m (Ω) denote the Sobolev space of functions with derivatives up to order m in L2 (Ω). We denote by (·, ·) the inner product in L2 (Ω` × Ωx ) and by k · k0 the corresponding L2 -norm defined by Z (v, w) = vw d`dx and kvk20 = (v, v). Ω` ×Ωx

To distinguish the inner products and the corresponding norms with respect to the internal coordinate and the space variable we need some more notations. For this, let us denote by (·, ·)` and k · kL2 (Ω` ) the L2 -inner product and the associated norm in Ω` , respectively, and by (·, ·)x and k · kL2 (Ωx ) the L2 -inner product and the associated norm in Ωx , i.e., Z
v, w ` = vw d` and kvk2L2 (Ω` ) = (v, v)` , Ω` Z
v, w x = vw dx and kvk2L2 (Ωx ) = (v, v)x . Ωx

The norm in the Sobolev space H m (Ωx ) is defined as  1/2 X kvkm =  kDα vk2L2 (Ω )  x

|α|≤m

where α = (α1 , α2 , · · · , αd ) is a multi-index. We also consider certain Bochner spaces. For this, let X be a Banach space with norm k · kX . Then we define n o C(Ω` ; X) = v : Ω` → X, v continuous , Z n o L2 (Ω` ; X) = v : Ω` → X, kv(`)k2X d` < ∞ , Ω`

n ∂j v H m (Ω` ; X) = v ∈ L2 (Ω` ; X) : j ∈ L2 (Ω` ; X), ∂`

o 1≤j≤m ,

{ss1_e1}

March 2, 2011

16:10

Operator splitting method for population balance equation

5

where the derivatives ∂ j v/∂`j are understood in the sense of distribution on Ω` . For spaces X and Y we use the short notation Y (X) := Y (Ω` ; X). The norms in the above defined spaces are given as follows Z 1/2 2 kvkC(X) = sup kv(`)kX , kvkL2 (X) = kv(`)kX d` , `∈Ω`

kvkH m (X)

Ω`

1/2

 Z =

m j 2 X

∂ v

j d` X Ω` j=0 ∂`

.

3. Operator splitting method The numerical method for solving (2.1) in d + 1 variable is based on an operator splitting with respect to (`, t) and (x, t) in Ω` and Ωx direction, respectively. We consider a uniform partition of the time interval τ = T /N , i.e. tn = τ n, n = 1, . . . , N . Then starting with u(t0 ) = z0 , two subproblems are sequentially solved on the sub-intervals (tn , tn+1 ], n = 0, 1, . . . , N − 1: Given u(tn ) find u ˜ : (tn , tn+1 ] × Ω` × Ωx → R such that  ∂u ˜  + Lx u ˜=f in (tn , tn+1 ] × Ω` × Ωx  ∂t (3.1) {spl1} u ˜=0 on (tn , tn+1 ] × Ω` × ∂Ωx   n+ n u ˜(t ) = u(t ). Find u : (tn , tn+1 ] × Ω` × Ωx → R such that  ∂u  + L` u = 0 in (tn , tn+1 ] × Ω` × Ωx  ∂t u|`min = zmin on (tn , tn+1 ] × Ωx   u(tn+ ) = u ˜(tn+1 ),

(3.2)

{spl2}

where L` z =

∂(Gz) , ∂`

Lx z = −ε∆x z + b · ∇x z.

(3.3) {ss1_e2}

This two-steps operator splitting scheme defines u(tn ), n = 1, . . . , N , as an approximation of z(tn ). In the framework of PBE, the first subproblem (3.1) is a time-dependent convection-diffusion equation posed on Ωx parameterized by the variable ` in internal coordinate and the second subproblem (3.2) is a one-dimensional transport problem on Ω` parameterized by the space variable x. Let us consider the spaces V = H01 (Ωx ) and W = H 1 (Ω` ). We introduce the space n o P = v ∈ L2 (Ω` × Ωx ) : v ∈ L2 (Ωx ; W ) ∩ L2 (Ω` ; V ) . A variational form of (3.1) and (3.2) reads as follows:

March 2, 2011

6

16:10

N. Ahmed, G. Matthies, L. Tobiska

First step: Find u ˜ : (tn , tn+1 ] → P with u ˜(tn+ ) = u(tn ) such that Z Z Z

a(˜ u, v) d` = ∀v ∈ P, u ˜t , v x d` + f, v x d` Ω`

Ω`

(3.4) {s2_e1}

Ω`

where the bilinear form a is defined as a(u, v) = ε(∇x u, ∇x v)x + (b · ∇x u, v)x .

{s2_e2}

Second step: Find u : (tn , tn+1 ] → P with u(tn+ ) = u ˜(tn+1 ) such that Z


∀v ∈ P, ut , v x d` + b u, v = (Gz)min , v(`min ) x

(3.5)

Ω`

where wmin = w(`min ) and the bilinear form b is defined as Z    ∂(Gu)  b(u, v) = , v d` + (Gu)(`min ), v(`min ) . ∂` x x Ω` Note that we have imposed the boundary condition (u|`min = zmin ) in `-direction in weak sense. After discretizing in time by the backward Euler method, the first order accurate implicit scheme is considered as two-step method: First step: Given un ∈ P, find u ˜n+1 ∈ P such that Z  n+1 Z Z u ˜ − un  {s2_e2a} , v d` + a(˜ un+1 , v) d` = (f n+1 , v)x d` τ x Ω` Ω` Ω`

(3.6)

for all v ∈ P.

{s2_e5}

Second step: Update u ˜n+1 from the first step and find the solution un+1 ∈ P such that Z  n+1   u −u ˜n+1  n+1 , v d` + b(un+1 , v) = Gmin zmin , v(`min ) (3.7) τ x x Ω` n+1 for all v ∈ P, where zmin = zmin (tn+1 , ·). The next paragraph gives the stability of the two-step method (3.6) and (3.7).

Lemma 3.1 (Stability). Assume that u ˜n , un , n = 1, 2 . . . , N , is the solution obtained from the two-step algorithm (3.6) and (3.7). If ∂` G ≥ 0 and τ ≤ 41 , then the following estimate shows the stability N −1 X

N 2

u + τ 0 n=0

{s2_e6}

Z Ω`

n

2 2ε u ˜n+1 H 1 (Ω

x)

o

2 + ∂` G un+1 L2 (Ω ) d` x

( N −1  X

2

2 1/2

2 ≤ exp(3T /2) u0 + τ 2 f n+1 + G z n+1 2 0

0

n=0

min min

L (Ωx )

) 

.

(3.8)

March 2, 2011

16:10

Operator splitting method for population balance equation

7

Proof. Setting v = u ˜n+1 in (3.6), yields Z Z Z n+1 n n+1 n+1 n+1 (˜ u − u ,u ˜ )x d` + τ a(˜ u ,u ˜ ) d` = τ (f n+1 , u ˜n+1 )x d`. Ω`

Ω` 2

Ω` 2

2

Using the relation 2(a − b)a = a − b + (a − b) , one can write Z 1 n+1 2 1 n 2 1 n+1 u k0 − ku k0 + k˜ u − un k20 . (˜ un+1 − un , u ˜n+1 )x d` = k˜ 2 2 2 Ω` Integrating by parts with respect to x the second term in the bilinear form a(·, ·), one obtains Z Z a(˜ un+1 , u ˜n+1 ) d` = ε k˜ un+1 k2H 1 (Ωx ) d` Ω`

Ω`

n+1

since u ˜ vanishes on the boundary ∂Ωx and ∇x · b = 0. Hence by using CauchySchwarz inequality for the right-hand side, we have for the first step Z k˜ un+1 k20 − kun k20 + k˜ un+1 − un k20 + 2τ ε k˜ un+1 k2H 1 (Ωx ) d` Ω`

≤ τ kf n+1 k20 + τ k˜ un+1 k20 .

(3.9)

{s2_e7}

Substituting v = un+1 in the second step (3.7) gives Z   n+1 (un+1 − u ˜n+1 , un+1 )x d` + τ b(un+1 , un+1 ) = τ Gmin zmin , un+1 (`min ) . (3.10) {s2_e7a} x

Ω`

Starting from n+1

b(u

n+1

,u

Z )= Ω`

 ∂(Gun+1 ) ∂`

, un+1

 x

  d` + Gmin un+1 (`min ), un+1 (`min )

x

an integration by parts twice with respect to ` gives Z

2

2 1 1 n+1 b(un+1 , un+1 ) = (`max ) L2 (Ω ) ∂` G un+1 L2 (Ωx ) d` + G1/2 max u x 2 Ω` 2

2 1 1/2 + Gmin un+1 (`min ) L2 (Ωx ) . 2 where Gmax = G(`max ). Cauchy-Schwarz inequality gives for the right-hand side in (3.10)  

2 1 1/2 n+1 1 1/2 n+1

2 2 Gmin zmin , un+1 (`min ) ≤ Gmin zmin + Gmin un+1 (`min ) L2 (Ω ) . L (Ωx ) x 2 2 x Combining these two results in (3.10) and using the same relation 2(a − b)a = a2 − b2 + (a − b)2 for first term, we get for second step Z

n+1 2 n+1 2 n+1

2 n+1 2

u

− u

˜ + u −u ˜ +τ ∂` G un+1 L2 (Ωx ) 0 0 0 Ω`

1/2 n+1 2

2 ≤ τ Gmin zmin . L (Ωx )

(3.11) {s2_e8}

March 2, 2011

8

16:10

N. Ahmed, G. Matthies, L. Tobiska

Adding (3.9) and (3.11), neglecting some contribution of positive terms, and summing over n = 0, . . . , N − 1, we obtain N −1 Z o n X

n+1 2

N 2 n+1 2

u

2

u + τ

u + ∂ G d` 2ε ˜ ` 1 0 H (Ω ) L (Ω ) n=0

x

Ω`

x

N −1 n N −1 o X X

2

n+1 2 1/2 n+1 2

n+1 2

+ G z

f

2

.

u ≤ u0 0 + τ + τ ˜ min min L (Ωx ) 0 0 n=0

n=0

From (3.9) we have

n+1 2

≤

u ˜ 0

{s2_e8a}

τ

f n+1 2 + 1 un 2 . 0 0 1−τ 1−τ

(3.12)

Using this estimate in the last inequality, we get N −1 Z n o X

2

2

N 2

u + τ 2ε u ˜n+1 H 1 (Ωx ) + ∂` G un+1 L2 (Ωx ) d` 0 n=0

Ω`

N −1 X

2 ≤ u0 0 + τ



n=0



4

f n+1 2 + G1/2 z n+1 2 2 min min L (Ωx ) 0 3

 +

N −1

4τ X

un 2 , 0 3 n=0

where we have used 1/(1 − τ ) ≤ 4/3 for τ ≤ 1/4. We conclude the statement by using Gronwall’s lemma. This completes the proof. The critical issue of the operator splitting method is the overall accuracy of the two-step method. Using Taylor series expansions first order accuracy of the twostep method (3.1) and (3.2) can be shown. A detail error analysis for the first order Lie operator splitting of the sum of two elliptic operators can be found in Ref. 16, 17. Unfortunately, we can’t use these results due to the hyperbolic nature of the operator L` . 4. Fully-discrete method In view of different properties of operator L` and Lx , the operator splitting technique allows us to use different types of discretization methods to solve the problems in Ω` and Ωx . Since the first subproblem (3.5) is convection-dominated, we use the local projection method to stabilize the space discretization. While the second subproblem (3.7) is a transport problem with pure advection, one suitable choice is the discontinuous Galerkin method for the discretization with respect to the internal coordinate. 4.1. Local projection stabilization in space In this subsection, we discretize the subproblem in space. For this, let us denote by {Th } a family of shape regular decompositions of Ωx into d-simplices, quadrilateral or hexahedra. The diameter of a cell K ∈ Th is denoted by hK and h describes

March 2, 2011

16:10

Operator splitting method for population balance equation

9

the maximum diameter of cells K. We will consider the one-level LPS in which the approximation and projection space live on the same mesh. For other variants of LPS we refer to Ref. 18, 28, 34, 40, 43, 44. Let Vh ⊂ V denote the standard finite element space of continuous, piecewise polynomials of degree r. The Galerkin discretization of the problem (3.5) is generally unstable due to dominating advection when the diffusion coefficient is very small ε  1. We will handle this difficulty by adding a stabilizing term based on local projection. Let Dh be the projection space of discontinuous, piecewise polynomials of degree r−1 with r ≥ 1. Let Dh (K) = {qh |K : qh ∈ Dh } be the local projection space and πK : L2 (K) → Dh (K) the local L2 -projection onto Dh (K). Define the global projection πh : L2 (Ωx ) → Dh by (πh v)|K := πK (v|K ). The fluctuation operator κh : L2 (Ωx ) → L2 (Ωx ) is given by κh := id − πh , where id : L2 (Ωx ) → L2 (Ωx ) is the identity mapping. We define the stabilizing term Sh   X Sh (uh , vh ) = µK κh (∇x uh ), κh (∇x vh ) K

K∈Th

with user chosen non-negative constant µK , K ∈ Th . It gives additional control over the fluctuations of gradients. Note that one can also replace the gradient ∇x wh by the derivative in streamline direction b · ∇x wh or (even better Ref. 27, 28) by bK · ∇x wh where bK is a piecewise constant approximation of b, which leads to similar results. The stabilized bilinear form is then defined as ah (uh , vh ) = a(uh , vh ) + Sh (uh , vh ).

(4.1)

{ss3_1e3}

The bilinear form ah is coercive on Vh with respect to the mesh dependent norm  1/2 X |||v||| := ε|v|2H 1 (Ωx ) + µK kκh (∇x v)k2L2 (K) ,

(4.2)

{ss3_1e4}

K∈Th

that is ah (vh , vh ) ≥ |||vh |||2 for all vh ∈ Vh . The stability and convergence properties of the LPS method are based on the following assumptions with respect to the pair (Vh , Dh ), see Ref. 34, 40. Assumption A1 : There is an interpolation operator jh : H 2 (Ω) → Vh such that the approximation properties, kv − jh vk0,K + hK |v − jh v|1,K ≤ Chlk kvkl,K

∀v ∈ H l (Ωx ), 2 ≤ l ≤ r + 1, (4.3)

{ss3_1e6}

for all K ∈ Th and the orthogonality (v − jh v, qh ) = 0 hold true.

∀qh ∈ Dh , ∀v ∈ H 2 (Ω)

(4.4) {ss3_1e7}

March 2, 2011

10

16:10

N. Ahmed, G. Matthies, L. Tobiska

Assumption A2 : The fluctuation operator κh satisfies the following approximation property kκh qk0,K ≤ ChlK |q|l,K

∀K ∈ Th , ∀q ∈ H l (K), 0 ≤ l ≤ r.

(4.5)

In numerical computations, we use mapped finite element spaces, see Ref. 8, where b the enriched spaces are given by on the reference cell K b = Pr (K) b + ˆb4 Pr−1 (K), b Prbubble (K)  b = Qr (K) b + span ˆb x Qbubble (K) ˆr−1 , i = 1, 2 . r i Here, ˆb4 and ˆb are the cubic bubble on the reference triangle and the biquadratic disc bubble on the reference square, respectively. The pairs (Prbubble , Pr−1 ), r ≥ 1, on bubble disc triangles and the pairs (Qr , Pr−1 ), r ≥ 1, on quadrilaterals fulfill assumptions A1 and A2. Further examples of spaces (Vh , Dh ) satisfying A1 and A2 are given in Ref. 34, 40. 4.2. Discontinuous Galerkin method in internal coordinate To discretize (3.5) and (3.7) in internal coordinate `, we apply a discontinuous Galerkin method. Let M > 0 be a given positive integer and `min = `0 < `1 < · · · < `M = `max is a partition of Ω` with Ii = (`i−1 , `i ], ki = `i − `i−1 , and k = max ki . i

Let us introduce the function space of discontinuous piecewise polynomials of degree q ≥ 1 as q n o X Skq = v : Ω` → R : v|Ii (`) = vj `j with vj ∈ R, j = 0, . . . , q . j=0 r,q Then we give the fully discrete space Sh,k as follows

{ss3_1e10}

r,q Sh,k = Vh ⊗ Skq q n X = v : Ω` × Ωx → R : v|Ii (`) = vj `j

with

o vj ∈ Vh , j = 0, . . . , q . (4.6)

j=0

The functions in these spaces are allowed to be discontinuous at the nodes `i , − i = 1, . . . , M − 1. The jumps across the nodes are defined by [φ]i = φ(`+ i ) − φ(`i ), where ± ϕ± m = ϕ(`m ) =

lim

`→`m ±0

ϕ(`).

In the next paragraph, we define the fully discrete scheme based on two-step method.

{ss3_1e11}

r,q r,q First step : For given unh,k ∈ Sh,k , find u ˜n+1 h,k ∈ Sh,k such that Z Z Z  n+1  u ˜h,k − unh,k , X d` + ah (˜ un+1 , X) d` = (f n+1 , X)x d` h,k τ x Ω` Ω` Ω` r,q r,q for all X ∈ Sh,k where u0h,k is a suitable approximation of z0 in Sh,k .

(4.7)

{ss3_1e8}

March 2, 2011

16:10

Operator splitting method for population balance equation

11

r,q n+1 Second step : Update the solution u ˜n+1 h,k from (4.7) and find uh,k ∈ Sh,k such that

{ss3_1e12}

Z  n+1 uh,k − u ˜n+1 h,k τ

Ω`

,X

 x

  n+1 + d` + B(un+1 , X) = G z , X(` ) min min,h 0 h,k

(4.8) x

r,q r,q n+1 n+1 for all X ∈ Sh,k , where zmin,h ∈ Sh,k is an approximation of zmin and the bilinear form B is defined as M −1  M Z     X X 
 ∂(Gu)  + , v d` + Gu i , v(`+ B(u, v) = + Gmin u(`+ . 0 ), v(`0 ) i ) ∂` x x x i=1 i=1 Ii

(4.9) {b_uv_1} Integrating by parts Z  Z     
+
∂v  ∂(Gu)  + − Gu (` ), v(` ) − Gu, ), v(` ) − , v d` = Gu (`− d` i−1 i−1 i i ∂` ∂` x x x x Ii Ii leads to the representation M −1   X 
   ∂v  − + Gmax u(`− . d`− u(`− i ), Gv i M ), v(`M ) ∂` x x x i=1 Ii i=1 (4.10) {b_uv_2} We introduce the mesh dependent norm

B(u, v) = −

kvk2dG

=

M Z  X

M Z X i=1

Ii

Gu,

1/2

2

∂` Gkvk2L2 (Ωx ) d` + Gmin v(`+ 0 ) L2 (Ω

x

+ )

M −1 X

 1/2  2

(G v) 2 i L (Ω

x)

i=1

− 2 + G1/2 max v(`M ) L2 (Ωx ) .

(4.11) {dg_norm}

Lemma 4.1. The bilinear form B is coercive with respect to the mesh dependent norm k · kdG , i.e., B(v, v) ≥

1 kvk2dG . 2

(4.12)

r,q holds for all v ∈ Sh,k .

Proof. Setting u = v in (4.9) and (4.10), then adding them together we conclude the statement of the lemma. The next lemma provides a stability result of the fully discrete two-step method (4.7) and (4.8).

{B_coer}

March 2, 2011

12

16:10

N. Ahmed, G. Matthies, L. Tobiska

Lemma 4.2 (Stability). Let ∂` G ≥ 0 and τ ≤ 1/2, then the solution u ˜nh,k and unh,k , n = 1, 2, . . . , N , of (4.7) and (4.8),respectively, satisfies N −1 Z N −1 X X

N 2 n+1 2

n+1 2

uh,k + 2τ u d` + τ

u

˜ h,k h,k dG 0 n=0

Ω`

n=0

( ) N −1  X

n+1 2 1/2 n+1 2

0 2 4

+ (G z

f

. ≤ exp(3T /2) uh,k 0 + τ min min,h ) L2 (Ωx ) 0 3 n=0 (4.13) Proof. Following the similar derivation steps as in Lemma 3.1, we get the proof of lemma. 5. Error analysis In this section, we derive the error estimates of the fully discrete two-step scheme (4.7) and (4.8). First we define a special interpolant Πk w(t, ·, x) ∈ Skq of a function w(t, `, x) by − Πk w(`− i ) = w(`i ),

{int_pi1}

Z {int_pi2}

(Πk w − w)`s d` = 0,

i = 1, . . . , M − 1, s ≤ q − 1, i ≥ 1,

(5.1) (5.2)

Ii

i.e., Πk w interpolates at the nodal points and the interpolation error is orthogonal to the space of polynomials of degree q − 1 on Ii . For this type of interpolant we have the following error estimates {int_pi3}

sup |Πk w(`) − w(`)|j ≤ Ck q+1 sup |w(q+1) (`)|j , 0≤`≤`M Z Z (i) (i) 2 2(q+1−i) |Πk w (`) − w (`)|j d` ≤ Ck |w(q+1) (`)|j d`,

j = 0, 1,

(5.3)

i, j = 0, 1,

(5.4)

0≤`≤`M

{int_pi4}

Ii

Ii

see Ref. 38, 42. In order to obtain the error estimate for the splitting method in space and internal coordinate, we define a projection operator Rh which maps onto r,q the tensor product space Sh,k . It is defined as follows Rh w = jh Πk w = Πk jh w

{R}

∀w ∈ P,

(5.5)

where jh is the special interpolant in space satisfying Assumption A1. In addition, we have the stability property of interpolant Πk given by Z Z

Πk u 2 r+1 {stb_jpi} d` ≤ C kuk2H r+1 (Ωx ) d` (5.6) H (Ω ) Ω`

x

Ω`

since Πk acts in `-direction and the norms are with respect to the space direction. Let us consider ξ n := u(tn ) − Rh u(tn ) and η n := Rh u(tn ) − unh,k . We also denote ξ˜n := u ˜(tn ) − Rh u ˜(tn ) and η˜n := Rh u ˜(tn ) − u ˜nh,k , then the error u(tn ) − unh,k can be decomposed as follows en = u(tn ) − unh,k = ξ n + η n

{ss3_1e19}

March 2, 2011

16:10

Operator splitting method for population balance equation

13

where unh,k is the solution for fully discrete scheme (4.7) and (4.8) and u(tn ) is the solution of (3.1) and (3.2). Furthermore, to obtain the separate estimates in space and internal coordinate we use the following decomposition of errors

Rh w − w = Rh w − Πk w + Πk w − w = ϑ + ϕ. (5.7) {sp_sp_int} Assumption A3 : Let u, ut , utt , u ˜, u ˜t , u ˜tt , zmin and z0 satisfy the following regularity conditions



u, u ˜ ∈ H 1 L2 (H r+1 ) ∩ H 1 H q+1 (H 1 ) , ut , u ˜t ∈ L2 L2 (H r+1 ) ∩ L2 H q+1 (L2 ) ,


utt , u ˜tt ∈ L2 L2 (L2 ) , z0 ∈ L2 Ω` ; H r+1 (Ωx ) ∩ H q+1 Ω` ; L2 (Ωx ) ,
zmin ∈ H 1 0, T ; H r+1 (Ωx ) .

Lemma 5.1. Let the assumptions A1-A3 be fulfilled. Then for all t ∈ (0, T ], we have the following estimates for the interpolation error   kϑ(t)kdG ≤ C hr+1 ku(t)kL2 (H r+1 ) + ku(t)kC(H r+1 ) , kϕ(t)kdG ≤ C k q+1/2 ku(t)kH q+1 (L2 ) . Proof. For simplicity we skip the dependency of t within the proof. Since for the interpolation error the jumps [jh u − u]i , i = 1, . . . , M − 1, vanishes due to the continuity of jh u in internal coordinate, we have from (4.11), the interpolation error estimates (4.3) and condition (5.6) M −1 Z X 1/2 − 2 2 1/2 2 kϑkdG ≤ ∂` Gkϑk2L2 (Ωx ) d` + kGmin ϑ(`+ 0 )kL2 (Ωx ) + kGmax ϑ(`M )kL2 (Ωx ) i=1

≤Ch

Ii

2r+2



kuk2L2 (H r+1 )

+

kuk2C(H r+1 )

 .

For the second estimate with respect to the internal coordinate, we use the definition of interpolant Πk u, i.e., the interpolation Πk u satisfies Πk u(`− i ) = u(`i ), i = 1, . . . , M , thus from the second representation (4.10) of the bilinear form B and interpolation estimates (5.3), (5.4), we have M Z  X ∂ϕ  d` kϕk2dG ≤ B(ϕ, ϕ) = − Gϕ, ∂` x i=1 Ii M Z X ≤ kGϕkL2 (Ωx ) k∂` ϕkL2 (Ωx ) d` i=1

≤Ck

Ii 2q+1

M Z X i=1

Ii

kuq+1 k2L2 (Ωx ) d` ≤ C k 2q+1 kuk2H q+1 (L2 )

which completes the proof of the lemma.

March 2, 2011

14

16:10

N. Ahmed, G. Matthies, L. Tobiska

Lemma 5.2. Let the assumptions A1-A3 be fulfilled and τK ∼ hK . Then for all t ∈ (0, T ], the following estimates hold Z


ah ϑ(t), η(t) d` ≤ C (ε1/2 + h1/2 )hr ku(t)kL2 (H r+1 )

Ω`

2

Ω`

+Ch Z

1/2 |||η(t)||| d`

Z

r+1

ku(t)kL2 (H r+1 ) kη(t)k0 ,


ah ϕ(t), η(t) d` ≤ C (ε1/2 + h1/2 ) k q+1 ku(t)kH q+1 (H 1 )

Ω`

(5.8) {c4:lem7:atheta} 1/2 2 |||η(t)||| d`

Z Ω`

+Ck

q+1

ku(t)kH q+1 (H 1 ) kη(t)k0 ,

(5.9)



Proof. For simplicity of the presentation we again skip the dependency of the time within the proof. From the definition of the stabilized bilinear form ah , we have Z Z Z Z



ah ϑ, η d` = ε ∇x ϑ, ∇x η x + b · ∇x ϑ, η x d` + Sh ϑ, η d` Ω`

Ω`

Ω`

Ω`

= I1 + I2 + I3 .

{c4:lem7:e1}

(5.12)

We start by estimating the first term on the right-hand side. Using Cauchy-Schwarz inequality, the interpolation estimates (4.3) of jh and the condition (5.6), it follows that Z Z |I1 | ≤ ε ||ϑ||H 1 (Ωx ) ||η||H 1 (Ωx ) d` ≤ C ε1/2 hr kΠk ukH r+1 |||η||| d` Ω`

≤ Cε1/2 hr

Ω`

Z

{c4:lem7:aphi}



B ϑ(t), η(t) ≤ C hr+1 ku(t)kH 1 (H r+1 ) kη(t)k0 + ku(t)kC(H r+1 ) kη(t)kdG , (5.10)
B ϕ(t), η(t) ≤ C k q+1 ku(t)kH q+1 (L2 ) kη(t)k0 . (5.11)


1/2 Z kuk2H r+1 d`

Ω`

≤ C ε1/2 hr kukL2 (H r+1 )

|||η|||2 d`

1/2

Ω`

Z

|||η|||2 d`

1/2 .

Ω`

Integrating I2 by parts with respect to the space variable x, using the orthogonality property of interpolant jh and Cauchy-Schwarz inequality to get Z Z

|I2 | = b · ∇x ϑ, η x d` = ϑ, b · ∇x η x d` Ω` Ω` Z
≤ ϑ, κh (b · ∇x η) x d` Ω` Z X ≤ kϑkL2 (K) kκh (b · ∇x η)kL2 (K) d`. Ω` K∈T h

Let b be the L2 -projection of b in the space of piecewise constant functions with respect to Th . Using the L2 -stability of the fluctuation operator κh , inverse inequality

{c4:lem7:btheta} {c4:lem7:bphi}

March 2, 2011

16:10

Operator splitting method for population balance equation

{C4_th1_14a}

15

and κh (b · ∇x )η = b · κh (∇x η), we get in a same fashion as in Ref. 34 the following estimate



κh (b · ∇x )η 2 ≤ C|b|1,∞,K kηkL2 (K) + kbk0,∞,K κh (∇x η) 2 . (5.13) L (K)

L (K)

Thus inserting this in the previous estimate, using (4.3), µK ∼ hK , and (5.6) to get Z X |I2 | ≤ C |b|1,∞,K kϑkL2 (K) kηkL2 (K) d` Ω` K∈T h

Z +C

X

Ω` K∈T h

|b|0,∞,K kϑkL2 (K) κh (∇x η) L2 (K) d`

Z kukH r+1 (Ωx ) kηkL2 (Ωx ) d` + C hr+1/2 kukH r+1 (Ωx ) |||η||| d` Ω` Ω`  1/2  Z kukL2 (H r+1 ) . ≤ C hr+1/2 h1/2 kηk0 + |||η|||2 d`

≤ C hr+1

Z

Ω`

For I3 , the Cauchy-Schwarz inequality and interpolation error estimates give Z Z

1/2
1/2 Sh ϑ, η d` ≤ Sh ϑ, ϑ Sh η, η d` |I3 | = Ω`

Ω`

r+1/2

Z

r+1/2

≤ Ch

kukH r+1 (Ωx ) |||η||| d` ≤ Ch

1/2 |||η||| d` .

Z

2

kukL2 (H r+1 )

Ω`

Ω`

Combining I1 , I2 and I3 , we get the desired estimate Z Z
ah ϑ, η d` ≤ C (ε1/2 + h1/2 )kukL2 (H r+1 ) Ω`

|||η|||2 d`

1/2

Ω`

+ C hr+1 kukL2 (H r+1 ) kηk0 . Next, we find the estimates in internal coordinate. From the definition, we have Z Z Z Z


Sh ϕ, η d` ah (ϕ, η) d` = ε ∇x ϕ, ∇x η x d` + b · ∇x ϕ, η x d` + Ω`

Ω`

Ω`

Ω`

= I4 + I5 + I6 . Then by using the Cauchy-Schwarz inequality, the stability property of the fluctuation operator κh , the approximation properties (5.3) of interpolant Πk and the parameter choice µK ∼ hK , we get for I4 , I5 , and I6 the following estimates Z |I4 | ≤ ε kΠk u − ukH 1 (Ωx ) kηkH 1 (Ωx ) d` Ω` Z ≤ ε1/2 kΠk u − ukH 1 (Ωx ) |||η||| d` Ω` Z 1/2  Z 1/2 ≤ ε1/2 kΠk u − uk2H 1 (Ωx ) d` |||η|||2 d` Ω` Ω` Z 1/2 ≤ C ε1/2 k q+1 kukH q+1 (H 1 ) |||η|||2 d` . Ω`

|I5 | ≤ C k q+1 kukH q+1 (H 1 ) kηk0 .

March 2, 2011

16

16:10

N. Ahmed, G. Matthies, L. Tobiska

Z

X

|I6 | ≤

Ω` K∈T h 1/2



µK κh (∇x (Πk u − u)) L2 (K) κh (∇x η) L2 (K) d`

Z



∇x (Πk u − u) 2 |||η||| d` L (Ωx ) Ω` Z  1/2 ≤ C h1/2 k q+1 kukH q+1 (H 1 ) |||η|||2 d` .

≤Ch

Ω`

Hence, combining these estimates we get the second statement of the lemma Z Z 1/2
ah ϕ, η d` ≤ C (ε1/2 + h1/2 ) k q+1 kukH q+1 (H 1 ) |||η|||2 d` Ω`

Ω`

+ C k q+1 kukH q+1 (H 1 ) kηk0 . To obtain the last two estimates, we use the two different representations (4.9) and (4.10) of B. Note that the jump terms [jh u − u]i , i = 1, . . . , M − 1, vanishes due to the continuity of the interpolant jh u in `-direction. We have from (4.9), (4.3), and (5.6) M Z    X ∂(Gϑ)  + , η d` + Gmin ϑ(`+ ), η(` ) 0 0 ∂` x x i=1 Ii Z M X

1/2

1/2

+

∂` (Gϑ) 2 kηkL2 (Ωx ) d` + Gmin ϑ(`+ ≤ 0 ) L2 (Ωx ) Gmin η(`0 ) L2 (Ωx ) L (Ωx )

B(ϑ, η) =

i=1

Ii

≤Ch

r+1



 kukH 1 (H r+1 ) kηk0 + kukC(H r+1 ) kηkdG .

The interpolation Πk u satisfies Πk u(`− i ) = u(`i ), i = 1, . . . , M . Hence, we get from (4.10) the relation B(ϕ, η) =

M Z X i=1

 ∂η  − Gϕ, d`. ∂` x Ii

2

Let Π0 G be the L -projection of G in a space of piecewise constant functions in `-direction. Using the orthogonality (5.2) of the interpolant Πk , we get B(ϕ, η) =

M Z  X i=1

≤

Ii

M Z X i=1

Ii

ϕ, (G − Π0 G)

∂η  d` ∂` x

kϕkL2 (Ωx ) (G − Π0 G)∂` η L2 (Ωx ) dx

≤ C k q+1 kukH q+1 (L2 ) kηk0 . Here, we used the Cauchy-Schwarz inequality, the inverse inequality and the interpolation error estimates (5.3). This complete the proof.

March 2, 2011

16:10

Operator splitting method for population balance equation

17

Theorem 5.1. Let u ˜(tn ), u(tn ) and u ˜nh,k , unh,k , be the solutions of two-step method (3.1), (3.2) and (4.7), (4.8), respectively. Under the assumptions A1-A3 and µK ∼ hK there holds for η n = Rh u(tn ) − unh,k and η˜n = Rh u(tn ) − unh,k N −1 X

N 2

η + τ 0

N −1 X n+1 2

n+1 2 d` + τ

η˜

η dG 2 Ω` n=0 i h

2 ≤ Cu e9T /2 Rh z0 − u0h,k 0 + τ 2 + (ε + h) h2r + k 2q+2

Z

n=0

{err_e23}

(5.14)

and for en = u(tn ) − unh,k and e˜n = u ˜(tn ) − u ˜nh,k N −1 X

N 2

e + τ 0

N −1 X n+1 2

n+1 2 e˜ d` + τ

e dG 2 n=0 n=0 Ω`  

2 ≤ Cu e9T /2 Rh z0 − u0h,k 0 + τ 2 + (ε + h) h2r + k 2q+1

Z

(5.15)

{err_e23a}

where Cu depends on u, ut , utt , u ˜, u ˜t , u ˜tt and zmin . Note that the error to the interpolant Rh u is superclose with respect to the internal coordinate (order k + 1 instead of k + 1/2). Proof. From the result of the Lemma 4.2, we can write for η n = Rh u(tn ) − unh,k

1

η N 2 − 0 2

N −1 X

1

η 0 2 + τ 0 2 n=0

N −1 X n+1 2

n+1 2 η˜ d` + τ

η

≤ T1 + T2 (5.16) {err_e24} dG 2 n=0 Ω`

Z

where N −1 Z X

 η˜n+1 − η n n+1  , η˜ + ah (˜ η n+1 , η˜n+1 ) d`, τ x n=0 Ω`   Z  n+1 N −1 X η − η˜n+1 n+1  n+1 n+1 ,η T2 = τ d` + B(η ,η ) . τ x Ω` n=0 T1 = τ



(5.17)

(5.18)

We first consider T1 . Using (4.7), we obtain N −1 Z X


Rh u ˜(tn+1 ) − Rh u(tn ) n+1  , η˜ + ah Rh u ˜(tn+1 ), η˜n+1 τ x n=0 Ω`  Z n  u ˜n+1 h,k − uh,k n+1 , η˜n+1 − ah (˜ un+1 , η ˜ ) d` − h,k τ x Ω`  N −1 Z X
Rh u ˜(tn+1 ) − Rh u(tn ) n+1  =τ , η˜ + ah Rh u ˜(tn+1 ), η˜n+1 τ x n=0 Ω` 
− f n+1 , η˜n+1 x d`.

T1 = τ



March 2, 2011

18

16:10

N. Ahmed, G. Matthies, L. Tobiska

For the last term on the right-hand side of the above equation, using (3.4) at t = tn+1 , we get for the first term N −1 Z X

R u  n+1 ) − Rh u(tn ) h ˜(t −u ˜t (tn+1 ), η˜n+1 d` τ x n=0 Ω` N −1 Z N −1 Z X X

+τ a Rh u ˜(tn+1 ) − u ˜(tn+1 ), η˜n+1 d` + τ Sh Rh u ˜(tn+1 ), η˜n+1 d`

T1 = τ

Ω`

n=0

Ω`

n=0

N −1 Z X

 R u n+1 ) − Rh u(tn ) h ˜(t −u ˜t (tn+1 ), η˜n+1 d` τ x n=0 Ω` N −1 Z N −1 Z X X

+τ ah Rh u ˜(tn+1 ) − u ˜(tn+1 ), η˜n+1 d` + τ Sh u ˜(tn+1 ), η˜n+1 d`

=τ

Ω`

n=0

n=0

Ω`

= T1,1 + T1,2 + T1,3 .

(5.19) {T_1}

We treat the contribution of the terms on the right-hand side of (5.19) separately. For the first term, using Cauchy-Schwarz inequality, the Young’s inequality and the initial condition u ˜(tn ) = u(tn ) for first step N −1 Z X

|T1,1 | ≤ τ

Ω`

n=0

≤

τ 2

n+1

Rh u ) − Rh u(tn ) n+1

˜(t −u ˜t (t )

τ

L2 (Ωx )

N −1 Z X Ω`

n=0

n+1

η˜ kL2 (Ωx ) d`

2

n+1

Rh u ) − Rh u(tn ) n+1

˜(t − u ˜ (t ) t

2

τ

d`

L (Ωx )

N −1 Z τ X k˜ η n+1 k2L2 (Ωx ) d` 2 n=0 Ω`

2 N −1 X

Rh u

˜(tn+1 ) − Rh u ˜(tn ) n+1

≤τ − Rh u ˜t (t )

τ 0 n=0

+

+τ

N −1 X

−1

2 τ NX

n+1 2

Rh u

η˜

. ˜t (tn+1 ) − u ˜t (tn+1 ) 0 + 0 2 n=0 n=0

For first term, applying Taylor’s theorem with integral remainder term and for second term the approximation properties of interpolant jh and Πk with the stability property (5.6) yields |T1,1 | ≤ τ

2

N −1 Z tn+1 X n=0

+ Cτ

tn

N −1  X

N −1

2 τ X

u

η˜n+1 2 ˜tt 0 dt + 0 2 n=0





n+1 2 2q+2 n+1 2 u ˜t (t ) L2 (H r+1 ) + k u ˜t (t ) H q+1 (L2 ) .

2r+2

h

n=0

To find the estimates for T1,2 , we use the decomposition (5.7) of errors into space

March 2, 2011

16:10

Operator splitting method for population balance equation

19

and internal coordinate and get  N −1  X

n+1 n+1 n+1 n+1 ˜ T1,2 = τ ah ϑ , η˜ + ah ϕ˜ , η˜ . n=0

Then from the results (5.8) and (5.9) of Lemma 5.2, we obtain  N −1  X

n+1 2

n+1 2 2r 2q+2

|T1,2 | ≤ C (ε + h) τ ˜(t u ˜(t h u ) L2 (H r+1 ) + k ) H q+1 (H 1 ) n=0

+Cτ

N −1  X

n+1 2

n+1 2 ˜(t ) H q+1 (H 1 ) ˜(t ) L2 (H r+1 ) + k 2q+2 u h2r+2 u



n=0

+

τ 4

N −1 Z X n=0

N −1 X n+1 2

n+1 2 d` + τ

. η˜

η˜ 0 2 n=0 Ω`

The estimate for T1,3 follows from the approximation properties of the fluctuation operator κh and the choice of the stabilizing parameter µK ∼ hK . We have N −1 Z N −1 Z X

τ X n+1 n+1 |T1,3 | ≤ τ Sh u ˜(t ), u ˜(t ) d` + Sh η˜n+1 , η˜n+1 d` 4 n=0 Ω` n=0 Ω` N −1 N −1 Z X

n+1 2 τ X 2r+1 n+1 2

η˜ d`. ≤Ch τ u ˜(t ) L2 (H r+1 ) + 4 n=0 Ω` n=0 Finally, by inserting the estimates T1,1 , T1,2 , and T1,3 into (5.19), we obtain N −1 Z tn+1 N −1 Z N −1 X X

2 n+1 2

n+1 2 τ X 2



η˜

u ˜tt 0 dt + |T1 | ≤ τ η˜ d` + τ 0 2 n=0 Ω` n n=0 t n=0  N −1  X

n+1 2

2 (ε + h) u ˜(t ) 2 r+2 + h2 u ˜t (tn+1 ) 2 r+1 + C h2r τ L (H

)

L (H

)

n=0

+Ck

2q+2

τ

N −1  X



n+1 2

n+1 2

(ε + h + 1) u ˜(t ) H q+1 (H 1 ) + u ˜t (t ) H q+1 (L2 ) .

n=0

(5.20) {T_1_n} Now we estimate the second term T2 . Using (4.8) and (3.5) we obtain the following error equation for second step N −1 Z  R u(tn+1 ) − R u  n+1 X ) h h ˜(t T2 = τ − ut (tn+1 ), η n+1 d` τ x n=0 Ω` +τ

N −1 X

B Rh u(tn+1 ) − u(tn+1 ), η n+1




n=0

−τ

N −1  X

 n+1 n+1 Gmin zmin − Gmin zmin,h , η n+1 (`+ ) 0

n=0

= T2,1 + T2,2 + T2,3 .

x

(5.21) {T_2}

March 2, 2011

20

16:10

N. Ahmed, G. Matthies, L. Tobiska

The estimates for the first term can be obtained by using the same procedure as for T1,1 and get

2 N −1 X

Rh u(tn+1 ) − Rh u

˜(tn+1 ) n+1

|T2,1 | ≤ τ − R u (t ) h t

τ 0 n=0 N −1 X

n+1 2



η

Rh ut (tn+1 ) − ut (tn+1 ) 2 + τ 0 0 2 n=0 n=0 N −1 Z tn+1 N −1 X X

2

n+1 2 2

utt dt + τ

η ≤τ 0 0 2 n t n=0 n=0  N −1  X



2r+2 n+1 2 2q+2 n+1 2 + Cτ ut (t ut (t h ) L2 (H r+1 ) + k ) H q+1 (L2 ) .

+τ

N −1 X

n=0

Note that in above estimates we have used the initial condition u(tn ) = u ˜(tn+1 ) from (3.2). The bounds on the second term T2,2 are obtained by using the error decomposition (5.7) and the estimates (5.10) and (5.11) N −1 Xn o

B ϑn+1 , η n+1 + B ϕn+1 , η n+1 |T2,2 | = τ n=0  N −1  X

2

n+1 2

u(t ) 1 r+1 + u(tn+1 ) ≤ C h2r+2 τ r+1 H (H

C(H

)

)

n=0

+ C k 2q+2 τ

N −1 X

N −1 N −1 X

n+1 2

n+1 2 τ X

u(t

η n+1 2 + τ

η

. ) H q+1 (L2 ) + 0 dG 2 8 n=0 n=0 n=0

Cauchy-Schwarz inequality and Young’s inequality give for T2,3 |T2,3 | ≤ τ

N −1 X

1/2

1/2 n+1 +

G zmin (tn+1 ) − G1/2 z n+1 2

(`0 ) L2 (Ω min min min,h L (Ω ) Gmin η

x)

x

n=0

≤ C h2r+2 τ

N −1

τ X

η n+1 2 .

zmin (tn+1 ) 2 r+1 + dG H (Ωx ) 8 n=0 n=0

N −1 X

Finally using these estimates in (5.21) we get for T2 N −1 Z tn+1 N −1 N −1 X X

2

τ X 2 n+1 2

utt dt + τ

η n+1 2 |T2 | ≤ τ kη k + 0 0 dG 4 n n=0 t n=0 n=0 N −1  X

n+1 2

2

2

u(t + Cτ h2r+2 ) 1 r+1 + zmin (tn+1 ) r+1 + ut (tn+1 ) 2 H (H

)

H

(Ωx )

n=0

2 + u(tn+1 ) C(H r+1 ) + Cτ k 2q+2

N −1  X





n+1 2

2

u(t ) H q+1 (L2 ) + ut (tn+1 ) H q+1 (L2 ) .

n=0

L (H r+1 )

March 2, 2011

16:10

Operator splitting method for population balance equation

21

Inserting T1 and T2 in (5.16), absorbing the triple norm and the dG norm contributions in the left-hand side and using (3.12), we get N −1 Z N −1 X X



n+1 2 n+1 2 1

η N 2 − 1 η 0 2 + τ d` + τ

η

η˜ 0 0 dG 2 2 2 n=0 n=0 Ω` N −1 Z tn+1 N −1 N −1 X X X

2

n 2

n+1 2

utt dt + τ

η + 2τ

f ≤ τ2 γ n 0 0 0 tn

n=0

+ C h2r τ

n=0

n=0

N −1  X



2

2 (ε + h) u(tn+1 ) H 1 (H r+1 ) + h2 ut (tn+1 ) L2 (H r+1 )

n=0

2

2

+ h zmin (tn+1 ) H r+1 (Ωx ) + u(tn+1 ) C(H r+1 ) 2

+ C k 2q+2 τ

N −1  X



2

2 (ε + h + 1) u(tn+1 ) H q+1 (H 1 ) + ut (tn+1 ) H q+1 (L2 )



n=0

where γ0 = 2, γN = 1 and γn = 3, n = 1, . . . , N − 1. We conclude by applying the Gronwall’s Lemma in the same fashion as in Lemma 3.1. 6. Implementation of numerical method This section indicates the implementation of the operator splitting method in the context of finite element methods. For more details, we refer to Ref. 14. Using the bases Vh = span{φi }, 1 ≤ i ≤ Nx ,

Skq = span{ψk }, 1 ≤ k ≤ N` ,

r,q the tensor product space Sh,k is defined as follows ) ( N` Nx X X r,q Sh,k = v = αik φi (x)ψk (`), αik ∈ R, 1 ≤ i ≤ Nx , 1 ≤ k ≤ N` . i=1 k=1

The finite element functions are represented as unh,k =

N` Nx X X

n ξik φi (x)ψk (`),

X=

i=1 k=1

N` Nx X X

xjl φj (x)ψl (`).

j=1 l=1

We define the matrices Mx , Tx , Dx , Sx ∈ RNx ×Nx by

(Mx )ij = φi (x), φj (x) x , (Dx )ij = ε ∇x φi (x), ∇x φj (x) x

(Tx )ij = b · ∇x φi (x), φj (x) x , (Sx )ij = Sh φi (x), φj (x) . Similarly we define the matrices M` , T` ∈ RN` ×N` as   (M` )kl = ψk (`), ψl (`) , `

(T` )kl =

N`  X i=1

N ` −1  X + + ∂` (Gψk (`)), ψl (`) + [Gψk (`)]i ψl (`+ i ) + Gψk (`0 )ψl (`0 ). Ii

i=1

March 2, 2011

22

16:10

N. Ahmed, G. Matthies, L. Tobiska

For the ease of presentation let us consider (2.1) with source term f = 0. Then the computing scheme for the operator splitting method described in (4.7) and (4.8) is as follows: Within each time interval (tn , tn+1 ], we begin with the x-direction step where we are looking for the solution of the time-dependent convection-diffusion equa˜ n+1 tion (4.7). We compute η ∈ RNx , j = 1, . . . , N` , by solving the linear systems j (Mx + τ Dx + τ Tx + τ Sx )˜ η n+1 = Mx η nj , j

j = 1, . . . , N` .

˜ n+1 η , j

With obtaining the solutions j = 1, . . . , N` , the x-direction step is completed. Then, we continue with the `-direction step where we update the solution from the first step and compute the solution of the one-dimensional transport problem (4.8) by a discontinuous Galerkin method. In this step we solve the linear systems ˜ n+1 (M` + τ T` )η n+1 = M` η , j j

j = 1, . . . , Nx ,

and the obtained solutions η n+1 , j = 1, . . . , N` , are used as input for the time j interval (tn+1 , tn+1 ]. 7. Numerical tests We report in this section the numerical computations illustrating the theoretical results obtained in the previous section. The two-step method (4.7) and (4.8) in the context of finite element method in space and discontinuous Galerkin method in internal coordinate is implemented in the finite element package MooNMD.19 The tests are made in two plus one dimensions, i.e, we consider Ωx = (0, 1)×(0, 1) as two-dimensional domain in space and Ω` = (0, 1) as one-dimensional domain in the internal coordinate. We consider the velocity field b as b1 = b2 = 0.1, the growth rate G(`) = 1 and two different choices for the diffusion coefficient ε, ε = 1 and ε  1. The source term f and the boundary and initial conditions are chosen such that the analytical solution of the problem (2.1) is z(t, `, x, y) = e−0.1t sin(π`) cos(πx) cos(πy). Let en := z(tn ) − unh,k , where z is the exact solution of (2.1) and the numerical solution unh,k is obtained by two-step method (4.7) and (4.8). We use the following notations !1/2 N N X X n 2 n 2 kek0 = τ ke kL2 (L2 ) + τ ke kdG ,

kek1 =

τ

n=1

n=1

N X

N X

ken k2L2 (H 1 ) + τ

n=1

kekDG =

τ

ken k2dG

,

n=1

N Z X n=1

!1/2

Ω`

N X n 2 e d` + τ ken k2dG n=1

!1/2 .

March 2, 2011

16:10

Operator splitting method for population balance equation

23

In order to illustrate the convergence order in time, internal coordinate and space, we use the well known strategy, i.e., the convergence order in time can be obtained by assuming that the mesh sizes k and h are small enough compared to the time-step size τ . In the numerical computations, we have used triangular and quadrilateral meshes which are generated by successive refinement starting from coarsest meshes (level 0) as in Fig. 1 for two-dimensional domain Ωx and a line divided into two cells for one-dimensional domain Ω` .

Fig. 1: Meshes for Ωx on level 0.

Case ε = 1: In this case, the Galerkin finite element method in space is combined with a discontinuous Galerkin method in internal coordinate. For time discretization, the backward Euler time stepping scheme is used with final time T = 1. One can expect a convergence for k · k0 -norm of order O(hr+1 ) and for k · k1 -norms of order O(hr ) using Qr and Pr finite elements in space with sufficiently small time step length τ and mesh size k. The results are presented in Tables 1–4. Tables 1 and 2 show the second order convergence in the k · k0 -norm and first order convergence in the k·k1 -norm for both Q1 and P1 finite elements in space with dG(1) in internal coordinate. The length of the time step was set to be τ = 10−3 and mesh size to k = 1/64. For Q2 and P2 finite elements in space with dG(2) in internal coordinate, the time step length was set to τ = 10−4 and mesh size k = 1/64. The results of Tables 3 and 4 show third order convergence for the k·k0 -norm and second order for the k · k1 -norm. In Tables 5 and 6, the errors and convergence orders for internal coordinate and time are listed. We expect a convergence of order O(k q+1/2 ) in the internal coordinate and a convergence of O(τ ) in time. The errors for dG(1) in internal coordinate with Q1 on level 7 and time step length τ = 2.5 · 10−4 are presented in Table 5. We see that the expected orders of convergence are achieved. The numerical errors and convergence orders in time are listed in Table 6 for dG(1) with k = 1/32 and Q1 on level 6. The theoretically predicted convergence order is achieved. Case ε = 10−9 : In the case of convection-dominated convection-diffusion, we

March 2, 2011

24

16:10

N. Ahmed, G. Matthies, L. Tobiska

Table 1: Errors and rate of convergence in space for Q1 and dG(1), k = 1/64 and τ = 10−3 . Level 0 1 2 3

kek0 error 1.719554e-01 4.746460e-02 1.206219e-02 3.167958e-03

order —— 1.8571 1.9764 1.9289

kek1 error 1.006185 4.892384e-01 2.412003e-01 1.201483e-01

order —— 1.0403 1.0203 1.0054

Table 2: Errors and rate of convergence in space for P1 and dG(1), k = 1/64 and τ = 10−3 . Level 0 1 2 3

kek0 error 2.353104e-01 7.412177e-02 1.981996e-02 5.144843e-03

order —— 1.6666 1.9029 1.9458

kek1 error 1.432599 7.996426e-01 4.113880e-01 2.072235e-01

order —— 0.8413 0.9589 0.9893

Table 3: Errors and rate of convergence in space for Q2 and dG(2), k = 1/64 and τ = 10−4 . Level 0 1 2

kek0 error 1.916287e-02 2.599528e-03 3.354662e-04

order —— 2.8820 2.9540

kek1 error 2.396151e-01 6.137457e-02 1.561139e-02

order —— 1.9650 1.9750

Table 4: Errors and rate of convergence in space for P2 and dG(2), k = 1/64 and τ = 10−3 . Level 0 1 2

kek0 error 3.511498e-02 4.796648e-03 6.138514e-04

order —— 2.8720 2.9661

kek1 error 5.583590e-01 1.526520e-01 3.929766e-02

order —— 1.8710 1.9577

consider local projection as stabilization in space. Discontinuous Galerkin methods of first and second order are used for the discretization in internal coordinate. For time discretization, the backward Euler time stepping scheme is used. The numerical tests are performed using for (Vh , Dh ) the pairs (P1bubble , P0disc ),

March 2, 2011

16:10

Operator splitting method for population balance equation

25

Table 5: Errors and rate of convergence in internal coordinate for dG(1), Q1 on level 6 and τ = 2.5 · 10−4 . k 1/2 1/4 1/8

kek0 6.696513e-02 1.829413e-02 6.521805e-03

—— 1.7398 1.4880

Table 6: Errors and rate of convergence in time for Q1 and dG(1) on level = 6 and k = 1/32. τ 1/10 1/20 1/40 1/80

kek0 error 1.815303e-01 9.577853e-02 4.983170e-02 2.567753e-02

order —— 0.9224 0.9427 0.9566

kek1 error 4.027364 2.170105 1.141479 5.869174e-01

order —— 0.8921 0.9269 0.9597

(P2bubble , P1disc ), (Qbubble , P0disc ), and (Qbubble , P1disc ). The stabilization parameters 1 2 µK have been chosen as µK := µ0 hK

∀K ∈ Th

where µ0 denotes a constant which will be given for each of the test calculations. In Tables 7 and 8 we show the convergence results for space in norm k · kDG . Table 7 shows the error in space with stabilizing parameter µ0 = 5, time step length τ = 10−3 and mesh size k = 1/64 for (Qbubble , P0 ) and (P1bubble , P0 ) with 1 dG(1) in internal coordinate. In Table 8, the convergence results for (Qbubble , P1disc ) 2 bubble disc and (P2 , P1 ) with dG(2) in internal coordinate with µ0 = 5, k = 1/64 and τ = 10−4 are listed. We see that the expected orders of convergence O(hr+1/2 ) are achieved for quadrangles. For smaller mesh size h, the convergence order starts to decrease for triangles. This is because the influence of the error in internal coordinate increases, i.e., the mesh size k is not small enough that one can see the corresponding convergence rate in space for higher order elements. The numerical errors and convergence orders in internal coordinate are listed in , P0 ) with µ0 = 5 on level 7 and τ = 2.5 · 10−4 . The Table 9 for dG(1) and (Qbubble 1 convergence order starts to decrease for small mesh size k since the errors in space have increasing influence. Finally, Table 10 shows the errors and convergence orders in time for (Qbubble , P0 ) on level 6 with µ0 = 2.5 and dG(1) with k = 1/32. We see 1 that the time stepping scheme is of first order convergent.

March 2, 2011

26

16:10

N. Ahmed, G. Matthies, L. Tobiska

Table 7: Errors and rate of convergence in space for (Qbubble , P0 ) and (P1bubble , P0 ) 1 and dG(1), k = 1/64, τ = 10−3 and µK = 5hK . Level 0 1 2 3

(Qbubble , P0 ) 1 ken kDG 1.756772 6.394630e-01 2.280495e-01 8.245890e-02

—— 1.4580 1.4875 1.4678

(P1bubble , P0 ) ken kDG 1.93314 7.247844e-01 2.661525e-01 1.086554e-01

—— 1.4153 1.4453 1.2925

Table 8: Errors and rate of convergence in space for (Qbubble , P1disc ) and 2 bubble disc −4 (P2 , P1 ) and dG(2), k = 1/64, τ = 10 and µK = 5hK . Level 0 1 2 3

(Qbubble , P1disc ) 2 ken kDG 1.272972 2.558153e-01 4.700162e-02 8.010563e-03

—— 2.3151 2.4443 2.5527

(P2bubble , P1disc ) ken kDG 1.234504 2.352103e-01 5.094834e-02 1.222369e-02

—— 2.3919 2.2069 2.0593

Table 9: Errors and rate of convergence in internal coordinate for dG(1) and (Qbubble , P0 ) on level 7 with µK = 5hK and τ = 2.5 · 10−4 . 1 k 1/2 1/4 1/8 1/16

ken kDG 2.493607e-01 9.283060e-02 3.425394e-02 1.446166e-02

—— 1.4256 1.4383 1.2441

Table 10: Errors and rate of convergence in time for dG(1) and (Qbubble , P0 ) on 1 level with µK = 2.5hK and k = 1/32. τ 1/10 1/20 1/40 1/80

ken kDG 8.017623e-01 4.318566e-01 2.270064e-01 1.166372e-01

—— 0.8926 0.9278 0.9607

8. Conclusion In this paper we have been concerned with the numerical solution of the population balance equation with one internal coordinate posed on the domain Ω` × Ωx in d + 1 dimension. We proposed an operator splitting method to reduce the original

March 2, 2011

16:10

Operator splitting method for population balance equation

27

problem into two subproblems. The method combines the continuous finite element method (and local projection stabilization) in space with discontinuous Galerkin method in internal coordinate. We have considered first order backward Euler time stepping scheme. Under certain regularity of exact solution, we have derived error estimates for the two-step method, i.e., using polynomials of degree r in space and of degree q in internal coordinate the error is of order O(τ + hr+1/2 + k q+1/2 ) when ε  1 and O(τ + hr + k q+1/2 ) when ε = 1. The application of discontinuous Galerkin method makes the mass matrix corresponding to the internal coordinate diagonal, which leads to the feasibility of the implementation without any projection between the two-steps in the computation process. Computational results shown in Section 7 confirms the theoretical prediction of error estimates.

References 1. N. Ahmed, G. Matthies, L. Tobiska, and H. Xie, Discontinuous Galerkin time stepping with local projection stabilization for transient convection-diffusion-reaction problems, Preprint 10-24, Fakult¨ at f¨ ur Mathematik, Otto-von-Guericke-Universit¨ at Magdeburg, 2010. 2. R. Becker and M. Braack, A finite element pressure gradient stabilization for the Stokes equations based on local projections, Calcolo, 38 (2001), pp. 173–199. 3. R. Becker and M. Braack, A finite element pressure gradient stabilization for the Stokes equations based on local projections, Calcolo, 38 (2001), pp. 173–199. 4. M. Braack and E. Burman, Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method, SIAM J. Numer. Anal., 43 (2006), pp. 2544–2566 (electronic). 5. E. Burman, A unified analysis for conforming and nonconforming stabilized finite element methods using interior penalty, SIAM J. Numer. Anal., 43 (2005), pp. 2012– 2033 (electronic). ´ ndez, Finite element methods with symmetric sta6. E. Burman and M. A. Ferna bilization for the transient convection-diffusion-reaction equation, Comput. Methods Appl. Mech. Engrg., 198 (2009), pp. 2508–2519. 7. E. Burman and P. Hansbo, The edge stabilization method for finite elements in CFD, in Numerical mathematics and advanced applications, Springer, Berlin, 2004, pp. 196–203. 8. P. G. Ciarlet, The finite element method for elliptic problems, North-Holland Publishing Co., Amsterdam, 1978. Studies in Mathematics and its Applications, Vol. 4. 9. B. Cockburn, An introduction to the discontinuous Galerkin method for convectiondominated problems, in Advanced numerical approximation of nonlinear hyperbolic equations (Cetraro, 1997), vol. 1697 of Lecture Notes in Math., Springer, Berlin, 1998, pp. 151–268. 10. R. Codina, Stabilization of incompressibility and convection through orthogonal subscales in finite element methods, Comput. Methods Appl. Mech. Engrg., 190 (2000), pp. 1579–1599. 11. R. Codina and J. Blasco, Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales, Comput. Vis. Sci., 4 (2002), pp. 167–174. ´ jek, and K. Svadlenka, Space-time discontinuous Galerkin 12. M. Feistauer, J. Ha

March 2, 2011

28

13.

14.

15. 16. 17. 18. 19. 20.

21. 22.

23.

24. 25. 26.

27.

28. 29.

30.

31.

16:10

N. Ahmed, G. Matthies, L. Tobiska

method for solving nonstationary convection-diffusion-reaction problems, Appl. Math., 52 (2007), pp. 197–233. S. Ganesan, Population balance equations, Streamline-Upwind Petrov-Galerkin finite element methods, operator-splitting method, backward Euler scheme, error analysis, Preprint 1531, WIAS, Berlin, 2010. S. Ganesan and L. Tobiska, Implementation of an operator splitting finite element method for high-dimensional parabolic problems, Preprint 11-04, Fakult¨ at f¨ ur Mathematik, Otto-von-Guericke-Universit¨ at Magdeburg, 2010. J.-L. Guermond, Stabilization of Galerkin approximations of transport equations by subgrid modeling, M2AN Math. Model. Numer. Anal., 33 (1999), pp. 1293–1316. E. Hansen and A. Ostermann, Dimension splitting for evolution equations, Numer. Math., 108 (2008), pp. 557–570. , Dimension splitting for time dependent operators, Discrete Contin. Dyn. Syst., (2009), pp. 322–332. L. He and L. Tobiska, The two-level local projection stabilization as an enriched one-level approach, Adv. Comput. Math. (to appear), (2011). V. John and G. Matthies, MooNMD—a program package based on mapped finite element methods, Comput. Vis. Sci., 6 (2004), pp. 163–169. V. John, M. Roland, T. Mitkova, K. Sundmacher, L. Tobiska, and A. Voigt, Simulations of population balance systems with one internal coordinate using finite element methods, Chem. Eng. Sci., 64 (2009), pp. 733–741. ¨ ranta, An analysis of the discontinuous Galerkin method C. Johnson and J. Pitka for a scalar hyperbolic equation, Math. Comp., 46 (1986), pp. 1–26. ˇur, B. Malengier, and M. Remeˇ ´ , Convergence of an operator splitJ. Kac s´ıkova ting method on a bounded domain for a convection-diffusion-reaction system, J. Math. Anal. Appl., 348 (2008), pp. 894–914. K. H. Karlsen, K. Brusdal, H. K. Dahle, S. Evje, and K.-A. Lie, The corrected operator splitting approach applied to a nonlinear advection-diffusion problem, Comput. Methods Appl. Mech. Engrg., 167 (1998), pp. 239–260. K. H. Karlsen and K.-A. Lie, An unconditionally stable splitting scheme for a class of nonlinear parabolic equations, IMA J. Numer. Anal., 19 (1999), pp. 609–635. K. H. Karlsen and N. H. Risebro, An operator splitting method for nonlinear convection-diffusion equations, Numer. Math., 77 (1997), pp. 365–382. ¨ li, A heterogeneous alternating-direction method for a D. J. Knezevic and E. Su micro-macro dilute polymeric fluid model, M2AN Math. Model. Numer. Anal., 43 (2009), pp. 1117–1156. P. Knobloch, On the application of local projection methods to convection-diffusionreaction problems, in BAIL 2008—boundary and interior layers, vol. 69 of Lect. Notes Comput. Sci. Eng., Springer, Berlin, 2009, pp. 183–194. , A generalization of the local projection stabilization for convection-diffusionreaction equations, SIAM J. Numer. Anal., 48 (2010), pp. 659–680. J. Koch, Effiziente Behandlung von Integraloperatoren bei populationsdynamischen Modellen, PhD thesis, Otto-von-Guericke-Universit¨ at Magdeburg, Fakult¨ at f¨ ur Mathematik, 2005. D. Lanser and J. G. Verwer, Analysis of operator splitting for advection-diffusionreaction problems from air pollution modelling, J. Comput. Appl. Math., 111 (1999), pp. 201–216. Numerical methods for differential equations (Coimbra, 1998). P. Lasaint and P.-A. Raviart, On a finite element method for solving the neutron transport equation, in Mathematical aspects of finite elements in partial differential equations (Proc. Sympos., Math. Res. Center, Univ. Wisconsin, Madison, Wis., 1974),

March 2, 2011

16:10

Operator splitting method for population balance equation

32. 33.

34.

35.

36. 37. 38.

39. 40.

41. 42. 43.

44.

45.

29

Math. Res. Center, Univ. of Wisconsin-Madison, Academic Press, New York, 1974, pp. 89–123. Publication No. 33. W. Layton, A connection between subgrid scale eddy viscosity and mixed methods, Appl. Math. Comput., 133 (2002), pp. 147–157. G. I. Marchuk, Splitting and alternating direction methods, in Handbook of numerical analysis, Vol. I, Handb. Numer. Anal., I, North-Holland, Amsterdam, 1990, pp. 197–462. G. Matthies, P. Skrzypacz, and L. Tobiska, A unified convergence analysis for local projection stabilisations applied to the Oseen problem, M2AN Math. Model. Numer. Anal., 41 (2007), pp. 713–742. , Stabilization of local projection type applied to convection-diffusion problems with mixed boundary conditions, Electron. Trans. Numer. Anal., 32 (2008), pp. 90– 105. A. Mersmann, Batch precipitation of barium carbonate, Chem. Eng. Process., 38 (1993), pp. 6177–6184. , Crystallization and precipitation, Chem. Eng. Process., 38 (1999), pp. 345–353. G. M. Phillips, Interpolation and approximation by polynomials, CMS Books in Mathematics/Ouvrages de Math´ematiques de la SMC, 14, Springer-Verlag, New York, 2003. W. H. Reed and T. R. Hill, Triangular mesh methods for the neutron transport equation, Tech. Report LA-UR-73-479, Los Alamos Scientic Laboratory, (1973). H.-G. Roos, M. Stynes, and L. Tobiska, Robust numerical methods for singularly perturbed differential equations, vol. 24 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, second ed., 2008. Convection-diffusion-reaction and flow problems. G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal., 5 (1968), pp. 506–517. ´e, Galerkin finite element methods for parabolic problems, vol. 25 of Springer V. Thome Series in Computational Mathematics, Springer-Verlag, Berlin, second ed., 2006. L. Tobiska, Recent results on local projection stabilization for convection-diffusion and flow problems, in BAIL 2008—boundary and interior layers, vol. 69 of Lect. Notes Comput. Sci. Eng., Springer, Berlin, 2009, pp. 55–75. L. Tobiska and C. Winkel, The two-level local projection type stabilization as an enriched one-level approach. A one-dimensional study, Int. J. Numer. Anal. Model., 7 (2010), pp. 520–534. N. N. Yaneko, The method of fractional steps, Springer, Berlin, 1971.