SE Method for ... - ACS Publications

Article pubs.acs.org/IECR

Application of Space−Time CE/SE Method for Solving Gas−Solid Reaction and Chemotaxis Models Shamsul Qamar*,†,‡ and Waqas Ashraf‡ †

Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstr. 1, 39106 Magdeburg, Germany Department of Mathematics, COMSATS Institute of Information Technology, Park Road, Chak Shahzad Islamabad, Pakistan

‡

ABSTRACT: This work is concerned with the simulation of two different convection−diffusion−reaction type partial differential equations. The first one is the one-dimensional heterogeneous gas−solid reaction model encountered in a variety of chemical engineering problems. The exothermic nature of such reactions enlarges the gradients of temperature and concentrations that generate steep reaction fronts. The second one is the two-dimensional Keller−Segel model of the chemotaxis that generates delta-type singularities in the solution during a finite time. In its simplest form, the model is a nonlinear-coupled system of the convection−diffusion equation for the cell density and the reaction−diffusion equation for the chemoattractant concentration. The space−time conservation element and solution element (CE/SE) method is proposed for approximating both systems. The method has capabilities to capture sharp discontinuities of the solutions without spurious oscillations. To validate the accuracy and efficiency of the method, several case studies are carried out. The results of the current method are compared with the staggered central schemes.

■

INTRODUCTION The convection−diffusion−reaction (CDR) type partial differential equations (PDEs) have various engineering applications, namely, the modeling of fluid dynamics, chemical tubular reactors, and aerodynamics. Two CDR systems are numerically investigated in this manuscript. The first CDR system simulates the heterogeneous gas−solid reaction (GSR) that can be used in the modeling of heavy oil recovery, pyrolysis of coal and biomass, roasting and reduction of ores, incineration of waste, solids combustion, acid gases absorption by lime, deposition of reactive vapor phase, ceramic materials production, gasification of coal, manufacturing of catalyst, and so on.1,2 The exothermic nature of GSR produces large gradients of temperature that generates steep reaction fronts in the concentration profile. The capturing of these fronts is very difficult for less accurate numerical methods. Therefore, efficient and accurate numerical methods are essential for producing physically realistic solutions. The second CDR system is a classical Keller−Segel (KS) model of the chemotaxis that admits solutions that generate delta-type singularities within a finite simulation time.3−6 Chemotaxis is a phenomenon in which cells change their state of movement in the presence of certain chemical substances and thus approach chemically favorable environments and avoid unfavorable ones. The new chemotaxis model introduced by Keller and Segel7 is considered as a regularized KS system. This regularization is based on a basic physical principle: boundedness of the chemotactic convective flux, which should depend on the gradient of the chemoattractant concentration in a nonlinear way. Solutions of the new system may develop spiky structures that model the concentration phenomenon.3,4 The simulation of the convection-dominated CDR system is a challenging task for the numerical schemes due to steep reaction gradients and narrow peaks in the solutions. In such situations, a better solution can be obtained by conducting a numerical simulation on a refined grid. However, such computations are © 2012 American Chemical Society

highly expensive at large scale. On the other hand, the numerical results on coarse grids have dissipated/smeared wave fronts that can make the interpretation of experimental data less effective and produces negative affects on the design of industrial process. For that reason, efficient and accurate numerical methods are needed to get desired results. During the past decade, researchers in this field have given their attention to introduce better numerical techniques for solving such CDR PDEs. Flux-limiting finite volume schemes were applied to approximate heterogeneous gas−solid reaction models.2 The semidiscrete centralupwind and the discontinuous Galerkin methods were implemented to solve the two-dimensional KS-system.3,4 Other articles on these models are also available in the literature.2,8−11 In this manuscript, the space−time CE/SE method is applied to solve these CDR systems. In the case of one-dimensional heterogeneous GSR model, the original CE/SE method is applied,12 while a variant CE/SE method on rectangular mesh elements is employed for approximating two-dimensional KSmodel.13 The variant CE/SE method is suitable because of its simplified formulation and implementation on regular rectangular meshes as compared to the original two-dimensional CE/SE method on unstructured triangular meshes.14 The main features of the CE/SE method include the unified space−time treatment, the use of conservation elements (CEs) and solution elements (SEs), as well as its ability to capture sharp discontinuities without employing Riemann solvers. The technique has already been implemented to complex problems in different areas, namely, problems related to unsteady flows,12,14,15 vortex dynamics in aeroacoustics,16 diffusion problems,17 viscous flows,10 inviscid and axisymmetric flows,18 supersonic jets,19 Received: Revised: Accepted: Published: 9173

March 1, 2012 June 12, 2012 June 12, 2012 June 12, 2012 dx.doi.org/10.1021/ie3005622 | Ind. Eng. Chem. Res. 2012, 51, 9173−9185

Industrial & Engineering Chemistry Research

Article

shallow water magnetohydrodynamics,20 and magnetohydrodynamics.21 The rest of the paper is arranged as follows. In Section 2, the one-dimensional heterogeneous GSR model is given, followed by the derivation of one-dimensional CE/SE method.12 In Section 3, the two-dimensional Keller−Segel (KS) model is given. Afterward, the two-dimensional CE/SE method on uniform regular rectangular grids is briefly described.13 Numerical test cases are presented in Section 4. This section contains several one- and two-dimensional numerical test problems that are validated against the staggered central schemes.22,23 Finally, conclusions and remarks are given in Section 5.

coordinate. Further, the subscripts g stands for the gas and s for the solid. The following substitutions are introduced to get the dimensionless form of the model 1−32,24 cg c E(T − T *) , cdg = , cds = s , xd = 2 c c RT * g,inlet s,inlet

ψ=

td =

Ek 0(−ΔH )cg

ε

= Dε

∂t

∂ 2cg ∂x 2

−u

∂cg ∂x

− k 0cgcse−E / RT

∂c (1 − ε) s = −k 0cgcse−E / RT ∂t

uρκp λt * , PeH = ρκp λ

γg =

ρκpRT *

εE(−ΔH )cg,inlet

(6)

∂cg ∂x

∂x

= u(cg

inlet

= 0 at x = ∞

− cg) at x = 0

, γs =

ρκpRT *

(1 − ε)E(−ΔH )cg,inlet

,β =

RT * E

These are standard group of parameters available in the literature of combustion to get dimensionless form. With the help of this group of scaling, the dimensionless equations are given as ⎛ ψ ⎞ ∂ψ ∂ 2ψ ∂ψ = − PeH + cdgcdsexp⎜ ⎟ 2 ∂td ∂xd ⎝ βψ + 1 ⎠ ∂xd ∂cdg

=

(15)

2 ⎛ ψ ⎞ Pe ∂cdg 1 ∂ cdg − γgcdgcdsexp⎜ − M ⎟ 2 Le ∂xd Le ∂xd ⎝ βψ + 1 ⎠

⎛ ψ ⎞ ∂cds = −γscdgcdsexp⎜ ⎟ ∂td ⎝ βψ + 1 ⎠

(3)

cs = cs0 , t = 0, x ∈ [0,∞]

− εD

2

(16)

(5)

∂cg

λt * λ , Le = ρκp ρκpD

(14)

(2)

cg = cg , t = 0, x ∈ [0,∞]

∂T = 0 at x = ∞ ∂x

λt * u , PeM = ρκp εD

(13)

∂td

(4)

∂T = ρκpu(Tinlet − T ) at x = 0 ∂x

ρκpRT * ⎛ −E ⎞ ⎛ E ⎞ ⎟t , t * = ⎟ exp⎜ exp⎜ ⎝ RT * ⎠ Ek 0(−ΔH )cg,inletcs,inlet ⎝ RT * ⎠

2

(1)

T = T0 , t = 0, x ∈ [0, ∞]

−λ

i

ρκpRT *

x* =

(17)

The dimensionless initial and boundary conditions are given as

The initial and boundary condition are given as2,24

0

2

2

2

∂cg

x

(12)

GOVERNING EQUATION OF GSR MODEL This model is based on the pseudohomogeneous one-phase model in one-space dimension.24 The reaction takes place between the solid phase and a gaseous oxidizer (nitrogen or oxygen), the product of which is again a gas or solid. In the solid phase, Fourier’s law of heat conduction holds. The sintering effects are neglected, and it is assumed that porosity of the system remains constant throughout the process. The effects of radiation are included into thermal conductivity, and a first-order reaction is assumed corresponding to solid and gas reactants. All physical properties of the system are considered constant. The governing differential equation describing propagation of reaction fronts form a set of nonlinear mixed parabolic-hyperbolic PDEs2,24 ∂T ∂T ∂T = λ 2 − uρκp + k 0( −ΔH )cgcse−E / RT ∂t ∂x ∂x

λt *

(11)

■

ρκp

ρκp

ψ = ψ0 , x ∈ [0,∞], t = 0

(18)

cdg = cdg , x ∈ [0,∞], t = 0

(19)

cds = cds0 , x ∈ [0,∞], t = 0

(20)

0

∂cdg

= −PeM (1 − cdg), t ≥ 0, x = 0

(7)

∂xd

(8)

∂ψ = −PeH (ψ − ψ0), t ≥ 0, x = 0 ∂xd ∂cdg

(9)

∂xd

= 0, x = ∞

(21)

(22)

(23)

This completes the one-dimensional formulation of gas−solid reactions model. The One-Dimensional CE/SE Method. The CE/SE method was developed to solve conservation laws governing the motion of fluids.12 The method was designed to overcome some of the shortcomings of well established methods. Let us first rewrite eqs 15−17 as

(10)

In the above model, ρκp = ρκpε + ρsκps(1 − ε) is the effective heat capacity, T denotes the temperature, c represents the reactant concentration, ρ denotes the density, κp stands for the heat-capacity, λ is the averaged thermal conductivity, k0 is the pre-exponential factor, ΔH is the reaction heat, E represents the activation energy, R denotes the universal gas constant, D is the coefficient of molecular diffusion, u is the velocity, ε is the porosity of solid material, t is the time, and x represents the space

∂w ∂w ∂ 2w +a − b 2 = Q(w) ∂t ∂x ∂x 9174

(24)

dx.doi.org/10.1021/ie3005622 | Ind. Eng. Chem. Res. 2012, 51, 9173−9185


Article

Figure 1. (a) Space−time staggered grid near SE(j,n). (b) CE−(j,n) and CE+(j,n), respectively. (c) CE(j,n).

the inner shaded space−time region in Figure 1. It includes segments of horizontal and vertical lines and their close neighborhood. The actual size of the neighborhood is not a concern. At point (t,x) in SE, the approximation of wk(t,x), f k(t,x) and Fk(t,x) are given as w*k(t,x;j,n),f *k(t,x;j,n) and F*2(t,x;j,n), respectively. Let

where ⎛ PeH 0 0 ⎞ ⎛1 ⎛ψ ⎞ ⎜ ⎟ ⎜ PeM ⎟ ⎜c ⎟ ⎜ w = ⎜ dg ⎟, a = ⎜ 0 0 ⎟, b = ⎜ 0 ⎜ Le ⎜ ⎟ ⎜⎜ ⎟⎟ ⎜ ⎝ cds ⎠ ⎝0 ⎝ 0 0 0⎠

⎛ ⎛ ψ ⎞ ⎞ ⎜ cdgcdsexp⎜ ⎟ ⎟ ⎝ βψ + 1 ⎠ ⎟ ⎜ ⎟ ⎜ ⎛ ψ ⎞⎟ ⎜ Q(w) = ⎜−γgcdgcdsexp⎜ ⎟ ⎝ βψ + 1 ⎠ ⎟⎟ ⎜ ⎜ ⎛ ⎞⎟ ⎜ −γ c c exp⎜ ψ ⎟ ⎟ ⎟ ⎜ s dg ds ⎝ βψ + 1 ⎠ ⎠ ⎝

0 0⎞ ⎟ 1 0⎟ Le ⎟ ⎟ 0 0⎠

(25)

w*k (t , x ; j , n) = (wk)nj + (wkt )nj (t − t n) + (wkx)nj (x − xj) (30)

where (wk)nj , (wkx)nj and (wkt)nj are constant in SE. Similarly, we can define f *k (t ,x ; j ,n) = (fk )nj + (fkt )nj (t − t n) + (fkx )nj (x − xj) (31)

By employing chain rule, we obtain 3

(26)

(fkx )nj =

Now, eq 24 can be rewritten as ∂w ∂ ⎛⎜ ∂w ⎞ + aw − b ⎟ = Q(w) ∂t ∂x ⎝ ∂x ⎠

m=1

∑ (fk ,m )nj (wmt )nj m=1

(33)

where fk ,m =

(28)

where f = aw − bwx. Let x1 = t and x2 = x be the coordinates of Euclidean space (E2) and let for k = 1,2,3, Fk = [wk,f k,−Qk]T. After applying the divergence theorem, eq 28 becomes

∮s(U) Fk × ds = 0, k = 1, 2, 3

(32)

3

(fkt )nj =

(27)

Therefore, the above equation has the following conservation form ∂w ∂f + =Q ∂t ∂x

∑ (fk ,m )nj (wmx)nj

∂fk ∂wm

, k ,m = 1,2,3

(34)

The Jacobian matrix formed by f k,m, k,m = 1,2,3, with k and m being the row and column indices, respectively, is given as

⎛ PeH 0 0 ⎞ ⎜ ⎟ PeM ⎟ ⎜ A=⎜ 0 0⎟ Le ⎜⎜ ⎟⎟ ⎝ 0 0 0⎠

(29)

where k stands for the total number of equations, and s(U) denotes the boundary of a space−time domain U. Equation 29 is restricted to the conservation element (CE) in the space−time domain that permits the flow variables discontinuities. The actual numerical integration is done directly on the solution elements (SEs). Let Ω denote the set of grid points (j,n) with n = 0, ± (1/2), ± 1, ± (3/2), ··· and for every n, j = n ± (1/2),n ± 1, n ± (3/2), ···. With each (j,n) ∈ Ω one SE is associated. Assume that SE(j,n) is

(35)

Because Fk = (f k,wk, −Qk), one can write F *k (t ,x:j ,n) = (f *k (t ,x:j ,n),w*k (t ,x:j ,n),−Q *k (t ,x:j ,n)) (36)

Because of are functions 9175

their definitions, (f k)nj are functions of (wk)nj , (f kx)nj of (wk)nj and (wkx)nj , and (f kt)nj of (wk)nj and (wkt)nj . dx.doi.org/10.1021/ie3005622 | Ind. Eng. Chem. Res. 2012, 51, 9173−9185


Article

Moreover, it is assumed that for all (t,x) in SE, w*k, f *k and Q*k satisfy the eq 28, that is, ∂w*k (t ,x ;j ,n) ∂f *k (t ,x ;j ,n) + = Q k(w*k ) ∂t ∂x

Now integrating over the mesh point (j + (1/2),n − (1/2)) we get

∫x

(37)

xj + 1/2

j

According to eq 28, eq 37 can be equivalently written as (wkt )nj = −(fkx )nj + (Q k)nj

− (38)

=−

(39)

Therefore, (wk)nj and (wkx)nj are the independent variables of the scheme at every grid point (j,n). In Figure 1, the rectangular nonoverlapped regions represent the CEs. The two CEs, CE−(j,n) and CE+(j,n), correspond to every inner grid point (j,n). Thus CE at point (j,n) is the combination of CE−(j,n) and CE+(j,n). Moreover, the CE−(j,n) boundary contains the subset of SE(j,n) and SE(j−(1/2), n−(1/2)), and that of CE+(j,n) is obtained from the subset of SE(j,n) and SE(j+(1/2),n−(1/2)). Consider the integral representation of eq 28:

∫∂Ω wkdx − fk dt + ∫Ω Q kdxdt = 0, k = 1,2,3

−

∫x

xj

∫t

=−

j − 1/2

=−

∫t

tn n − 1/2

∫x

xj

n − 1/2

f *k (t ,xj − 1/2)dt

Q k dx dt

j

∫t

(44)

(1 − νk) (1 + νk) 1/2 1/2 (wk)nj −−1/2 (wk)nj +−1/2 + 2 2

+

Δt [(Q k)nj − 1/2 + (Q k)nj +−11/2 ] 4

(46)

and the subtraction of eqs 43 and 45 implies + n (wkx )j =

f *k (t ,xj − 1/2)dt

Q k dx dt

(45)

+ n − 1/2 + n − 1/2 ) j − 1/2 − (wkx ) j + 1/2 ] + (1 − νk2 − ξk)[(wkx

(41)

tn n − 1/2

ΔxΔt (Q k)nj +−11/2 4

(wk)nj =

f *k (t ,xj)dt

j − 1/2

− (1 − νk2)

1/2 1/2 − (wk)nj +−1/2 [(wk)nj −−1/2 ] 2(1 − νk2 + ξk) −1 + n − 1/2 + (1 − νm2 − ξk)[(1 − νk)(wkx ) j − 1/2 2(1 − νk2 + ξk) Δt + n − 1/2 + (1 + νk)(wkx ) j + 1/2 ] − 8(1 − νk2 + ξk)

×[(1 − νk)(Q k)nj − 1/2 − (1 + νk)(Q k)nj +−11/2 ]

(42)

(47)

Here

Substituting the values of w* and f * and using f k = akwk − bkwkx and f kt = akwkt = ak(−akwkx) = −a2k kkx, we get after simplification

1 − νk2 + ξk ≠ 0

⎛ 2⎞ Δx ⎞⎟ n ⎛ Δt ⎞ 1/2 ⎛ 1/2 Δx ⎜ + (wkx)nj −−1/2 (wk)nj −−1/2 ⎟ − ak (wk) j ⎜ ⎟ ⎜ ⎝ 2 ⎠ ⎝ 2 ⎠ ⎝ 8 ⎠

(48)

and νk =

⎛ Δt 2 ⎞ ⎛ Δt ⎞ n ⎛ Δx ⎞ + bk (wkx)nj ⎜ ⎟ − ak2(wkx)nj ⎜ ⎟ − (wk) j ⎜ ⎟ ⎝ 2 ⎠ ⎝ 2 ⎠ ⎝ 8 ⎠

ak Δt 4b Δt Δx + n , ξk = k 2 , (wkx )j = (wkx)nj 4 Δx (Δx)

(49)

This completes the derivation of the one-dimensional CE/SE method.

⎛ Δx 2 ⎞ n − 1/2 ⎛ Δt ⎞ + (wkx)nj ⎜ ⎟ + a(wk) j − 1/2 ⎜ ⎟ ⎝ 2 ⎠ ⎝ 8 ⎠

■

THE KELLER−SEGEL (KS) MODEL The dimensionless KS system is given as3,4

⎛ 2⎞ Δt ⎟⎞ 1/2 ⎛ 1/2 Δt ⎜ − bk (wkx)nj −−1/2 − ak2(wkx)nj −−1/2 ⎜ ⎟ ⎝ 2 ⎠ ⎝ 8 ⎠ ΔxΔt =− (Q k)nj − 1/2 4

xj + 1/2

n − 1/2

The summation of eqs 43 and 45 gives

tn

w*k (t n ,x)dx +

∫x

tn

⎛ Δt 2 ⎞ ⎛ Δt ⎞ − bk (wkx)nj ⎜ ⎟ + ak2(wkx)nj ⎜ ⎟ ⎝ 2 ⎠ ⎝ 8 ⎠

where, ∂Ω is the boundary of Ω. On integrating the rectangular region over the mesh point (j ± (1/2),n − (1/2)), we obtain (wk)nj and (wkx)nj . Integration over the mesh point (j − (1/2), n − (1/2)) gives ,x)dx −

n − 1/2

∫t

⎛ Δx 2 ⎞ ⎛ Δx ⎞ n ⎛ Δt ⎞ − (wk)nj ⎜ ⎟ − (wkx)nj ⎜ ⎟ + ak (wk) j ⎜ ⎟ ⎝ 2 ⎠ ⎝ 2 ⎠ ⎝ 8 ⎠

(40)

w*k (t

w*k (t n ,x)dx +

f *k (t ,xj + 1/2)dt

⎛ 2⎞ Δt ⎞⎟ 1/2 ⎛ 1/2 Δt ⎜ + bk (wkx)nj +−1/2 + ak2(wkx)nj +−1/2 ⎜ ⎟ ⎝ 2 ⎠ ⎝ 8 ⎠

w*k (t ,x ;j ,n) = (wk)nj + (wkx)nj ((x − xj) − ak (t − t n))

j − 1/2

∫t

tn

n − 1/2

⎛ 2⎞ Δx ⎞⎟ n − 1/2 ⎛ Δt ⎞ 1/2 ⎛ 1/2 Δx ⎜ − (wkx)nj +−1/2 (wk)nj +−1/2 ⎟ − ak (wk) j + 1/2 ⎜ ⎟ ⎜ ⎝ 2 ⎠ ⎝ 2 ⎠ ⎝ 8 ⎠

Thus, diffusion in eq 24 puts no impact in eq 39. Resultantly, there is no effect of the diffusion term on the time variation of w*k(t,x;j,n) within SE. Thus, one gets from eqs 30 and 39

∫x

xj + 1/2

tn

After some simplification, we obtain

(wkt )nj = −ak (wkx)nj + (Q k)nj

n − 1/2

∫x

∫t

j

Equations 30 and 31 are first-order Taylor’s expansions, therefore eq 27 implies for f kx = akwkx

xj

w*k (t n − 1/2 ,x)dx −

⎧ ⎪ ρt + ∇ × (ξρ∇c ) = Δρ ⎨ ⎪ ⎩ ct = Δc − c + ρ

(43) 9176

(50)



Article

Table 1. Data Used in Numerical Simulation parameters

γs

γg

β

ψ0

cdg0

cds0

PeH

PeM

Le

problem 1 problem 2 problem 3

2 0.214 2.0

0.76 0.1 0.1

0.15 0.0793 0.1

−0.5 −4.667 −1.0

0 0 0.2

1 1 10

100 50 10

200 500 200

1 10 1.2

Figure 3. Results of problem 1 at td = 0.6.

Here, ρ(t,x,y) represent the cell density, c(t,x,y) denotes the concentration of chemoattractant, and ξ is a sensitivity constant of chemotactic. The above model can be replaced by a simpler model if it is assumed that the concentration c changes at smaller time scale compared to the density ρ. Then the equation for c becomes an elliptic equation3,4

Figure 2. Problem 1: Errors at different numbers of grid points. 9177



Article


ρt + (ξρu)x + (ξρv)y = ρxx + ρyy


ut − ρx = uxx + uyy − u

ρt + ∇ × (ξρ∇c) = Δρ Δc − c + ρ = 0

vt − ρy = vxx + vyy − v (51)

This system can be expressed in the form of a CDR system as

Let us consider the two-space dimensions, and let (u,v) = (cx, cy). On differentiating the second equation of eq 50 with respect to x and y and on rewriting the system in equivalent two-space dimensions, we obtain

wt + f(w)x + g(w)y = Δw + Q(w)

(54)

where w = (ρ, u, v)T, f = (ξρu, − ρ,0)T, g = (ξρv, 0, −ρ)T, and Q = (0, −u, −v)T. The Jacobian matrices of f and g are given as

ρt + (ξρu)x + (ξρv)y = ρxx + ρyy

⎛ ξu ξρ 0 ⎞ ⎛ ξv 0 ξρ ⎞ ∂g ⎜ ⎟ ⎟ ⎜ ∂f = ⎜−1 0 0 ⎟ and =⎜0 0 0⎟ ∂w ∂w ⎜ ⎟ ⎟ ⎜ ⎝ 0 0 0⎠ ⎝− 1 0 0 ⎠

(cx)t = (cx)xx + (cx)yy − cx + ρx (cy)t = (cy)xx + (cy)yy − cy + ρy

(53)

(52)

(55)

The Two-Dimensional CE/SE Method. The complete derivations of this method are given in the article by

After substituting u = cx and v = cy, we get 9178



Article

Table 2. Comparison of Numerical Errors (fraction) at Different Grid Points N = 200

N = 103

N = 500

CPU (s)

methods

ψ

cdg

cds

ψ

cdg

cds

ψ

cdg

cds

N = 200

CE/SE Superbee Koren Central Upwind

0.120 0.124 0.129 0.130 0.170

0.030 0.039 0.040 0.041 0.110

0.008 0.009 0.009 0.009 0.012

0.079 0.083 0.087 0.090 0.130

0.014 0.021 0.023 0.025 0.091

0.003 0.003 0.003 0.004 0.005

0.017 0.020 0.024 0.026 0.087

0.009 0.016 0.017 0.018 0.069

0.001 0.001 0.001 0.001 0.002

0.57 0.63 0.56 0.55 0.44

Figure 6. Results of problem 4 at t = 1 × 10−4. 9179



Article

Figure 7. Results of problem 4 at t = 4.4 × 10−6, t = 4.4 × 10−5 and t = 4.4 × 10−4 (from top to bottom).

Zhang et al.13 Here, we give the final formulation of the scheme on regular rectangular grid. The current CE/SE method is suitable because of its simplified formulation and because of its easy implementation on regular rectangular meshes as compared to the original two-dimensional CE/SE method on unstructured triangular meshes.14,15 Let N and M denote the larger integers in the space directions x and y, respectively. Consider a Cartesian grid on a rectangular domain covered by cells Ωij: = [xi−(1/2), xi+(1/2)] × [yj−(1/2), yj+(1/2)] for i = 1, 2, ···, N and j = 1, 2, ···, M. The representative coordinates of Ωij are given by (xi, yj). The scheme has the following formulation for eq 5413

w in, j =

1 n − 1/2 n − 1/2 n − 1/2 (w i − 1, j + w in+−1,1/2 j + w i,j−1 + w i,j+1 4 n − 1/2 n − 1/2 n − 1/2 + (sx)in−−1,1/2 j − (sx )i + 1, j + (sy)i , j − 1 − (sy)i , j + 1 )

λx ((wx)in,−j 1/2 − (wx)in−−1,1/2 j ) 2 λy Δt n − 1/2 Q + ((wy)in,−j 1/2 − (wy)in,−j −1/2 1 ) + 2 2 i,j +

(56)

where λx = Δt/Δx, λy = Δt/Δy and (sx)in,−j 1/2 = 9180

Δx (wx)in,−j 1/2 + λx(f in,−j 1/2 + (f t)in,−j 1/2 ) 2

(57)



Article

Figure 8. Results of problem 5 at t = 5 × 10−6, t = 1 × 10−4 and t = 4.5 × 10−4 (from top to bottom).

(sy)in,−j 1/2 =

Δy (wy)in,−j 1/2 + λy(g in,−j 1/2 + (gt)in,−j 1/2 ) 2

To obtain

(wx)i,jn

and

(wy)i,jn ,

Similarly (58)

a central difference-type

⎛ Δt n⎞ ⎟ u1,y = u4,y = ⎜w in,−j +1/2 (wt )in,−j +1/2 1 + 1 − w i , j /Δy ⎝ ⎠ 2

(61)

⎛ ⎞ Δt ⎟ u 2,y = u3,y = ⎜w in,j − w in,−j −1/2 (wt )in,−j −1/2 1 − 1 /Δy ⎝ ⎠ 2

(62)

reconstruction procedure is employed. Let u1,x = u 2,x

⎛ ⎞ Δt ⎟/Δx (wt )in−−1,1/2 = ⎜w in,j − w in−−1,1/2 j − j ⎝ ⎠ 2

(59)

Then ⎛ Δt n⎞ ⎟ u3,x = u4,x = ⎜w in+−1,1/2 (wt )in+−1,1/2 j + j − w i , j /Δx ⎝ ⎠ 2

ϕl =

(60) 9181

(u i ,x)2 + (u i ,y)2 , for i = 1,...,4

(63)



Article

Figure 9. Results of problem 6 at t = 3 × 10−4.

where 1 ≤ θ ≤ 2. This concludes the derivation of twodimensional CE/SE method on rectangular grids.

and ⎛ϕ ϕ ϕ ⎞ ⎛ ω1 ⎞ ⎜ 2 3 4 ⎟ ⎜ω ⎟ ⎜ ϕ ϕ ϕ ⎟ ⎜ 2⎟ = ⎜ 1 3 4 ⎟ ⎜ ω3 ⎟ ⎜ ϕ ϕ ϕ ⎟ ⎜ ⎟ ⎜ 1 2 4⎟ ⎝ ω4 ⎠ ⎜ ϕ ϕ ϕ ⎟ ⎝ 1 2 3⎠

■

NUMERICAL CASE STUDIES In this section, numerical case studies of the one-dimensional gas−solid reaction model and two-dimensional Keller−Segel model are presented. One-Dimensional Test Problems for Gas−Sold Reaction Model. Here, three case studies are considered, and the numerical results of CE/SE method are compared with the staggered central scheme.22 The initial data of all test

(64)

Using the above relations, we finally obtain 4

(wα)in, j =

∑i = 1 u i ,α|ωi|θ 4

∑i = 1 |ωi|θ

, α ∈ {x ,y}

(65) 9182



Article

Figure 10. Results of problem 6 at t = 1.5 × 10−5, t = 4.5 × 10−5 and t = 3.0 × 10−4 (from top to bottom). 2 ⎧ ⎫1/2 ⎪ [∫ (φ − φ ) dx ] ⎪ ref ⎬ Error = ⎨ ⎪ ⎪ [∫ φ 2dx] ⎩ ⎭

problems are given in Table 1. In all test problems, the computational domain 0 ≤ x ≤ 200 is discretized into 500 grid points. The first two test problems were also considered by Hassanzadeh et al.2 Figure 2 shows errors in the schemes at different numbers of grid points. Figure 3 gives the results of problem 1 at td = 0.6. The numerical results of both schemes are in good agreement with each other. However, CE/SE method has better resolved the peak in the temperature profile as compared to the central scheme. The numerical error is calculated by using the following expression

(66)

where φ can be either temperature or concentration and the subscript ref denotes reference solution. The reference solution was obtained at 5 × 103 grid points by using CE/SE method. Table 2 shows a comparison of numerical errors at different numbers of grid points and CPU times at N = 200. In this table, CE/SE method is compared with the Superbee, Koren, central, and first-order upwind schemes.25,26 It can be observed that the proposed CE/SE method has produced less errors as compared 9183



Article

method captures the height of density profile much better than the central scheme. The numerical results are comparable to those obtained by Epshteyn.4

to other schemes at low computational cost. It can be observed that in the CE/SE method errors reduce faster on increasing the numbers of grid points. The numerical results are in good agreement with available results in the literature.2 The results of problem 2 at td = 2.5 are shown in Figure 4. The reaction peak in the gas concentration obtained by CE/SE method has larger height than those of central scheme and other methods presented in the article of Hassanzadeh et al.2 Similarly, results of problem 3 at td = 1.0 are given in Figure 5. An overshoot can be observed in the temperature profile of the central scheme which is very small in the results of CE/SE method. In all test problems, the results of both schemes show good agreement. However, the CE/SE method has better resolved the sharp discontinuities and peaks. Two-Dimensional Test Problems for KS Model. In this section, three two-dimensional test problems are considered for the current KS-model. Problem 4. The initial data inside a square region [(−1/2), (1/2)] × [(−1/2),(1/2)] with ξ = 1 are taken as3 2

ρ(0,x ,y) = 103exp−10 (x

2

+ y2 )

= 5 × 102exp−50(x

2

■

CONCLUSIONS AND REMARKS The space time CE/SE method was implemented for solving one- and two-dimensional convection−diffusion−reaction type partial differential equations. The major difficulties associated with the current systems include the formation of steep reaction fronts and blowup of the solution in a finite time. The traditional numerical approaches are not capable to capture such rapid variation of the physical quantities. The CE/SE was successfully implemented in one- and two-space dimensions, and numerical results of the method were found in good agreement with those available in the literature and with central schemes on staggered grids. The staggered central schemes were also implemented for the first time to these models. Both methods captured the rapid variations in the solutions more efficiently and accurately. However, CE/SE method performed better for the solutions with sharp discontinuities and narrow peaks. The proposed method produced less errors in the results as compared to other available methods.

, c(0,x ,y) + y2 )

■

(67)

A uniform grid of 100 × 100 points with outflow boundary conditions is used. The cell densities computed at t = 1 × 10−4 are displayed in Figure 6. For comparison, the same problem was also solved by the central scheme,23 and the results of both schemes are compared through the one-dimensional plots along y = 0.0 in Figure 6. Further, the cell densities are computed at different times, and their comparison along y = 0.0 are presented in Figure 7. The blowup of density forming a narrow spike is visible in Figures 6 and 7. The figures show that density peaks become narrow and longer with the passage of time. Moreover, it can be observed that both schemes give comparable results. However, CE/SE method resolves the peaks better. The numerical results are also in good agreement with the available results in the literature and the peaks heights are comparable.3 Problem 5. In this problem, ξ = 1 and the same square region [(−1/2),(1/2)] × [(−1/2),(1/2)] is considered. The initial data are 2

ρ(0,x ,y) = 103exp−10 (x

2

+ y2 )

, c(0,x ,y) = 0

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This work was partially supported by the Higher Education Commission (HEC) of Pakistan.

2

+ y2 )

, c(0,x ,y) = 420exp−42(x

(68)

2

REFERENCES

(1) Alhumaizi, K. Flux-limiting solution techniques for simulation of reaction diffusion−convection system. Commun. Nonlinear Sci. Num. Simul. 2007, 12, 953−965. (2) Hassanzadeh, H.; Adebi, J.; Darvish, M. P. A comparative study of flux limiting methods for numerical simulation of gas solid reaction with arrhenius type kinetic. Comput. Chem. Eng. 2009, 33, 133−134. (3) Chertock, A.; Kurganov, A. A second-order positivity preserving central-upwind scheme for chemotaxis and haptotaxis models. Numer. Math. 2008, 111, 169−205. (4) Epshteyn, Y. Discontinuous Galerkin methods for the chemotaxis and haptotaxis models. J. Comput. Appl. Math. Phys. 2009, 224, 168− 181. (5) Hortsmann, D. From 1970 until now: The Kellar−Segal model in chemotaxis and its consequence I. Jahresber DMV 2003, 105, 103−165. (6) Hortsmann, D. From 1970 until now: The Kellar−Segal model in chemotaxis and its consequence II. Jahresber DMV 2004, 106, 51−69. (7) Keller, E. F.; Segel, L. A. Model for chemotaxis. J. Theor. Biol. 1971, 30, 225−234. (8) Alhumaizi, K.; Henda, H.; Soliman, M. Numerical analysis of reaction−diffusion−convection system. Comput. Chem. Eng. 2003, 27, 579−594. (9) Filbet, F. A finite volume scheme for the Patlak−Keller−Segel chemotaxis model. Numer. Math. 2006, 104, 457−488. (10) Loh, C. Y.; Zaman, K. B. M. Q. Numerical investigation of transonic resonance with a convergent−divergent nozzle. AIAA J. 2002, 40, 2393−2401. (11) Tyson, R.; Stern, L. G.; LeVeque, R. J. Fractional step methods applied to a chemotaxis model. J. Math. Biol. 2000, 41, 455−475. (12) Chang, S. C. The method of space time conservation element and solution element: A new approach for solving the Navier Stokes and Euler equations. J. Comput. Phys. 1995, 119, 295−234.

Figure 8 shows the results at different times. In this problem, the density ρ blows up at later time than in problem 4. Resultantly, the spreading of cell density occur and the maximum value is lesser than that in problem 4. The figures show the spreading of density profiles and the hight of peaks reduces with the passage of time. The number of grid points is the same as in problem 4. Once again, the results of CE/SE method are superior to the central scheme. Problem 6. This problem considers more initial data with ξ = 1 as given below ρ(0,x ,y) = 840exp−84(x

AUTHOR INFORMATION

+ y2 )

(69) −4

Figure 9 shows the results at t = 3 × 10 . In this problem, ρ and c blowup at the origin in a finite time. A fixed mesh with 100 × 100 points is considered. Both schemes have comparable results, but CE/SE method has better resolved the density peak. This can be further observed in the density plots obtained at different times in Figure 10. The figures show that density peaks become narrow and longer with the passage of time. The CE/SE 9184



Article

(13) Zhang, M.; Yu, S. T.; Chang, S. C. A space−time conservation element and solution element method for solving the two-dimensional un-steady Euler equations using quadrilateral and hexhedral meshes. J. Comput. Phys. 2000, 175, 168−199. (14) Chang, S. C.; Wang, X. Y.; Chow, C. Y. The space−time conservation element and solution element method. A new high resolution and genuinely multidimensional paradigm for solving conservation laws. J. Comput. Phys. 1999, 156, 89−136. (15) Chang, S. C.; Wang, X. Y.; Chow, C. Y. New Developments in the Method of Space−time Conservation Element and Solution Element: Applications to Two-Dimensional Time-Marching Problems; NASA TM 106758, 1994. (16) Liu, M.; Wang, J. B.; Wu, K.-Q. The direct aero-acoustics simulation of flow around a square cylinder using the CE/SE scheme. J. Algorithms Comput. Technol. 2009, 1, 525−537. (17) Chang, S. C.; Wang, X. Y.; To, W. M. Application of the spacetime conservation element and solution element method to onedimensional convection−diffusion problems. J. Comput. Phys. 2000, 165, 189−215. (18) Liu, N. S.; Chen, K. H. Noise Computation of a Shock-Containing Supersonic Axisymmetric Jet by The CE/SE Method; AIAA Paper 2000− 0475; 38th AIAA Aerospace Sciences Meeting, January 10−13, Reno, NV, 2000. (19) Liu, N. S.; Chen, K. H. Wave computation in incompressible flow using the space−time conservation element and solution element method. AIAA J. 2001, 39, 794−801. (20) Qamar, S.; Warnecke, G. Application of space−time CE/SE method to shallow water magnetohydrodynamics equations. J. Comput. Appl. Math. 2006, 196, 132−149. (21) Qamar, S.; Mudasser, S. On the application of a variant CE/SE method for solving two-dimensional MHD equations. J. Appl. Numer. Math. 2010, 60, 587−696. (22) Nessyahu, H.; Tadmor, E. Non-oscillatory Central differencing for hyperbolic conservation laws. J. Comput. Phys. 1990, 87, 408−463. (23) Jaing, G.-S.; Tadmor, E. Non-oscillatory central schemes for multidimensional hyperbolic conservation laws. SIAM J. Sci. Comput. 1998, 19, 1892−1917. (24) Rajaiah, J.; Dandekar, H.; Hlavacek, V. Modeling of exothermic gas−solid noncatalytic reactions. Ind. Eng. Chem. Res. 1987, 26, 1424− 1434. (25) Koren, B. A Robust Upwind Discretization Method for Advection, Diffusion and Source Terms. In Numerical Methods for Advection−Diffusion Problems; Vreugdenhil, C. B., Koren, B., Eds.; Notes on Numerical Fluid Mechanics 45; Vieweg Verlag: Braunschweig, 1993; pp 117−138. (26) Leveque, R. J. Finite Volume Methods for Hyperbolic Systems; Cambridge University Press: Cambridge, 2003.

9185


SE Method for ... - ACS Publications

Recommend Documents