- Email: [email protected]
A finite volume difference scheme for a model of settling particle dispersion from an elevated source in an open-channel flow
Computers and Mathematics with Applications (
)
–
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
A finite volume difference scheme for a model of settling particle dispersion from an elevated source in an open-channel flow T. Chernogorova a,∗ , L. Vulkov b a
FMI, University of Sofia, 5, J. Bourchier Blvd., 1164 Sofia, Bulgaria
b
FNSE, University of Rousse, 8, Studentska str., 7017 Rousse, Bulgaria
article
info
Article history: Available online xxxx Keywords: Dispersion Settling velocity Degenerate parabolic equation Finite volume difference scheme Conservation of non-negativity
abstract A finite volume fitted difference scheme is constructed to solve the unsteady convectivediffusion equation transformed on a finite domain, modeling longitudinal dispersion of suspended particles with settling velocity in a turbulent shear flow over a rough-bed surface. First we discuss the well-posedness of the differential problem and the non-negativity of its solution. Then, to overcome the degeneracy at the part of the boundary, starting from the divergent form of the equation we perform a local fitted space discretization. This approximation is determined by a set of two-point boundary value problems. Non-negativity of the numerical concentration of suspended fine particles is proved. Some results from computational experiments are presented to illustrate the properties of the constructed scheme. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction We discuss a model due to Mondal and Mazumder [1–3], see also [4], for a dispersion of fine particles in an oscillatory turbulent flow. The model is based on the time-dependent advection–diffusion equation
∂C ∂C ∂ 2C ∂ ∂C ∂C + u(z ) − ws = kx (z ) 2 + kz (z ) + f (x, z ) ∂t ∂x ∂z ∂x ∂x ∂z on (−∞, ∞) × (z 0 , 1) and subjected to boundary and initial conditions ∂C C (±∞, z , t ) = 0, kz (z ) + ws C = 0, C (x, z , 0) = 0, ∂z z = z 0 ,1
(1)
(2)
where z 0 is the boundary roughness height, ws is the setting velocity and f (x, z ) = δ(x)δ(z − zp ), i.e., a point source at (x, z ) = (0, zp ). Here the problem has been expressed in terms of dimensionless variables and it is used for fully developed homogeneous turbulent flow [4,1,2]. In our present study, the velocity distribution is 1 z u(z ) = ln 0 + W (z ), k z where k > 0 is the von-Karman constant, z 0 is the equivalent bed roughness, W (z ) = 2kΠ sin2 π2 z , and Π is the wakestrength parameter.
∗
Corresponding author. Tel.: +359 897846436. E-mail addresses: [email protected] (T. Chernogorova), [email protected] (L. Vulkov).
http://dx.doi.org/10.1016/j.camwa.2014.03.012 0898-1221/© 2014 Elsevier Ltd. All rights reserved.
2
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
When the Reynolds decomposition is performed on the laminar advection–diffusion equation by decomposing the velocity (u = u + u′ ) and concentration (C = C + C ′ ) into the sum of their mean and fluctuating parts, the term (∂ u′ C ′ /∂ xj ) represents the transport of concentration C due to turbulent fluctuations. Because the molecular diffusivity (kmd ) is usually a very small quantity, u′ C ′ ≫ kmd ∂ C /∂ xj , the effects for molecular diffusion are often negligible compared to the effects of turbulence, and so omitted from the model. However, in this model the molecular diffusion is an important mechanism for mixing of the concentration to the flow at the smallest scales, and so the authors of [2,3] take the diffusivity profile as kx (x) = kmd + ked (z ), where ked (z ) is the eddy-diffusivity as proposed in [5,6]: ked (z ) = k(1 − z )
1 z
−1 + Π π sin(π z ) .
Since our problem is posed in a fully developed homogeneous turbulent flow, we take kx (z ) = kz (z ) and neglect the crossstream diffusion terms. The studies of dispersion phenomena of settling particles released from an elevated continuous line-source in a turbulent open-channel flow are important to understand the process of contamination in the environment [1–3,5]. Mazumder and Bandyopadhyay [1] presented a numerical solution to the convection–diffusion equation to study the stream-wise dispersion of contaminant due to turbulent shear flow over a gravel-bed surface, employing a modified ‘‘Logwake law’’ and corresponding diffusivity due to Nezu and Nakagawa [5]. Their method provides a more general numerical scheme for the unsteady convection–diffusion equation than the theoretical study of Sullivan and Yip [6] based on the small time asymptotic solution. A combined scheme of central difference for the diffusion terms and a four-point upwind scheme for convective terms are used to solve the convection–diffusion equation for the steady problem. The alternating direction implicit (ADI) method is adopted to generate the numerical solution of a two-dimensional unsteady convection–diffusion equation and is developed in [2,3]. The authors use this scheme for numerical simulations to examine the stream-wise dispersion of suspended fine particles with settling velocities in oscillatory turbulent shear flow when these particles are being released from an elevated source. They, in order to incorporate the zero Dirichlet boundary condition of (2), transform the unbounded region to a bounded one. We also will work with transformed problem, namely (3a)–(3d), in the next section. However, as we shall see, the features apparent in the model problem (3a)–(3d) present difficulties which require new numerical approaches. These difficulties include the degeneracy of the elliptic part of the parabolic equation (3a), noncoercivity of the bilinear elliptic form, corner points of the domain, etc. We provide a theoretical analysis in appropriate weighted Sobolev spaces of the continuous problem on bounded domain. Here we propose a new approach for the construction of an efficient difference scheme. This is a fitted technique based on the old idea of Allen and Southwell well described in [7] and that was extended to different problems by several authors [8–10]. In a previous work of the first author [8] numerical experiments show higher accuracy of the present construction in comparison with the other known difference schemes near the degeneracy. The rest of the paper is organized as follows. In Section 2 we present the differential problem on bounded domain and discuss our basic assumptions and some properties of the solution. The space discretization is developed in Section 3. Section 4 is devoted to the system matrix after the time discretization. We show that the maximum principle is satisfied and thus the discrete solution is non-negative in correspondence to the continuous solution. In Section 5 we present some results from computational experiments. In Section 6 we formulate some conclusions. 2. The differential problem on bounded domain To compute an accurate numerical solution to the model given above, several issues need to be addressed. First, the problem is posed on a semi-infinite domain. Numerical solutions are possible only on finite domains, so the domain must either be truncated or transformed. Here we take the approach, discussed in [4,1,2], since it has proved successful in related studies. In [2,3] a tanh-transformation from domain (x, z ) ∈ (−∞, ∞) × [z 0 , 1] to (ξ , z ) ∈ (−1, 1) × [z 0 , 1] is made using x=
1 2a
log
1+ξ 1−ξ
,
where a is a ‘‘stretching’’ factor, usually chosen based on the computational experience. Following the transformation the model becomes
∂C ∂C ∂ 2C ∂C + u(z )a(1 − ξ 2 ) − a2 kξ (z )(1 − ξ 2 ) (1 − ξ 2 ) 2 − 2ξ ∂t ∂ξ ∂ξ ∂ξ ∂ ∂C ∂C = kz (z ) + ws + δ(ξ )δ(z − zp ), ∂z ∂z ∂z l± (C ) ≡ C (±1, z , t ) = 0,
(3a) (3b)
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
lz 0 ,1 (C ) ≡ kz (z )
∂C + ws C ∂z
)
–
3
= 0,
z = z 0 ,1
(3c)
C (ξ , z , 0) = C0 (ξ , z ) = 0.
(3d)
The parabolic equation (3a) belongs to the second order PDEs with the non-negative characteristic form. The main character of such kind of equations is the degeneracy. It can be easily seen that at ξ = ∓1 Eq. (3a) degenerates into two parabolic one-dimensional equations:
∂ ∂C = ∂t ∂z
kz (z )
∂C ∂z
+ ws
∂C , ∂z
(z , t ) ∈ (z 0 , 1) × (0, T ].
(4)
By the well-known Fichera’s theory (see [11]) for degenerate parabolic equations, we have that at the degenerate boundaries ξ = −1 and ξ = 1, the boundary conditions should not be given. Therefore, the parabolic problem (3a), (4), (3c) is well defined. But from physical motivation we have imposed the boundary conditions (3b) and the initial condition (3d). It is easy to check that (4) follows from (3b), (3d). Therefore the boundary conditions (3b) are the particular case of (4). As discussed in [12], the Poisson equation with the Dirac-delta function on the right-hand side does not have an H 1 (Ω )solution. For this and also for the reason of computations (see Section 5) we will assume that the right hand side is already regularized and at least the regularized function belongs to L2 (Ω ), for example, see [13]. Also, this will simplify the next theoretical considerations. We rewrite Eq. (3a) in divergent form:
∂C − LC + p(ξ , z )C = fR (ξ , z ) for (ξ , z , t ) ∈ QT = Ω × (0, T ], ∂t
(5)
where LC = ∇ · (A∇ C + bC ) ,
A=
a11 0
0 a22
,
Ω = (−1, 1) × (z 0 , 1),
a11 = a2 kξ (z )(1 − ξ 2 )2 ,
a22 = kz (z ),
b1 (ξ , z ) (1 − ξ 2 )a 2ξ akξ (z ) − u(z ) , = b2 ws
p(ξ , z ) = 2a2 kξ (z )(1 − 3ξ 2 ) + 2ξ au(z ), and fR (ξ , z ) is a regularization of the delta function δ(ξ )δ(z − zp ). Further, to handle the degeneracy in Eq. (5), we introduce the weighted inner product and corresponding norm on L2,w (Ω ) by
(u, v)w :=
Ω
(1 − ξ 2 )2 uv dξ dz ,
∥v∥0,w =
(v, v)w =
Ω
(1 − ξ 2 )2 v 2 dξ dz
1/2
.
The space of all weighted square-integrable function is defined as L2,w (Ω ) := {v : ∥v∥0,w < ∞}. Using a standard argument it is easy to show that the pair (L2,w (Ω ), (·, ·)w ) is a Hilbert space (cf., for example [14]). Using L2 (Ω ) and L2,w (Ω ) we define the following weighted Sobolev space 1 Hw (Ω ) :=
∂v ∂v ∈ L2,w (Ω ) and ∈ L2 ( Ω ) v ∈ L2 (Ω ), ∂ξ ∂z
1 with the corresponding inner product. It is also easy to prove that Hw (Ω ) is a Hilbert space with the norm
∥v∥1,w :=
Ω
v + (1 − ξ ) 2
2 2
∂v ∂ξ
2
+
∂v ∂z
2
1/2 dξ dz
.
1 (Ω ) define the bilinear form Next, for C , η ∈ Hw
A(C , η) ≡
Ω
a11
∂ C ∂η ∂ C ∂η ∂η ∂η + a22 + b1 C + b2 C + pC η dξ dz . ∂ξ ∂ξ ∂z ∂z ∂ξ ∂z
Using this weighted Sobolev space we define the following variational problem corresponding to (5) and (3b)–(3d): find 1 C ∈ Hw (Ω ) satisfying the initial condition (3d) such that for all η ∈ Hw1 (Ω )
Ω
∂C ηdξ dz + A(C , η) = ∂t
Ω
fR ηdξ dz ,
a.e. in (0, T ].
1 Theorem 1. There exists a unique solution C ∈ Hw (Ω ) to problem (5), (3b)–(3d).
4
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
Proof. By changing the variables C := eρ t C , ρ = const > 0, we can replace the coefficient p in (5) by p + ρ (let us note that we will use C for the new function as well). Therefore, we will assume that
∂C − LC + (p + ρ)C = e−ρ t fR (ξ , z ), ∂t ∂ C ∂η ∂ C ∂η ∂η ∂η + a22 + b1 C + b2 C + (p + ρ)C η dξ dz . Aρ (C , η; t ) = a11 ∂ξ ∂ξ ∂z ∂z ∂ξ ∂z Ω For any ε > 0, there exists a constant K (ε) > 0 such that b1 C
2 2 ∂C ∂C ∂C ∂C ∂C ∂C + b2 C ≥ − b1 C + b2 C ≥ −ε(1 − ξ 2 )2 −ε − K (ε)C 2 . ∂ξ ∂z ∂ξ ∂z ∂ξ ∂z
Therefore, for sufficiently small ε and large ρ
2 2 ∂C ∂C 2 2 2 Aρ (C , C ; t ) ≥ (a kξ (z ) − ε)(1 − ξ ) + (kz (z ) − ε) ∂ξ ∂z Ω + (p + ρ − K (ε))C 2 dξ dz ≥ K (ε, ρ)∥C ∥2H 1 (Ω ) , w
where K (ε, ρ) > 0 is a constant. Now, in a standard way, using the Cauchy–Schwarz inequality one can prove that the bilinear form Aρ (C , η; t ) is continuous. Finally by Lemma 1 and Theorem 1.33 of [15], we have that the problem (5), (3b)–(3d) has a unique solution and thus we complete the proof. Considering the settling particles dispersion, the concentration C (x, z , t ) of the suspended particles cannot be negative if the concentration was non-negative in the initial state and on the boundary of the channel and in addition, the unity concentration is maintained at the steady line-source (x = 0, z = zp ) at all times. This property is called non-negativity preservation (NP). Conservation of this property is very important at the adequate discretizations of PDEs [16]. Further we analyze the NP for problem (3a)–(3d). 1 For a function v in Hw (Ω ) we denote the positive and negative parts of v by v + and v − . That is v = v + + v − ; v + ≥ 0 − and v ≤ 0. By [17, Lemma 7.6, p. 150], we have that Dv + =
Dv, 0,
if v > 0 , if v ≤ 0
Dv − =
Dv, 0,
if v < 0 , if v ≥ 0
∂v or Dv = Dz v = ∂v . Also v + v − = Dv + Dv − = Dv + v − = v + Dv − = 0, (ξ , z ) ∈ Ω . where Dv = Dξ v = ∂ξ ∂z We consider the following inequalities
∂C − LC + p(ξ , z )C ≥ 0 in QT , ∂t l± (C ) ≥ 0
in (z 0 , 1) × (0, T ]
(6)
and lz 0 ,1 (C ) ≥ 0
in (−1, 1) × (0, T ].
(7)
Formally multiplication (6) by a non-negative smooth function η and integration by parts lead to the following definition. 1 ,0 Let V2 (QT ) (see [18]) be the Banach space consisting of all functions in Hw (QT ) = C ((0, T ); Hw1 (Ω )) having finite norm
∥v∥V2 (QT ) = sup ∥v(·, ·; t )∥L2 (Ω ) + ∥∇v∥L2,w (QT ) , 0 ≤t ≤T
where T
∥∇v∥L2,w (QT ) =
0
∥v∥2H 1 (Ω ) (t )dt
1/2
w
.
1,1 (QT ) = Following [18], a function v ∈ V2 (QT ) is said to satisfy weakly (6) and (7) if for any non-negative η ∈ Hw 1 H 1 ((0, T ); Hw (Ω ))
Ω
C (ξ , z , t )η(ξ , z , t )dξ dz −
Qt
C ηt dξ dzdτ +
t
A(C , η)dτ ≥ 0
(8)
0
for almost every t ∈ (0, T ]. Here Qt = Ω × (0, t ). It is easy to see that if C ∈ C 2,1 (Q T ) satisfies (6), (7) point-wise then (8) holds for C . Conversely, if C ∈ V2 (QT ) and satisfies (8) and C is sufficiently smooth then it also verifies (6), (7) point-wise. We will now prove the following assertion. Theorem 2. Let C (ξ , z , t ) be in V2 (QT ) and satisfy (6), (7). Then C (ξ , z , t ) ≥ 0 in QT .
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
5
Proof. In fact, we will work again with Aρ (C , η). Using the Steklov average and taking the limit (see [18], Chapter 3) we can formally take η = −C − ≥ 0 in (8) to obtain 1
2
Ω
(C − (ξ , z , t ))2 dξ dz
+ Qt
∂C ∂C− ∂C ∂C− ∂C− ∂C− − a11 + a22 + b1 C + b2 C + (p + ρ)CC dξ dzdτ ≤ 0. ∂ξ ∂ξ ∂z ∂z ∂ξ ∂z
(9)
We have DC .DC − = (DC − )2 , CC − = (C − )2 . Let us introduce Kz = minz 0 ≤z ≤1 kz (z ), A11 = a2 Kz . Next, using the estimate in the proof of Theorem 1 for the expression b1 C with ε =
∂C− ∂C− + b2 C ∂ξ ∂z 1 2
min(Kz , A11 ), we obtain from (9)
− 2 − 2 ∂C ∂C 2 2 (C (ξ , z , t )) dzdξ + min(Kz , A11 ) + (1 − ξ ) dξ dzdτ 2 Ω ∂ξ ∂z Qt ≤ [K (ε) − (p + ρ)](C − )2 dξ dzdτ . 1
−
2
Qt
Now, we choose ρ sufficiently large such that K (ε) − (p + ρ) ≤ 0. The above yields 1 2
(C (ξ , z , t )) dξ dz + min(Kz , A11 ) −
Ω
2
(1 − ξ )
2 2
Qt
for a. e. t ∈ (0, T ]. This implies that C (ξ , z , t ) ≥ 0 for a. e. t ∈ (0, T ].
1 2
Ω
∂C− ∂ξ
2
+
∂C ∂z
2
dξ dzdτ ≤ 0,
(C − (ξ , z , t ))2 dξ dz = 0 and therefore C − (ξ , z , t ) ≡ 0. We then conclude that
3. Space discretization In this section we propose a finite volume method with a fitting technique for the numerical solution of problem (5), (3b)–(3d). We have chosen this method due to the fact that it takes into account the degeneracy of the equation and for the problem investigated in [8] gives more accurate results near the boundaries of degeneration compared to the commonly used finite difference methods. We introduce the meshes
−1 = ξ1 < ξ2 < · · · < ξNξ +1 = 1,
z 0 = z1 < z2 < · · · < zNz +1 = 1
and put hξi = ξi+1 − ξi , i = 1, 2, . . . , Nξ , hzj = zj+1 − zj , j = 1, 2, . . . , Nz . We also introduce the auxiliary meshes
ξi−1/2 =
ξi−1 + ξi 2
,
i = 2, 3, . . . , Nξ + 1,
zj−1/2 =
zj−1 + zj 2
,
j = 2, 3, . . . , Nz + 1
and let ξ1/2 = ξ1 = −1, ξNξ +3/2 = ξNξ +1 = 1, z1/2 = z1 = z 0 , zNz +3/2 = zNz +1 = 1, mi = ξi+1/2 − ξi−1/2 , lj = zj+1/2 − zj−1/2 , Li,j = mi lj , i = 1, 2, . . . , Nξ + 1, j = 1, 2, . . . , Nz + 1. A. Internal grid nodes We integrate Eq. (5) on the cell ℜi,j = [ξi−1/2 , ξi+1/2 ] × [zj−1/2 , zj+1/2 ], i = 2, 3, . . . , Nξ ; j = 2, 3, . . . , Nz :
ℜi,j
∂C dξ dz = ∂t
∇ · (k(C ))dξ dz −
ℜi,j
pCdξ dz +
ℜi,j
fR dξ dz ,
(10)
ℜi,j
where k(C ) = A∇ C + bC . For all integrals with the exception of the second one in (10) we use the quadrature formula of the central rectangles, to get
∂ Ci,j Li,j = ∂t
∇ · (k(C ))dξ dz − Li,j pi,j Ci,j + Li,j fR,i,j ,
(11)
ℜi,j
where
∂ C ∂ Ci,j = , ∂ t (ξi ,zj ,t ) ∂t
p|(ξi ,zj ,t ) = pi,j ,
C |(ξi ,zj ,t ) = Ci,j ,
fR |(ξi ,zj ,t ) = fR,i,j .
6
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
Using the Gauss–Ostrogradski formula we find for the second integral in (10)
∇ · (k(C ))dξ dz =
ℜi,j
k · n ds =
(ξi−1/2 ,zj+1/2 )
− (ξi−1/2 ,zj−1/2 ) (ξi+1/2 ,zj−1/2 )
∂C + b1 C dz ∂ξ (ξi+1/2 ,zj−1/2 ) (ξi+1/2 ,zj+1/2 ) ∂C ∂C a22 + b1 C dz + + b2 C dξ a11 ∂ξ ∂z (ξi−1/2 ,zj+1/2 ) ∂C a22 + b2 C dξ = I1 − I2 + I3 − I4 . ∂z
a11
∂ℜi,j
(ξi+1/2 ,zj+1/2 )
− (ξi−1/2 ,zj−1/2 )
(12)
Let us consider the first integral in (12). From the expressions for a11 and b1 we get a11
∂C + b1 C = (1 − ξ 2 )ρ(C ), ∂ξ
where the flux is
ρ(C ) = a2 kξ (z )(1 − ξ 2 )
∂C + a 2ξ akξ (z ) − u(z ) C . ∂ξ
Further we will follow the basic idea of [10]. To approximate the flux ρ(C ) we solve the boundary value problem
′
ci+1/2,j (1 − ξ 2 )v ′ + bi+1/2,j v ξ = 0,
(13a)
v(ξi , zj ) = Ci,j ,
v(ξi+1 , zj ) = Ci+1,j , (13b) where c = a kξ (z ), b = a 2ξ akξ (z ) − u(z ) , ci+1/2,j = c (ξi+1/2 , zj ), bi+1/2,j = b(ξi+1/2 , zj ) and Ci,j , Ci+1,j are the values of the approximate solution C at the nodes (ξi , zj ) and (ξi+1 , zj ). Let us note that the coefficient c does not depend on ξ and in this case we have ci+1/2,j = c (zj ). We introduce the notations 2
αi,j =
bi+1/2,j ci+1/2,j
,
∆i,j (ξi ) =
1 + ξi
α2i,j
1 − ξi
.
Then for i = 2, 3, . . . , Nξ − 1, j = 2, 3, . . . , Nz we obtain I1 ≈ lj (1 − ξ 2 )ρ(C ) (ξ
i+1/2, zj )
= lj (1 − ξi2+1/2 )bi+1/2,j
∆i,j (ξi+1 ) Ci+1,j − ∆i,j (ξi ) Ci,j . ∆i,j (ξi+1 ) − ∆i,j (ξi )
(14)
As noted in [10] for i = Nξ the above analysis does not apply to the approximation of the flux because (13a) is degenerate and this leads to the solution of (13a) which can never satisfy both conditions in (13b) unless CNξ ,j = CNξ +1,j . To solve this difficulty let us re-consider (13a)–(13b) with an extra degree of freedom, solving the problem c¯Nξ +1/2,j (1 − ξ )v ′ + bNξ +1/2,j v
v(ξNξ , zj ) = CNξ ,j ,
′
ξ
= M2 ,
v(ξNξ +1 , zj ) = CNξ +1,j = 0,
where c¯ = (1 + ξ )a2 kξ (z ) and find for the first integral in (12) the approximation I1 ≈ lj (1 − ξN2ξ +1/2 )
1 CNξ +1,j c¯Nξ +1/2,j + bNξ +1/2,j − CNξ ,j c¯Nξ +1/2,j − bNξ +1/2,j . 2
(15)
By analogy with the investigations above we approximate the second integral in (12) and find for j = 2, 3, . . . , Nz – for 3 ≤ i ≤ Nξ :
I2 ≈ lj (1 − ξ 2 )ρ(C ) (ξ
i−1/2, zj )
∆i−1,j (ξi ) Ci,j − ∆i−1,j (ξi−1 ) Ci−1,j ; ∆i−1,j (ξi ) − ∆i−1,j (ξi−1 )
(16)
1 c˜3/2,j + b3/2,j C2,j − c˜3/2,j − b3/2,j C1,j ,
(17)
= lj (1 − ξi2−1/2 )bi−1/2,j
– for i = 2: I2 ≈ lj (1 − ξ32/2 ) ρ(C )|(ξ3/2, zj ) = lj (1 − ξ32/2 ) where c˜ = (1 − ξ )a2 kξ (z ).
2
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
7
The problem is non-degenerate with respect to the variable z and for the third integral in (12), following [19], we apply a classical approximation: (ξi+1/2 ,zj+1/2 )
∂C ∂C + b2 C dξ ≈ kz (z ) + ws C mi ∂z ∂z (ξi−1/2 ,zj+1/2 ) (ξi ,zj+1/2 ) ws di,j+1 ws di,j+1 + − Ci,j+1 − Ci,j , = mi
I3 =
a22
hzj
2
hzj
2
(18)
where
di,j+1 =
zj+1
1 hzj
−1
1 kz (z )
zj
≈ kz (zj+1/2 ).
dz
Analogically for the fourth integral in (12) we get
I4 ≈ m i
di,j
+
hzj−1
ws 2
C i ,j −
di,j hzj−1
−
ws
Ci,j−1 .
2
(19)
Replacing in (12) the obtained approximations for the integrals from (14)–(19) and then replacing the derived expression in (11), we obtain a semi-discretization of Eq. (5) in the internal grid points for i = 2, 3, . . . , Nξ , j = 2, 3, . . . , Nz . B. Boundary grid nodes For the left vertical boundary ξ = −1, z ∈ [z 0 , 1] and for the right vertical one ξ = 1, z ∈ [z 0 , 1] we have respectively C1,j = 0,
CNξ +1,j = 0,
j = 1, 2, . . . , Nz + 1.
For the down horizontal boundary ξ ∈ (−1, 1), z = z 0 we proceed as for the internal grid nodes, but integrating Eq. (5) in the region ℜi,1 = [ξi−1/2 , ξi+1/2 ] × [z1/2 , z3/2 ], i = 2, 3, . . . , Nξ (i.e., in the semi-interval by z). As in (11) we get
∂ Ci,1 L i ,1 = ∂t
∇ · (k(C ))dξ dz − Li,1 pi,1 Ci,1 + Li,1 fR,i,1 .
(20)
ℜi,1
For the integral in (20) we have (ξi+1/2 ,z3/2 )
∂C + b1 C dz ∇ · (k(C ))dξ dz = k · n ds = a11 ∂ξ ℜi,1 ∂ℜi,1 (ξi+1/2 ,z1 ) (ξi−1/2 ,z3/2 ) (ξi+1/2 ,z3/2 ) ∂C ∂C − a11 + b1 C d ξ + a22 + b2 C d ξ ∂ξ ∂z (ξi−1/2 ,z1 ) (ξi−1/2 ,z3/2 ) (ξi+1/2 ,z1 ) ∂C − a22 + b2 C dξ = I1d − I2d + I3d − I4d . ∂z (ξi−1/2 ,z1 )
(21)
We find approximations for I1d , I2d and I3d by proceeding analogically as for the internal grid nodes. In fact these approximations can be obtained from (14) to (18) for j = 1. For I4d we have I4d = 0. For the upper horizontal boundary we accomplish analogical operations integrating Eq. (5) in the region ℜi,Nz +1 = [ξi−1/2 , ξi+1/2 ] × [zNz +1/2 , zNz +3/2 ], i = 2, 3, . . . , Nξ . Using all approximations for the internal and boundary grid nodes we obtain the following semi-discrete problem (note that Ci,j = Ci,j (t ) for i = 1, 2, . . . , Nξ + 1, j = 1, 2, . . . , Nz + 1):
∂ Ci,j i ,j i,j i ,j i,j i ,j Li,j = −ei,j Ci,j + ei+1,j Ci+1,j + ei−1,j Ci−1,j + ei,j+1 Ci,j+1 + ei,j−1 Ci,j−1 + Li,j fR,i,j , ∂t i = 2, 3, . . . , Nξ , j = 2, 3, . . . , Nz ; ∂ Ci,1 i ,1 i,1 i,1 i,1 Li,1 = −ei,1 Ci,1 + ei+1,1 Ci+1,1 + ei−1,1 Ci−1,1 + ei,2 Ci,2 + Li,1 fR,i,1 , ∂t ∂ Ci,Nz +1 i ,N + 1 i,N +1 Li,Nz +1 = −ei,Nzz +1 Ci,Nz +1 + ei+1z,Nz +1 Ci+1,Nz +1 ∂t + eii,−N1z,+N1z +1 Ci−1,Nz +1 + eii,,NNzz +1 Ci,Nz + Li,Nz +1 fR,i,Nz +1 , i = 2, 3, . . . , Nξ ; C1,j = 0, CNξ +1,j = 0, j = 1, 2, . . . , Nz + 1.
(22)
8
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
The expressions for the coefficients of the system (22) are given below: 2,1 e2,1
=
l1 (1 − ξ52/2 )b5/2,1 ∆2,1 (ξ2 )
l1
(1 − ξ
2 3/2
) c˜3/2,1 + b3/2,1 + m2
d2,2
+ L2,1 p2,1 ; l1 (1 − ξN2ξ −1/2 )bNξ −1/2,1 ∆Nξ −1,1 ξNξ l1 = (1 − ξN2ξ +1/2 ) c¯Nξ +1/2,1 − bNξ +1/2,1 + 2 ∆Nξ −1,1 ξNξ − ∆Nξ −1,1 ξNξ −1 dNξ ,2 ws + mNξ − + LNξ ,1 pNξ ,1 ;
Nξ ,1
eNξ ,1
∆2,1 (ξ3 ) − ∆2,1 (ξ2 )
hz 1
2
hz 1
2
2
lNz +1 (1 − ξ52/2 )b5/2,Nz +1 ∆2,Nz +1 (ξ2 )
lN +1 + z (1 − ξ32/2 ) c˜3/2,Nz +1 + b3/2,Nz +1 ∆2,Nz +1 (ξ3 ) − ∆ 2 ) (ξ 2 2,Nz +1 d2,Nz +1 ws + m2 + + L2,Nz +1 p2,Nz +1 ;
2,N +1
e2,Nzz +1 =
hzNz
Nξ ,Nz +1 eNξ ,Nz +1
+
−
ws
2
(1 − ξN2ξ +1/2 ) c¯Nξ +1/2,Nz +1 − bNξ +1/2,Nz +1 + dNξ ,Nz +1 ws + + mNξ + LNξ ,Nz +1 pNξ ,Nz +1 ;
=
lNz +1 (1 − ξN2ξ −1/2 )bNξ −1/2,Nz +1 ∆Nξ −1,Nz +1 ξNξ
lNz +1
2
hzNz
∆Nξ −1,Nz +1 ξNξ − ∆Nξ −1,Nz +1 ξNξ −1
2
– for i = 3, 4, . . . , Nξ − 1: l1 (1 − ξi2+1/2 )bi+1/2,1 ∆i,1 (ξi )
i,1
ei,1 =
∆i,1 (ξi+1 ) − ∆i,1 (ξi )
i,N +1
ei,Nzz +1 =
+
l1 (1 − ξi2−1/2 )bi−1/2,1 ∆i−1,1 (ξi )
∆i−1,1 (ξi ) − ∆i−1,1 (ξi−1 )
lNz +1 (1 − ξi2+1/2 )bi+1/2,Nz +1 ∆i,Nz +1 (ξi )
∆i,Nz +1 (ξi+1 ) −∆i,Nz +1 (ξi ) ws di,Nz +1 + + Li,Nz +1 pi,Nz +1 ; + mi hzNz
+
+ mi
di,2 hz1
−
ws
2
+ Li,1 pi,1 ;
lNz +1 (1 − ξi2−1/2 )bi−1/2,Nz +1 ∆i−1,Nz +1 (ξi )
∆i−1,Nz +1 (ξi ) − ∆i−1,Nz +1 (ξi−1 )
2
– for j = 2, 3, . . . , Nz : 2,j e2,j
=
lj (1 − ξ52/2 )b5/2,j ∆2,j (ξ2 )
∆2,j (ξ3 ) − ∆2,j (ξ2 )
=
lj 2
di,j+1
+
hz j
di,j
hzj−1
2 3/2
2
∆ i,j (ξi+1 ) − ∆i,j (ξ i )
+ mi Nξ ,j eNξ ,j
+ (1 − ξ
lj (1 − ξi2+1/2 )bi+1/2,j ∆i,j (ξi )
i ,j
ei,j =
lj
+
) c˜3/2,j + b3/2,j + m2
+
(1 − ξ
+ Li,j pi,j ,
i = 3, 4, . . . , Nξ − 1;
hz j
+
d N ξ ,j
hzj−1
lj (1 − ξN2ξ −1/2 )bNξ −1/2,j ∆Nξ −1,j ξNξ
∆Nξ −1,j ξNξ − ∆Nξ −1,j ξNξ −1
+ LNξ ,j pNξ ,j ;
i ,j
lj (1 − ξi2+1/2 )bi+1/2,j ∆i,j (ξi+1 )
i ,j
lj (1 − ξi2−1/2 )bi−1/2,j ∆i−1,j (ξi−1 )
, i = 2, 3, . . . , Nξ − 1; ∆i,j (ξi+1 ) − ∆i,j (ξi ) N ξ ,j eNξ +1,j = 0.5lj (1 − ξN2ξ +1/2 ) c¯Nξ +1/2,j + bNξ +1/2,j ; 2 ,j e1,j = 0.5lj (1 − ξ32/2 ) c˜3/2,j − b3/2,j ; ∆i−1,j (ξi ) − ∆i−1,j (ξi−1 )
+ L2,j p2,j ;
hzj−1
∆i−1,j (ξi ) − ∆i−1,j (ξi−1 )
) c¯Nξ +1/2,j − bNξ +1/2,j +
dNξ ,j+1
lj (1 − ξi2−1/2 )bi−1/2,j ∆i−1,j (ξi )
– for j = 1, 2, . . . , Nz + 1:
ei−1,j =
hzj
d2,j
2 Nξ +1/2
+ mNξ
ei+1,j =
d2,j+1
,
i = 3, 4, . . . , Nξ ;
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
9
– for i = 2, 3, . . . , Nξ : i,j ei,j−1
i,j ei,j+1
= mi
di,j hzj−1
= mi
−
ws
+
ws
di,j+1 hzj
,
2
j = 2, 3, . . . , Nz + 1;
,
2
j = 1, 2, . . . , Nz .
4. Non-negativity preserving of the scheme Following Theorem 2, in this section we will discuss the discrete analog of the NP. Namely, we will show that the system matrix of the full discretization is an M-matrix so that the discrete maximum principle is satisfied by the discretization and the discrete solution is non-negative. To discretize the system of the ODEs (22) we introduce the time mesh
w τ = {0 = t0 < t1 < · · · < tM = T } and the notations τk = tk+1 − tk , τ = max0≤k≤M −1 τk . Then, we apply the two-level implicit-stepping θ -method 0 < θ ≤ 1 to the system of the ODEs (22) and obtain a five-diagonal linear system of equations for the unknowns on the new, (k + 1)-th 1 time level Ci,j ≡ Cik,+ j , k = 0, 1, . . . , M − 1: i,j
i ,j
i,j
i ,j
i,j
e˜ i,j Cˆ i,j − e˜ i+1,j Cˆ i+1,j − e˜ i−1,j Cˆ i−1,j − e˜ i,j+1 Cˆ i,j+1 − e˜ i,j−1 Cˆ i,j−1 = Fi,j , i,1
i ,1
i,1
i,1
i = 2, 3, . . . , Nξ , j = 2, 3, . . . , Nz ;
e˜ i,1 Cˆ i,1 − e˜ i+1,1 Cˆ i+1,1 − e˜ i−1,1 Cˆ i−1,1 − e˜ i,2 Cˆ i,2 = Fi,1 , i ,N + 1
i,N +1
i ,N + 1
i,N +1
e˜ i,Nzz +1 Cˆ i,Nz +1 − e˜ i+1z,Nz +1 Cˆ i+1,Nz +1 − e˜ i−1z,Nz +1 Cˆ i−1,Nz +1 − e˜ i,Nzz
Cˆ i,Nz = Fi,Nz +1 ,
i = 2, 3, . . . , Nξ ;
Cˆ i,j = 0, i = 1, Nξ + 1; j = 1, 2, . . . , Nz + 1.
(23)
The coefficients and the right hand side in the system (23) for i = 2, 3, . . . , Nξ are as follows: i,j
e˜ i,j =
Li,j
+ θ eii,,jj ,
τk
i,j
i,j
i ,j
i ,j
e˜ i+1,j = θ ei+1,j ,
i ,j
i,j
e˜ i−1,j = θ ei−1,j , i,j
j = 1, 2, . . . , Nz + 1;
i,j
e˜ i,j+1 = θ ei,j+1 , j = 1, 2, . . . , Nz ; e˜ i,j−1 = θ ei,j−1 , j = 2, 3, . . . , Nz + 1; Li,1 − (1 − θ )eii,,11 Ci,1 + (1 − θ ) eii,+11,1 Ci+1,1 + eii,−11,1 Ci−1,1 + eii,,12 Ci,2 + Li,1 fR,i,1 , Fi,1 = τ k
Fi,j =
Li,j
τk
i,j ei,j
− (1 − θ )
i,j
i ,j
i,j
i ,j
Ci,j + (1 − θ ) ei+1,j Ci+1,j + ei−1,j Ci−1,j + ei,j+1 Ci,j+1 + ei,j−1 Ci,j−1 + Li,j fR,i,j ,
j = 2, 3, . . . , Nz ;
Fi,Nz +1 =
Li,Nz +1
i,Nz +1 i,Nz +1 1 i,Nz +1 − (1 − θ )eii,,NNzz + +1 Ci,Nz +1 + (1 − θ ) ei+1,Nz +1 Ci+1,Nz +1 + ei−1,Nz +1 Ci−1,Nz +1 + ei,Nz Ci,Nz
τk + Li,Nz +1 fR,i,Nz +1 .
Since the unknown function C (ξ , z , t ) is the concentration of a substance, the important property of the discretization is its monotonicity. This provides the positivity (non-negativity) of the numerical solution with respect to time if the initial concentration is non-negative, Theorem 2. Let us denote by A the matrix of the system (23). Theorem 3. For τ sufficiently small, matrix A is an M-matrix. Proof. As is well known, to prove that matrix A is an M-matrix, it is sufficient to show that its diagonal elements are positive, the off-diagonal elements are non-positive and matrix A is strictly diagonally dominant. i,j From the expressions for e˜ i,j for i = 2, 3, . . . , Nξ , j = 1, 2, . . . , Nz + 1 it is easy to see that the diagonal elements are always positive since τk is small. For the off-diagonal elements i = 3, 4, . . . , Nξ − 1, j = 2, 3, . . . , Nz we have i,j
e˜ i+1,j = θ lj (1 − ξi2+1/2 )ci+1/2,j
αi,j 1 − r¯i
αi,j 2
> 0,
10
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
Fig. 1. Exact solution at t = T .
where 0 < r¯i =
1+ξi 1−ξi+1 1+ξi+1 1−ξi i,j
sign of αi,j . Analogically e˜ i−1,j
αi,j
< 1. That is because θ lj (1 − ξi2+1/2 )ci+1/2,j > 0 and the sign of 1 − r¯i 2 is the same as the d d > 0. The coefficient e˜ ii,,jj+1 = θ mi ih,jz+1 + w2s > 0 and e˜ ii,,jj−1 = θ mi hz i,j − w2s > 0 when 2 ,j
2,j
2 ,j
N ξ ,j
N ξ ,j
N ξ ,j
2 ,N + 1
2,N +1
j
2,1
2,1
Nξ ,1
j−1 N ξ ,1
i,1
i,1
i ,1
hzj−1 is sufficiently small. For e˜ 3,j , e˜ 2,j+1 , e˜ 2,j−1 , e˜ Nξ −1,j , e˜ Nξ ,j+1 , e˜ Nξ ,j−1 , j = 2, . . . , Nz , e˜ 3,1 , e˜ 2,2 , e˜ Nξ −1,1 , e˜ Nξ ,2 , e˜ i+1,1 , e˜ i−1,1 , e˜ i,2 , i,N +1
i ,N + 1
i ,N + 1
e˜ i+1z,Nz +1 , e˜ i−1z,Nz +1 , e˜ i,Nzz 2 ,j e1,j ,
˜
N ξ ,j eNξ +1,j ,
˜
, i = 3, . . . , Nξ − 1, e˜ 3,Nzz +1 , e˜ 2,Nzz
Nξ ,Nz +1
Nξ ,Nz +1
, e˜ Nξ −1,Nz +1 , e˜ Nξ ,Nz
arguments are the same. The coefficients
j = 1, 2, . . . , Nz + 1, are multiplied by the right hand side of the boundary conditions (3b) and are transferred
to the right hand side of the equations. algebraic N +1 i,j Nξ +1 Nz +1 i,j ˜ p,l ≡ e˜ ii,,jj − p=ξ 1p̸=i Nl=z 1+l̸=1 j e˜ ip,,jl . Let us denote by Di,j = e˜ i,j − p=1 l=1l̸=j −e p̸=i For i = 2, 3, . . . , Nξ we get Di,j = Di,1 = Di,j =
Li,j
τk L i ,1
τk Li,j
τk
+ θ lj (1 − ξi2−1/2 )bi−1/2,j − lj (1 − ξi2+1/2 )bi+1/2,j + Li,j pi,j ,
j = 2, 3, . . . , Nz ;
+ θ l1 (1 − ξi2−1/2 )bi−1/2,1 − l1 (1 − ξi2+1/2 )bi+1/2,1 − mi ws + Li,1 pi,1 ; + θ lj (1 − ξi2−1/2 )bi−1/2,j − lj (1 − ξi2+1/2 )bi+1/2,j + mi ws + Li,j pi,j ,
j = Nz + 1.
For i = 1, Nξ + 1; j = 1, 2, . . . , Nz + 1 we have Di,j = 1. As is seen from the expressions for Di,j , when the maximum time step τ is small enough, matrix A is strictly diagonally dominant. In fact, when the meshes on the spatial variables are regular, one can show that τ ≤ C min(hξ , hz ), where C does not depend on the discretization parameters. Thus, matrix A is an M-matrix. From the expressions for Fi,j it follows that Fi,j ≥ 0 if the maximum time step τ is sufficiently small. As matrix A is an M-matrix, we can conclude that the numerical solution is non-negative, if the initial data are non-negative. 5. Numerical experiments Numerical experiments presented in this section illustrate the properties of the constructed scheme. We perform numerical experiments of two kinds. (1) For the first kind we have used the following values of the coefficients in the problem under consideration: a = 10, ws = 0.1, k = 0.38, Π = 0.09, kmd = 10−2 , z 0 = 0.2, zp = 0.2, T = 1. In order to investigate numerically the convergence and the accuracy constructed schemes, we approximately solve the model problem with known of the analytical solution Can = 0.008 1 − ξ 2 z −3 t. We choose this function because its feature (Fig. 2) is similar to that of the exact solution to the problem under consideration (Fig. 1), but the analytical solution is differentiable. In the computations below τk = const = 0.001 is used, because all of the experiments with different meshes in space show that the further decreasing of τk changes only the digits after the seventh one of the mesh C -norm of the error. The rate of convergence (RC) is calculated using the double mesh principle RC = log2 (ERNξ ,Nz /ER2Nξ ,2Nz ),
Nξ ,Nz
ERNξ ,Nz = ∥Can
N ,N z
ξ − Cnum ∥.
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
11
Fig. 2. Analytical solution at t = T . Table 1 C -norm and L2 -norm of the error, θ = 0.5. hξ = hz =
C -norm of the error
RC
L2 -norm of the error
RC
0.1 0.05 0.025 0.0125 0.00625
3.974 E−1 1.203 E−1 3.128 E−2 7.524 E−3 1.674 E−3
– 1.72 1.94 2.06 1.80
8.249 E−2 1.932 E−2 4.213 E−3 9.238 E−4 2.283 E−4
– 1.71 1.84 1.83 1.60
Table 2 C -norm and L2 -norm of the error, θ = 1.
Nξ ,Nz
hξ = hz =
C -norm of the error
RC
L2 -norm of the error
RC
0.1 0.05 0.025 0.0125 0.00625
3.977 E−1 1.208 E−1 3.175 E−2 8.013 E−3 2.166 E−3
– 1.20 1.49 1.57 1.43
8.256 E−2 1.942 E−2 4.301 E−3 1.002 E−3 2.904 E−4
– 1.68 1.81 1.72 1.29
Nξ ,Nz
Here Can and Cnum are respectively the analytical solution and the numerical solution computed at the mesh with Nξ and Nz subintervals, ∥.∥ is the mesh C -norm or L2 -norm of the error, where N ,Nz
∥Canξ
N ,Nz
∥Canξ
N ,Nz
ξ − Cnum ∥C =
N ξ ,N z
Can max 1 ≤ i ≤ Nξ + 1 1 ≤ j ≤ Nz + 1
i ,j
N ,N ξ z − Cnum
, i,j
Nξ +1 N +1 N ,N 2 z Nξ ,Nz Nξ ,Nz ξ z − Cnum ∥L2 = mi lj Can . − Cnum i =1
j =1
i,j
i,j
Results from computational experiments, concerning the error and the rate of convergence with respect to space variables at fixed value of t = T , are presented in Table 1 for θ = 0.5 and in Table 2 for θ = 1 respectively. (2) For the second kind of numerical experiments we solve the real problem. We have used the same values of the coefficients a, ws , k, Π , kmd , T , but take z 0 = 0.0002, zp = 0.1. The exact solution at these data is shown in Fig. 3. In this case we use the Runge method for the practical estimation of the rate of convergence s of the scheme with respect to the space variables at a fixed value of t = T . In the case when the exact solution is not known the formula for s is
Ch − Ch/2 s = ln C −C h/2
h/4
ln 2,
where Ch , Ch/2 and Ch/4 are computed values of the numerical solution at grid nodes with step h, h/2 and h/4 respectively. We use three inserted grids with 40 × 20, 80 × 40 and 160 × 80 subintervals respectively by ξ and z in the domain (ξ , z ) ∈ [−1, 1] × [z 0 , 1] and τ = τk = 0.0001. A part of the results from the calculations for the rate of convergence is presented in Table 3. It can be seen from Table 3 that the rate of convergence is very high when the nodes are not near to the point (0, zp ) where the Dirac-delta functions are concentrated.
12
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
Fig. 3. Exact solution at t = T . Table 3 Runge method for the rate of convergence, θ = 0.5. z
0.0002 0.05 0.10 0.15 0.20 0.25 0.30
ξ −0.20
−0.15
−0.10
−0.05
0
7.12 6.91 0.48 1.23 3.44 5.43 7.38
6.91 4.73 0.16 1.37 3.63 5.56 7.34
6.71 3.16 −0.16 1.59 3.85 5.70 7.39
6.50 2.67 −0.46 1.89 4.05 5.82 7.43
6.30 2.92 −0.74 2.23 4.16 5.86 7.44
0.05
0.10
0.15
0.20
6.08 2.42 −0.46 1.89 4.05 5.82 7.41
5.86 2.02 −0.16 1.59 3.84 5.70 7.33
5.64 1.72 0.16 1.35 3.62 5.53 7.22
5.42 1.48 0.51 1.16 3.39 5.34 7.07
6. Conclusions Here we deal with a domain transformation method, often used in computational mechanics. We study the transformed degenerate parabolic problem modeling settling particles dispersion for non-negativity. We construct and discuss a fitted finite volume difference scheme for the problem transformed on the finite domain Ω = (−1, 1) × (z 0 , 1). We may use any interval (−l, l), l > 0 instead of (−1, 1). This leads to insignificant computational cost of the solution. We have shown that the numerical scheme results in a monotone numerical approximation, i.e., it conserves the NP of the continuous solution. Experiments demonstrate the efficiency of our scheme. Besides the full time-dependent model presented here, it is interesting to consider the analogous steady problem. We intend to do this in the future. Acknowledgments The authors are grateful to the referees for their comment and suggestions which helped improve the quality of the manuscript. The first author is supported by the Sofia University Foundation under Grant No. 106/2013. The second author is supported by the Bulgarian Fund of Sciences under Grant No. DID 02/37/09. References [1] B.S. Mazumder, S. Bandyopadhyay, On solute dispersion from an elevated line source in an open-channel flow, J. Eng. Math. 40 (2001) 197–209. [2] K.K. Mondal, B.S. Mazumder, On dispersion of settling particles from an elevated source in an open-channel flow, J. Comput. Appl. Math. 193 (2006) 22–37. [3] K.K. Mondal, B.S. Mazumder, Dispersion of fine settling particles from an elevated source in an oscillatory turbulent flow, Eur. J. Mech. B Fluids 27 (6) (2008) 707–725. [4] N. Madden, K.K. Mondal, Improved mathematical and numerical modelling of dispersion of a solute from a continuous source, in: Lect. Notes Comp. Sci. Enj., vol. 81, Springer, 2011, pp. 177–185. [5] I. Nezu, H. Nakagawa, Turbulence in open-chanel flows, in: IAHR Monograph, A.A. Balkema Publishers, Rotterdam, 1993, p. 281. [6] P.J. Sullivan, H. Yip, Near-field contaminant dispersion from an elevated line-source, Z. Angew. Math. Phys. 38 (1987) 409–423. [7] K.W. Morton, Numerical Solution of Convection–Diffusion Problems, Chapman and Hall, London, 1996. [8] T. Chernogorova, R. Valkov, Finite volume difference scheme for a degenerate parabolic equation in the zero-coupon bond pricing, Math. Comput. Modelling 54 (2011) 2659–2671. [9] J.J.H. Miller, S. Wang, An exponentially fitted element volume method for the numerical solution of 2D unsteady incompressible flow problems, J. Comput. Phys. 115 (1994) 56–64. [10] S. Wang, A novel finite volume method for Black–Scholes equation governing option pricing, IMA J. Numer. Anal. 24 (2004) 699–720.
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
13
[11] O.A. Oleinik, E.V. Radkevich, Second Order Equations with Nonnegative Characteristic Form, Plenum, Press, New York, 1973, MR0457908. [12] T. Apel, O. Benedix, D. Sirch, B. Vexler, A priori mesh grading for an elliptic problem with Dirac right hand side, SIAM J. Numer. Anal. 49 (3) (2011) 992–1005. [13] X. Yang, X. Zhang, Z. Li, Guo-Wei He, A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations, J. Comput. Phys. 228 (2009) 7821–7836. [14] A. Kufner, Weighted Sobolev Spaces, John Wiley and Sons, New York, 1985. [15] J. Haslinger, M. Miettinen, P.D. Panagiotopoulos, Finite Element Method for Hemivariational Inequalities, Kluwer, Dordrecht, The Netherlands, 1999. [16] I. Farago, I. Horvath, Discrete maximum principle and adequate discretizations of linear parabolic problems, SIAM Sci. Comput. 28 (2006) 2313–2336. [17] D. Gilbert, N. Trudinger, Elliptic Partial Differential Equations of Second Order, 32nd ed., Springer, 2011. [18] O.A. Ladyzenskaja, V.A. Solonnikov, N.N. Ural’tseva, Linear and Quasilinear Equations of Parabolic Type, in: Amer. Math. Soc. Transl. Monographs, vol. 23, 1968. [19] A.A. Samarskii, Theory of Finite Difference Schemes, Marcel Decker, New York, 2003.
)
–
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
A finite volume difference scheme for a model of settling particle dispersion from an elevated source in an open-channel flow T. Chernogorova a,∗ , L. Vulkov b a
FMI, University of Sofia, 5, J. Bourchier Blvd., 1164 Sofia, Bulgaria
b
FNSE, University of Rousse, 8, Studentska str., 7017 Rousse, Bulgaria
article
info
Article history: Available online xxxx Keywords: Dispersion Settling velocity Degenerate parabolic equation Finite volume difference scheme Conservation of non-negativity
abstract A finite volume fitted difference scheme is constructed to solve the unsteady convectivediffusion equation transformed on a finite domain, modeling longitudinal dispersion of suspended particles with settling velocity in a turbulent shear flow over a rough-bed surface. First we discuss the well-posedness of the differential problem and the non-negativity of its solution. Then, to overcome the degeneracy at the part of the boundary, starting from the divergent form of the equation we perform a local fitted space discretization. This approximation is determined by a set of two-point boundary value problems. Non-negativity of the numerical concentration of suspended fine particles is proved. Some results from computational experiments are presented to illustrate the properties of the constructed scheme. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction We discuss a model due to Mondal and Mazumder [1–3], see also [4], for a dispersion of fine particles in an oscillatory turbulent flow. The model is based on the time-dependent advection–diffusion equation
∂C ∂C ∂ 2C ∂ ∂C ∂C + u(z ) − ws = kx (z ) 2 + kz (z ) + f (x, z ) ∂t ∂x ∂z ∂x ∂x ∂z on (−∞, ∞) × (z 0 , 1) and subjected to boundary and initial conditions ∂C C (±∞, z , t ) = 0, kz (z ) + ws C = 0, C (x, z , 0) = 0, ∂z z = z 0 ,1
(1)
(2)
where z 0 is the boundary roughness height, ws is the setting velocity and f (x, z ) = δ(x)δ(z − zp ), i.e., a point source at (x, z ) = (0, zp ). Here the problem has been expressed in terms of dimensionless variables and it is used for fully developed homogeneous turbulent flow [4,1,2]. In our present study, the velocity distribution is 1 z u(z ) = ln 0 + W (z ), k z where k > 0 is the von-Karman constant, z 0 is the equivalent bed roughness, W (z ) = 2kΠ sin2 π2 z , and Π is the wakestrength parameter.
∗
Corresponding author. Tel.: +359 897846436. E-mail addresses: [email protected] (T. Chernogorova), [email protected] (L. Vulkov).
http://dx.doi.org/10.1016/j.camwa.2014.03.012 0898-1221/© 2014 Elsevier Ltd. All rights reserved.
2
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
When the Reynolds decomposition is performed on the laminar advection–diffusion equation by decomposing the velocity (u = u + u′ ) and concentration (C = C + C ′ ) into the sum of their mean and fluctuating parts, the term (∂ u′ C ′ /∂ xj ) represents the transport of concentration C due to turbulent fluctuations. Because the molecular diffusivity (kmd ) is usually a very small quantity, u′ C ′ ≫ kmd ∂ C /∂ xj , the effects for molecular diffusion are often negligible compared to the effects of turbulence, and so omitted from the model. However, in this model the molecular diffusion is an important mechanism for mixing of the concentration to the flow at the smallest scales, and so the authors of [2,3] take the diffusivity profile as kx (x) = kmd + ked (z ), where ked (z ) is the eddy-diffusivity as proposed in [5,6]: ked (z ) = k(1 − z )
1 z
−1 + Π π sin(π z ) .
Since our problem is posed in a fully developed homogeneous turbulent flow, we take kx (z ) = kz (z ) and neglect the crossstream diffusion terms. The studies of dispersion phenomena of settling particles released from an elevated continuous line-source in a turbulent open-channel flow are important to understand the process of contamination in the environment [1–3,5]. Mazumder and Bandyopadhyay [1] presented a numerical solution to the convection–diffusion equation to study the stream-wise dispersion of contaminant due to turbulent shear flow over a gravel-bed surface, employing a modified ‘‘Logwake law’’ and corresponding diffusivity due to Nezu and Nakagawa [5]. Their method provides a more general numerical scheme for the unsteady convection–diffusion equation than the theoretical study of Sullivan and Yip [6] based on the small time asymptotic solution. A combined scheme of central difference for the diffusion terms and a four-point upwind scheme for convective terms are used to solve the convection–diffusion equation for the steady problem. The alternating direction implicit (ADI) method is adopted to generate the numerical solution of a two-dimensional unsteady convection–diffusion equation and is developed in [2,3]. The authors use this scheme for numerical simulations to examine the stream-wise dispersion of suspended fine particles with settling velocities in oscillatory turbulent shear flow when these particles are being released from an elevated source. They, in order to incorporate the zero Dirichlet boundary condition of (2), transform the unbounded region to a bounded one. We also will work with transformed problem, namely (3a)–(3d), in the next section. However, as we shall see, the features apparent in the model problem (3a)–(3d) present difficulties which require new numerical approaches. These difficulties include the degeneracy of the elliptic part of the parabolic equation (3a), noncoercivity of the bilinear elliptic form, corner points of the domain, etc. We provide a theoretical analysis in appropriate weighted Sobolev spaces of the continuous problem on bounded domain. Here we propose a new approach for the construction of an efficient difference scheme. This is a fitted technique based on the old idea of Allen and Southwell well described in [7] and that was extended to different problems by several authors [8–10]. In a previous work of the first author [8] numerical experiments show higher accuracy of the present construction in comparison with the other known difference schemes near the degeneracy. The rest of the paper is organized as follows. In Section 2 we present the differential problem on bounded domain and discuss our basic assumptions and some properties of the solution. The space discretization is developed in Section 3. Section 4 is devoted to the system matrix after the time discretization. We show that the maximum principle is satisfied and thus the discrete solution is non-negative in correspondence to the continuous solution. In Section 5 we present some results from computational experiments. In Section 6 we formulate some conclusions. 2. The differential problem on bounded domain To compute an accurate numerical solution to the model given above, several issues need to be addressed. First, the problem is posed on a semi-infinite domain. Numerical solutions are possible only on finite domains, so the domain must either be truncated or transformed. Here we take the approach, discussed in [4,1,2], since it has proved successful in related studies. In [2,3] a tanh-transformation from domain (x, z ) ∈ (−∞, ∞) × [z 0 , 1] to (ξ , z ) ∈ (−1, 1) × [z 0 , 1] is made using x=
1 2a
log
1+ξ 1−ξ
,
where a is a ‘‘stretching’’ factor, usually chosen based on the computational experience. Following the transformation the model becomes
∂C ∂C ∂ 2C ∂C + u(z )a(1 − ξ 2 ) − a2 kξ (z )(1 − ξ 2 ) (1 − ξ 2 ) 2 − 2ξ ∂t ∂ξ ∂ξ ∂ξ ∂ ∂C ∂C = kz (z ) + ws + δ(ξ )δ(z − zp ), ∂z ∂z ∂z l± (C ) ≡ C (±1, z , t ) = 0,
(3a) (3b)
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
lz 0 ,1 (C ) ≡ kz (z )
∂C + ws C ∂z
)
–
3
= 0,
z = z 0 ,1
(3c)
C (ξ , z , 0) = C0 (ξ , z ) = 0.
(3d)
The parabolic equation (3a) belongs to the second order PDEs with the non-negative characteristic form. The main character of such kind of equations is the degeneracy. It can be easily seen that at ξ = ∓1 Eq. (3a) degenerates into two parabolic one-dimensional equations:
∂ ∂C = ∂t ∂z
kz (z )
∂C ∂z
+ ws
∂C , ∂z
(z , t ) ∈ (z 0 , 1) × (0, T ].
(4)
By the well-known Fichera’s theory (see [11]) for degenerate parabolic equations, we have that at the degenerate boundaries ξ = −1 and ξ = 1, the boundary conditions should not be given. Therefore, the parabolic problem (3a), (4), (3c) is well defined. But from physical motivation we have imposed the boundary conditions (3b) and the initial condition (3d). It is easy to check that (4) follows from (3b), (3d). Therefore the boundary conditions (3b) are the particular case of (4). As discussed in [12], the Poisson equation with the Dirac-delta function on the right-hand side does not have an H 1 (Ω )solution. For this and also for the reason of computations (see Section 5) we will assume that the right hand side is already regularized and at least the regularized function belongs to L2 (Ω ), for example, see [13]. Also, this will simplify the next theoretical considerations. We rewrite Eq. (3a) in divergent form:
∂C − LC + p(ξ , z )C = fR (ξ , z ) for (ξ , z , t ) ∈ QT = Ω × (0, T ], ∂t
(5)
where LC = ∇ · (A∇ C + bC ) ,
A=
a11 0
0 a22
,
Ω = (−1, 1) × (z 0 , 1),
a11 = a2 kξ (z )(1 − ξ 2 )2 ,
a22 = kz (z ),
b1 (ξ , z ) (1 − ξ 2 )a 2ξ akξ (z ) − u(z ) , = b2 ws
p(ξ , z ) = 2a2 kξ (z )(1 − 3ξ 2 ) + 2ξ au(z ), and fR (ξ , z ) is a regularization of the delta function δ(ξ )δ(z − zp ). Further, to handle the degeneracy in Eq. (5), we introduce the weighted inner product and corresponding norm on L2,w (Ω ) by
(u, v)w :=
Ω
(1 − ξ 2 )2 uv dξ dz ,
∥v∥0,w =
(v, v)w =
Ω
(1 − ξ 2 )2 v 2 dξ dz
1/2
.
The space of all weighted square-integrable function is defined as L2,w (Ω ) := {v : ∥v∥0,w < ∞}. Using a standard argument it is easy to show that the pair (L2,w (Ω ), (·, ·)w ) is a Hilbert space (cf., for example [14]). Using L2 (Ω ) and L2,w (Ω ) we define the following weighted Sobolev space 1 Hw (Ω ) :=
∂v ∂v ∈ L2,w (Ω ) and ∈ L2 ( Ω ) v ∈ L2 (Ω ), ∂ξ ∂z
1 with the corresponding inner product. It is also easy to prove that Hw (Ω ) is a Hilbert space with the norm
∥v∥1,w :=
Ω
v + (1 − ξ ) 2
2 2
∂v ∂ξ
2
+
∂v ∂z
2
1/2 dξ dz
.
1 (Ω ) define the bilinear form Next, for C , η ∈ Hw
A(C , η) ≡
Ω
a11
∂ C ∂η ∂ C ∂η ∂η ∂η + a22 + b1 C + b2 C + pC η dξ dz . ∂ξ ∂ξ ∂z ∂z ∂ξ ∂z
Using this weighted Sobolev space we define the following variational problem corresponding to (5) and (3b)–(3d): find 1 C ∈ Hw (Ω ) satisfying the initial condition (3d) such that for all η ∈ Hw1 (Ω )
Ω
∂C ηdξ dz + A(C , η) = ∂t
Ω
fR ηdξ dz ,
a.e. in (0, T ].
1 Theorem 1. There exists a unique solution C ∈ Hw (Ω ) to problem (5), (3b)–(3d).
4
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
Proof. By changing the variables C := eρ t C , ρ = const > 0, we can replace the coefficient p in (5) by p + ρ (let us note that we will use C for the new function as well). Therefore, we will assume that
∂C − LC + (p + ρ)C = e−ρ t fR (ξ , z ), ∂t ∂ C ∂η ∂ C ∂η ∂η ∂η + a22 + b1 C + b2 C + (p + ρ)C η dξ dz . Aρ (C , η; t ) = a11 ∂ξ ∂ξ ∂z ∂z ∂ξ ∂z Ω For any ε > 0, there exists a constant K (ε) > 0 such that b1 C
2 2 ∂C ∂C ∂C ∂C ∂C ∂C + b2 C ≥ − b1 C + b2 C ≥ −ε(1 − ξ 2 )2 −ε − K (ε)C 2 . ∂ξ ∂z ∂ξ ∂z ∂ξ ∂z
Therefore, for sufficiently small ε and large ρ
2 2 ∂C ∂C 2 2 2 Aρ (C , C ; t ) ≥ (a kξ (z ) − ε)(1 − ξ ) + (kz (z ) − ε) ∂ξ ∂z Ω + (p + ρ − K (ε))C 2 dξ dz ≥ K (ε, ρ)∥C ∥2H 1 (Ω ) , w
where K (ε, ρ) > 0 is a constant. Now, in a standard way, using the Cauchy–Schwarz inequality one can prove that the bilinear form Aρ (C , η; t ) is continuous. Finally by Lemma 1 and Theorem 1.33 of [15], we have that the problem (5), (3b)–(3d) has a unique solution and thus we complete the proof. Considering the settling particles dispersion, the concentration C (x, z , t ) of the suspended particles cannot be negative if the concentration was non-negative in the initial state and on the boundary of the channel and in addition, the unity concentration is maintained at the steady line-source (x = 0, z = zp ) at all times. This property is called non-negativity preservation (NP). Conservation of this property is very important at the adequate discretizations of PDEs [16]. Further we analyze the NP for problem (3a)–(3d). 1 For a function v in Hw (Ω ) we denote the positive and negative parts of v by v + and v − . That is v = v + + v − ; v + ≥ 0 − and v ≤ 0. By [17, Lemma 7.6, p. 150], we have that Dv + =
Dv, 0,
if v > 0 , if v ≤ 0
Dv − =
Dv, 0,
if v < 0 , if v ≥ 0
∂v or Dv = Dz v = ∂v . Also v + v − = Dv + Dv − = Dv + v − = v + Dv − = 0, (ξ , z ) ∈ Ω . where Dv = Dξ v = ∂ξ ∂z We consider the following inequalities
∂C − LC + p(ξ , z )C ≥ 0 in QT , ∂t l± (C ) ≥ 0
in (z 0 , 1) × (0, T ]
(6)
and lz 0 ,1 (C ) ≥ 0
in (−1, 1) × (0, T ].
(7)
Formally multiplication (6) by a non-negative smooth function η and integration by parts lead to the following definition. 1 ,0 Let V2 (QT ) (see [18]) be the Banach space consisting of all functions in Hw (QT ) = C ((0, T ); Hw1 (Ω )) having finite norm
∥v∥V2 (QT ) = sup ∥v(·, ·; t )∥L2 (Ω ) + ∥∇v∥L2,w (QT ) , 0 ≤t ≤T
where T
∥∇v∥L2,w (QT ) =
0
∥v∥2H 1 (Ω ) (t )dt
1/2
w
.
1,1 (QT ) = Following [18], a function v ∈ V2 (QT ) is said to satisfy weakly (6) and (7) if for any non-negative η ∈ Hw 1 H 1 ((0, T ); Hw (Ω ))
Ω
C (ξ , z , t )η(ξ , z , t )dξ dz −
Qt
C ηt dξ dzdτ +
t
A(C , η)dτ ≥ 0
(8)
0
for almost every t ∈ (0, T ]. Here Qt = Ω × (0, t ). It is easy to see that if C ∈ C 2,1 (Q T ) satisfies (6), (7) point-wise then (8) holds for C . Conversely, if C ∈ V2 (QT ) and satisfies (8) and C is sufficiently smooth then it also verifies (6), (7) point-wise. We will now prove the following assertion. Theorem 2. Let C (ξ , z , t ) be in V2 (QT ) and satisfy (6), (7). Then C (ξ , z , t ) ≥ 0 in QT .
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
5
Proof. In fact, we will work again with Aρ (C , η). Using the Steklov average and taking the limit (see [18], Chapter 3) we can formally take η = −C − ≥ 0 in (8) to obtain 1
2
Ω
(C − (ξ , z , t ))2 dξ dz
+ Qt
∂C ∂C− ∂C ∂C− ∂C− ∂C− − a11 + a22 + b1 C + b2 C + (p + ρ)CC dξ dzdτ ≤ 0. ∂ξ ∂ξ ∂z ∂z ∂ξ ∂z
(9)
We have DC .DC − = (DC − )2 , CC − = (C − )2 . Let us introduce Kz = minz 0 ≤z ≤1 kz (z ), A11 = a2 Kz . Next, using the estimate in the proof of Theorem 1 for the expression b1 C with ε =
∂C− ∂C− + b2 C ∂ξ ∂z 1 2
min(Kz , A11 ), we obtain from (9)
− 2 − 2 ∂C ∂C 2 2 (C (ξ , z , t )) dzdξ + min(Kz , A11 ) + (1 − ξ ) dξ dzdτ 2 Ω ∂ξ ∂z Qt ≤ [K (ε) − (p + ρ)](C − )2 dξ dzdτ . 1
−
2
Qt
Now, we choose ρ sufficiently large such that K (ε) − (p + ρ) ≤ 0. The above yields 1 2
(C (ξ , z , t )) dξ dz + min(Kz , A11 ) −
Ω
2
(1 − ξ )
2 2
Qt
for a. e. t ∈ (0, T ]. This implies that C (ξ , z , t ) ≥ 0 for a. e. t ∈ (0, T ].
1 2
Ω
∂C− ∂ξ
2
+
∂C ∂z
2
dξ dzdτ ≤ 0,
(C − (ξ , z , t ))2 dξ dz = 0 and therefore C − (ξ , z , t ) ≡ 0. We then conclude that
3. Space discretization In this section we propose a finite volume method with a fitting technique for the numerical solution of problem (5), (3b)–(3d). We have chosen this method due to the fact that it takes into account the degeneracy of the equation and for the problem investigated in [8] gives more accurate results near the boundaries of degeneration compared to the commonly used finite difference methods. We introduce the meshes
−1 = ξ1 < ξ2 < · · · < ξNξ +1 = 1,
z 0 = z1 < z2 < · · · < zNz +1 = 1
and put hξi = ξi+1 − ξi , i = 1, 2, . . . , Nξ , hzj = zj+1 − zj , j = 1, 2, . . . , Nz . We also introduce the auxiliary meshes
ξi−1/2 =
ξi−1 + ξi 2
,
i = 2, 3, . . . , Nξ + 1,
zj−1/2 =
zj−1 + zj 2
,
j = 2, 3, . . . , Nz + 1
and let ξ1/2 = ξ1 = −1, ξNξ +3/2 = ξNξ +1 = 1, z1/2 = z1 = z 0 , zNz +3/2 = zNz +1 = 1, mi = ξi+1/2 − ξi−1/2 , lj = zj+1/2 − zj−1/2 , Li,j = mi lj , i = 1, 2, . . . , Nξ + 1, j = 1, 2, . . . , Nz + 1. A. Internal grid nodes We integrate Eq. (5) on the cell ℜi,j = [ξi−1/2 , ξi+1/2 ] × [zj−1/2 , zj+1/2 ], i = 2, 3, . . . , Nξ ; j = 2, 3, . . . , Nz :
ℜi,j
∂C dξ dz = ∂t
∇ · (k(C ))dξ dz −
ℜi,j
pCdξ dz +
ℜi,j
fR dξ dz ,
(10)
ℜi,j
where k(C ) = A∇ C + bC . For all integrals with the exception of the second one in (10) we use the quadrature formula of the central rectangles, to get
∂ Ci,j Li,j = ∂t
∇ · (k(C ))dξ dz − Li,j pi,j Ci,j + Li,j fR,i,j ,
(11)
ℜi,j
where
∂ C ∂ Ci,j = , ∂ t (ξi ,zj ,t ) ∂t
p|(ξi ,zj ,t ) = pi,j ,
C |(ξi ,zj ,t ) = Ci,j ,
fR |(ξi ,zj ,t ) = fR,i,j .
6
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
Using the Gauss–Ostrogradski formula we find for the second integral in (10)
∇ · (k(C ))dξ dz =
ℜi,j
k · n ds =
(ξi−1/2 ,zj+1/2 )
− (ξi−1/2 ,zj−1/2 ) (ξi+1/2 ,zj−1/2 )
∂C + b1 C dz ∂ξ (ξi+1/2 ,zj−1/2 ) (ξi+1/2 ,zj+1/2 ) ∂C ∂C a22 + b1 C dz + + b2 C dξ a11 ∂ξ ∂z (ξi−1/2 ,zj+1/2 ) ∂C a22 + b2 C dξ = I1 − I2 + I3 − I4 . ∂z
a11
∂ℜi,j
(ξi+1/2 ,zj+1/2 )
− (ξi−1/2 ,zj−1/2 )
(12)
Let us consider the first integral in (12). From the expressions for a11 and b1 we get a11
∂C + b1 C = (1 − ξ 2 )ρ(C ), ∂ξ
where the flux is
ρ(C ) = a2 kξ (z )(1 − ξ 2 )
∂C + a 2ξ akξ (z ) − u(z ) C . ∂ξ
Further we will follow the basic idea of [10]. To approximate the flux ρ(C ) we solve the boundary value problem
′
ci+1/2,j (1 − ξ 2 )v ′ + bi+1/2,j v ξ = 0,
(13a)
v(ξi , zj ) = Ci,j ,
v(ξi+1 , zj ) = Ci+1,j , (13b) where c = a kξ (z ), b = a 2ξ akξ (z ) − u(z ) , ci+1/2,j = c (ξi+1/2 , zj ), bi+1/2,j = b(ξi+1/2 , zj ) and Ci,j , Ci+1,j are the values of the approximate solution C at the nodes (ξi , zj ) and (ξi+1 , zj ). Let us note that the coefficient c does not depend on ξ and in this case we have ci+1/2,j = c (zj ). We introduce the notations 2
αi,j =
bi+1/2,j ci+1/2,j
,
∆i,j (ξi ) =
1 + ξi
α2i,j
1 − ξi
.
Then for i = 2, 3, . . . , Nξ − 1, j = 2, 3, . . . , Nz we obtain I1 ≈ lj (1 − ξ 2 )ρ(C ) (ξ
i+1/2, zj )
= lj (1 − ξi2+1/2 )bi+1/2,j
∆i,j (ξi+1 ) Ci+1,j − ∆i,j (ξi ) Ci,j . ∆i,j (ξi+1 ) − ∆i,j (ξi )
(14)
As noted in [10] for i = Nξ the above analysis does not apply to the approximation of the flux because (13a) is degenerate and this leads to the solution of (13a) which can never satisfy both conditions in (13b) unless CNξ ,j = CNξ +1,j . To solve this difficulty let us re-consider (13a)–(13b) with an extra degree of freedom, solving the problem c¯Nξ +1/2,j (1 − ξ )v ′ + bNξ +1/2,j v
v(ξNξ , zj ) = CNξ ,j ,
′
ξ
= M2 ,
v(ξNξ +1 , zj ) = CNξ +1,j = 0,
where c¯ = (1 + ξ )a2 kξ (z ) and find for the first integral in (12) the approximation I1 ≈ lj (1 − ξN2ξ +1/2 )
1 CNξ +1,j c¯Nξ +1/2,j + bNξ +1/2,j − CNξ ,j c¯Nξ +1/2,j − bNξ +1/2,j . 2
(15)
By analogy with the investigations above we approximate the second integral in (12) and find for j = 2, 3, . . . , Nz – for 3 ≤ i ≤ Nξ :
I2 ≈ lj (1 − ξ 2 )ρ(C ) (ξ
i−1/2, zj )
∆i−1,j (ξi ) Ci,j − ∆i−1,j (ξi−1 ) Ci−1,j ; ∆i−1,j (ξi ) − ∆i−1,j (ξi−1 )
(16)
1 c˜3/2,j + b3/2,j C2,j − c˜3/2,j − b3/2,j C1,j ,
(17)
= lj (1 − ξi2−1/2 )bi−1/2,j
– for i = 2: I2 ≈ lj (1 − ξ32/2 ) ρ(C )|(ξ3/2, zj ) = lj (1 − ξ32/2 ) where c˜ = (1 − ξ )a2 kξ (z ).
2
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
7
The problem is non-degenerate with respect to the variable z and for the third integral in (12), following [19], we apply a classical approximation: (ξi+1/2 ,zj+1/2 )
∂C ∂C + b2 C dξ ≈ kz (z ) + ws C mi ∂z ∂z (ξi−1/2 ,zj+1/2 ) (ξi ,zj+1/2 ) ws di,j+1 ws di,j+1 + − Ci,j+1 − Ci,j , = mi
I3 =
a22
hzj
2
hzj
2
(18)
where
di,j+1 =
zj+1
1 hzj
−1
1 kz (z )
zj
≈ kz (zj+1/2 ).
dz
Analogically for the fourth integral in (12) we get
I4 ≈ m i
di,j
+
hzj−1
ws 2
C i ,j −
di,j hzj−1
−
ws
Ci,j−1 .
2
(19)
Replacing in (12) the obtained approximations for the integrals from (14)–(19) and then replacing the derived expression in (11), we obtain a semi-discretization of Eq. (5) in the internal grid points for i = 2, 3, . . . , Nξ , j = 2, 3, . . . , Nz . B. Boundary grid nodes For the left vertical boundary ξ = −1, z ∈ [z 0 , 1] and for the right vertical one ξ = 1, z ∈ [z 0 , 1] we have respectively C1,j = 0,
CNξ +1,j = 0,
j = 1, 2, . . . , Nz + 1.
For the down horizontal boundary ξ ∈ (−1, 1), z = z 0 we proceed as for the internal grid nodes, but integrating Eq. (5) in the region ℜi,1 = [ξi−1/2 , ξi+1/2 ] × [z1/2 , z3/2 ], i = 2, 3, . . . , Nξ (i.e., in the semi-interval by z). As in (11) we get
∂ Ci,1 L i ,1 = ∂t
∇ · (k(C ))dξ dz − Li,1 pi,1 Ci,1 + Li,1 fR,i,1 .
(20)
ℜi,1
For the integral in (20) we have (ξi+1/2 ,z3/2 )
∂C + b1 C dz ∇ · (k(C ))dξ dz = k · n ds = a11 ∂ξ ℜi,1 ∂ℜi,1 (ξi+1/2 ,z1 ) (ξi−1/2 ,z3/2 ) (ξi+1/2 ,z3/2 ) ∂C ∂C − a11 + b1 C d ξ + a22 + b2 C d ξ ∂ξ ∂z (ξi−1/2 ,z1 ) (ξi−1/2 ,z3/2 ) (ξi+1/2 ,z1 ) ∂C − a22 + b2 C dξ = I1d − I2d + I3d − I4d . ∂z (ξi−1/2 ,z1 )
(21)
We find approximations for I1d , I2d and I3d by proceeding analogically as for the internal grid nodes. In fact these approximations can be obtained from (14) to (18) for j = 1. For I4d we have I4d = 0. For the upper horizontal boundary we accomplish analogical operations integrating Eq. (5) in the region ℜi,Nz +1 = [ξi−1/2 , ξi+1/2 ] × [zNz +1/2 , zNz +3/2 ], i = 2, 3, . . . , Nξ . Using all approximations for the internal and boundary grid nodes we obtain the following semi-discrete problem (note that Ci,j = Ci,j (t ) for i = 1, 2, . . . , Nξ + 1, j = 1, 2, . . . , Nz + 1):
∂ Ci,j i ,j i,j i ,j i,j i ,j Li,j = −ei,j Ci,j + ei+1,j Ci+1,j + ei−1,j Ci−1,j + ei,j+1 Ci,j+1 + ei,j−1 Ci,j−1 + Li,j fR,i,j , ∂t i = 2, 3, . . . , Nξ , j = 2, 3, . . . , Nz ; ∂ Ci,1 i ,1 i,1 i,1 i,1 Li,1 = −ei,1 Ci,1 + ei+1,1 Ci+1,1 + ei−1,1 Ci−1,1 + ei,2 Ci,2 + Li,1 fR,i,1 , ∂t ∂ Ci,Nz +1 i ,N + 1 i,N +1 Li,Nz +1 = −ei,Nzz +1 Ci,Nz +1 + ei+1z,Nz +1 Ci+1,Nz +1 ∂t + eii,−N1z,+N1z +1 Ci−1,Nz +1 + eii,,NNzz +1 Ci,Nz + Li,Nz +1 fR,i,Nz +1 , i = 2, 3, . . . , Nξ ; C1,j = 0, CNξ +1,j = 0, j = 1, 2, . . . , Nz + 1.
(22)
8
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
The expressions for the coefficients of the system (22) are given below: 2,1 e2,1
=
l1 (1 − ξ52/2 )b5/2,1 ∆2,1 (ξ2 )
l1
(1 − ξ
2 3/2
) c˜3/2,1 + b3/2,1 + m2
d2,2
+ L2,1 p2,1 ; l1 (1 − ξN2ξ −1/2 )bNξ −1/2,1 ∆Nξ −1,1 ξNξ l1 = (1 − ξN2ξ +1/2 ) c¯Nξ +1/2,1 − bNξ +1/2,1 + 2 ∆Nξ −1,1 ξNξ − ∆Nξ −1,1 ξNξ −1 dNξ ,2 ws + mNξ − + LNξ ,1 pNξ ,1 ;
Nξ ,1
eNξ ,1
∆2,1 (ξ3 ) − ∆2,1 (ξ2 )
hz 1
2
hz 1
2
2
lNz +1 (1 − ξ52/2 )b5/2,Nz +1 ∆2,Nz +1 (ξ2 )
lN +1 + z (1 − ξ32/2 ) c˜3/2,Nz +1 + b3/2,Nz +1 ∆2,Nz +1 (ξ3 ) − ∆ 2 ) (ξ 2 2,Nz +1 d2,Nz +1 ws + m2 + + L2,Nz +1 p2,Nz +1 ;
2,N +1
e2,Nzz +1 =
hzNz
Nξ ,Nz +1 eNξ ,Nz +1
+
−
ws
2
(1 − ξN2ξ +1/2 ) c¯Nξ +1/2,Nz +1 − bNξ +1/2,Nz +1 + dNξ ,Nz +1 ws + + mNξ + LNξ ,Nz +1 pNξ ,Nz +1 ;
=
lNz +1 (1 − ξN2ξ −1/2 )bNξ −1/2,Nz +1 ∆Nξ −1,Nz +1 ξNξ
lNz +1
2
hzNz
∆Nξ −1,Nz +1 ξNξ − ∆Nξ −1,Nz +1 ξNξ −1
2
– for i = 3, 4, . . . , Nξ − 1: l1 (1 − ξi2+1/2 )bi+1/2,1 ∆i,1 (ξi )
i,1
ei,1 =
∆i,1 (ξi+1 ) − ∆i,1 (ξi )
i,N +1
ei,Nzz +1 =
+
l1 (1 − ξi2−1/2 )bi−1/2,1 ∆i−1,1 (ξi )
∆i−1,1 (ξi ) − ∆i−1,1 (ξi−1 )
lNz +1 (1 − ξi2+1/2 )bi+1/2,Nz +1 ∆i,Nz +1 (ξi )
∆i,Nz +1 (ξi+1 ) −∆i,Nz +1 (ξi ) ws di,Nz +1 + + Li,Nz +1 pi,Nz +1 ; + mi hzNz
+
+ mi
di,2 hz1
−
ws
2
+ Li,1 pi,1 ;
lNz +1 (1 − ξi2−1/2 )bi−1/2,Nz +1 ∆i−1,Nz +1 (ξi )
∆i−1,Nz +1 (ξi ) − ∆i−1,Nz +1 (ξi−1 )
2
– for j = 2, 3, . . . , Nz : 2,j e2,j
=
lj (1 − ξ52/2 )b5/2,j ∆2,j (ξ2 )
∆2,j (ξ3 ) − ∆2,j (ξ2 )
=
lj 2
di,j+1
+
hz j
di,j
hzj−1
2 3/2
2
∆ i,j (ξi+1 ) − ∆i,j (ξ i )
+ mi Nξ ,j eNξ ,j
+ (1 − ξ
lj (1 − ξi2+1/2 )bi+1/2,j ∆i,j (ξi )
i ,j
ei,j =
lj
+
) c˜3/2,j + b3/2,j + m2
+
(1 − ξ
+ Li,j pi,j ,
i = 3, 4, . . . , Nξ − 1;
hz j
+
d N ξ ,j
hzj−1
lj (1 − ξN2ξ −1/2 )bNξ −1/2,j ∆Nξ −1,j ξNξ
∆Nξ −1,j ξNξ − ∆Nξ −1,j ξNξ −1
+ LNξ ,j pNξ ,j ;
i ,j
lj (1 − ξi2+1/2 )bi+1/2,j ∆i,j (ξi+1 )
i ,j
lj (1 − ξi2−1/2 )bi−1/2,j ∆i−1,j (ξi−1 )
, i = 2, 3, . . . , Nξ − 1; ∆i,j (ξi+1 ) − ∆i,j (ξi ) N ξ ,j eNξ +1,j = 0.5lj (1 − ξN2ξ +1/2 ) c¯Nξ +1/2,j + bNξ +1/2,j ; 2 ,j e1,j = 0.5lj (1 − ξ32/2 ) c˜3/2,j − b3/2,j ; ∆i−1,j (ξi ) − ∆i−1,j (ξi−1 )
+ L2,j p2,j ;
hzj−1
∆i−1,j (ξi ) − ∆i−1,j (ξi−1 )
) c¯Nξ +1/2,j − bNξ +1/2,j +
dNξ ,j+1
lj (1 − ξi2−1/2 )bi−1/2,j ∆i−1,j (ξi )
– for j = 1, 2, . . . , Nz + 1:
ei−1,j =
hzj
d2,j
2 Nξ +1/2
+ mNξ
ei+1,j =
d2,j+1
,
i = 3, 4, . . . , Nξ ;
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
9
– for i = 2, 3, . . . , Nξ : i,j ei,j−1
i,j ei,j+1
= mi
di,j hzj−1
= mi
−
ws
+
ws
di,j+1 hzj
,
2
j = 2, 3, . . . , Nz + 1;
,
2
j = 1, 2, . . . , Nz .
4. Non-negativity preserving of the scheme Following Theorem 2, in this section we will discuss the discrete analog of the NP. Namely, we will show that the system matrix of the full discretization is an M-matrix so that the discrete maximum principle is satisfied by the discretization and the discrete solution is non-negative. To discretize the system of the ODEs (22) we introduce the time mesh
w τ = {0 = t0 < t1 < · · · < tM = T } and the notations τk = tk+1 − tk , τ = max0≤k≤M −1 τk . Then, we apply the two-level implicit-stepping θ -method 0 < θ ≤ 1 to the system of the ODEs (22) and obtain a five-diagonal linear system of equations for the unknowns on the new, (k + 1)-th 1 time level Ci,j ≡ Cik,+ j , k = 0, 1, . . . , M − 1: i,j
i ,j
i,j
i ,j
i,j
e˜ i,j Cˆ i,j − e˜ i+1,j Cˆ i+1,j − e˜ i−1,j Cˆ i−1,j − e˜ i,j+1 Cˆ i,j+1 − e˜ i,j−1 Cˆ i,j−1 = Fi,j , i,1
i ,1
i,1
i,1
i = 2, 3, . . . , Nξ , j = 2, 3, . . . , Nz ;
e˜ i,1 Cˆ i,1 − e˜ i+1,1 Cˆ i+1,1 − e˜ i−1,1 Cˆ i−1,1 − e˜ i,2 Cˆ i,2 = Fi,1 , i ,N + 1
i,N +1
i ,N + 1
i,N +1
e˜ i,Nzz +1 Cˆ i,Nz +1 − e˜ i+1z,Nz +1 Cˆ i+1,Nz +1 − e˜ i−1z,Nz +1 Cˆ i−1,Nz +1 − e˜ i,Nzz
Cˆ i,Nz = Fi,Nz +1 ,
i = 2, 3, . . . , Nξ ;
Cˆ i,j = 0, i = 1, Nξ + 1; j = 1, 2, . . . , Nz + 1.
(23)
The coefficients and the right hand side in the system (23) for i = 2, 3, . . . , Nξ are as follows: i,j
e˜ i,j =
Li,j
+ θ eii,,jj ,
τk
i,j
i,j
i ,j
i ,j
e˜ i+1,j = θ ei+1,j ,
i ,j
i,j
e˜ i−1,j = θ ei−1,j , i,j
j = 1, 2, . . . , Nz + 1;
i,j
e˜ i,j+1 = θ ei,j+1 , j = 1, 2, . . . , Nz ; e˜ i,j−1 = θ ei,j−1 , j = 2, 3, . . . , Nz + 1; Li,1 − (1 − θ )eii,,11 Ci,1 + (1 − θ ) eii,+11,1 Ci+1,1 + eii,−11,1 Ci−1,1 + eii,,12 Ci,2 + Li,1 fR,i,1 , Fi,1 = τ k
Fi,j =
Li,j
τk
i,j ei,j
− (1 − θ )
i,j
i ,j
i,j
i ,j
Ci,j + (1 − θ ) ei+1,j Ci+1,j + ei−1,j Ci−1,j + ei,j+1 Ci,j+1 + ei,j−1 Ci,j−1 + Li,j fR,i,j ,
j = 2, 3, . . . , Nz ;
Fi,Nz +1 =
Li,Nz +1
i,Nz +1 i,Nz +1 1 i,Nz +1 − (1 − θ )eii,,NNzz + +1 Ci,Nz +1 + (1 − θ ) ei+1,Nz +1 Ci+1,Nz +1 + ei−1,Nz +1 Ci−1,Nz +1 + ei,Nz Ci,Nz
τk + Li,Nz +1 fR,i,Nz +1 .
Since the unknown function C (ξ , z , t ) is the concentration of a substance, the important property of the discretization is its monotonicity. This provides the positivity (non-negativity) of the numerical solution with respect to time if the initial concentration is non-negative, Theorem 2. Let us denote by A the matrix of the system (23). Theorem 3. For τ sufficiently small, matrix A is an M-matrix. Proof. As is well known, to prove that matrix A is an M-matrix, it is sufficient to show that its diagonal elements are positive, the off-diagonal elements are non-positive and matrix A is strictly diagonally dominant. i,j From the expressions for e˜ i,j for i = 2, 3, . . . , Nξ , j = 1, 2, . . . , Nz + 1 it is easy to see that the diagonal elements are always positive since τk is small. For the off-diagonal elements i = 3, 4, . . . , Nξ − 1, j = 2, 3, . . . , Nz we have i,j
e˜ i+1,j = θ lj (1 − ξi2+1/2 )ci+1/2,j
αi,j 1 − r¯i
αi,j 2
> 0,
10
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
Fig. 1. Exact solution at t = T .
where 0 < r¯i =
1+ξi 1−ξi+1 1+ξi+1 1−ξi i,j
sign of αi,j . Analogically e˜ i−1,j
αi,j
< 1. That is because θ lj (1 − ξi2+1/2 )ci+1/2,j > 0 and the sign of 1 − r¯i 2 is the same as the d d > 0. The coefficient e˜ ii,,jj+1 = θ mi ih,jz+1 + w2s > 0 and e˜ ii,,jj−1 = θ mi hz i,j − w2s > 0 when 2 ,j
2,j
2 ,j
N ξ ,j
N ξ ,j
N ξ ,j
2 ,N + 1
2,N +1
j
2,1
2,1
Nξ ,1
j−1 N ξ ,1
i,1
i,1
i ,1
hzj−1 is sufficiently small. For e˜ 3,j , e˜ 2,j+1 , e˜ 2,j−1 , e˜ Nξ −1,j , e˜ Nξ ,j+1 , e˜ Nξ ,j−1 , j = 2, . . . , Nz , e˜ 3,1 , e˜ 2,2 , e˜ Nξ −1,1 , e˜ Nξ ,2 , e˜ i+1,1 , e˜ i−1,1 , e˜ i,2 , i,N +1
i ,N + 1
i ,N + 1
e˜ i+1z,Nz +1 , e˜ i−1z,Nz +1 , e˜ i,Nzz 2 ,j e1,j ,
˜
N ξ ,j eNξ +1,j ,
˜
, i = 3, . . . , Nξ − 1, e˜ 3,Nzz +1 , e˜ 2,Nzz
Nξ ,Nz +1
Nξ ,Nz +1
, e˜ Nξ −1,Nz +1 , e˜ Nξ ,Nz
arguments are the same. The coefficients
j = 1, 2, . . . , Nz + 1, are multiplied by the right hand side of the boundary conditions (3b) and are transferred
to the right hand side of the equations. algebraic N +1 i,j Nξ +1 Nz +1 i,j ˜ p,l ≡ e˜ ii,,jj − p=ξ 1p̸=i Nl=z 1+l̸=1 j e˜ ip,,jl . Let us denote by Di,j = e˜ i,j − p=1 l=1l̸=j −e p̸=i For i = 2, 3, . . . , Nξ we get Di,j = Di,1 = Di,j =
Li,j
τk L i ,1
τk Li,j
τk
+ θ lj (1 − ξi2−1/2 )bi−1/2,j − lj (1 − ξi2+1/2 )bi+1/2,j + Li,j pi,j ,
j = 2, 3, . . . , Nz ;
+ θ l1 (1 − ξi2−1/2 )bi−1/2,1 − l1 (1 − ξi2+1/2 )bi+1/2,1 − mi ws + Li,1 pi,1 ; + θ lj (1 − ξi2−1/2 )bi−1/2,j − lj (1 − ξi2+1/2 )bi+1/2,j + mi ws + Li,j pi,j ,
j = Nz + 1.
For i = 1, Nξ + 1; j = 1, 2, . . . , Nz + 1 we have Di,j = 1. As is seen from the expressions for Di,j , when the maximum time step τ is small enough, matrix A is strictly diagonally dominant. In fact, when the meshes on the spatial variables are regular, one can show that τ ≤ C min(hξ , hz ), where C does not depend on the discretization parameters. Thus, matrix A is an M-matrix. From the expressions for Fi,j it follows that Fi,j ≥ 0 if the maximum time step τ is sufficiently small. As matrix A is an M-matrix, we can conclude that the numerical solution is non-negative, if the initial data are non-negative. 5. Numerical experiments Numerical experiments presented in this section illustrate the properties of the constructed scheme. We perform numerical experiments of two kinds. (1) For the first kind we have used the following values of the coefficients in the problem under consideration: a = 10, ws = 0.1, k = 0.38, Π = 0.09, kmd = 10−2 , z 0 = 0.2, zp = 0.2, T = 1. In order to investigate numerically the convergence and the accuracy constructed schemes, we approximately solve the model problem with known of the analytical solution Can = 0.008 1 − ξ 2 z −3 t. We choose this function because its feature (Fig. 2) is similar to that of the exact solution to the problem under consideration (Fig. 1), but the analytical solution is differentiable. In the computations below τk = const = 0.001 is used, because all of the experiments with different meshes in space show that the further decreasing of τk changes only the digits after the seventh one of the mesh C -norm of the error. The rate of convergence (RC) is calculated using the double mesh principle RC = log2 (ERNξ ,Nz /ER2Nξ ,2Nz ),
Nξ ,Nz
ERNξ ,Nz = ∥Can
N ,N z
ξ − Cnum ∥.
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
11
Fig. 2. Analytical solution at t = T . Table 1 C -norm and L2 -norm of the error, θ = 0.5. hξ = hz =
C -norm of the error
RC
L2 -norm of the error
RC
0.1 0.05 0.025 0.0125 0.00625
3.974 E−1 1.203 E−1 3.128 E−2 7.524 E−3 1.674 E−3
– 1.72 1.94 2.06 1.80
8.249 E−2 1.932 E−2 4.213 E−3 9.238 E−4 2.283 E−4
– 1.71 1.84 1.83 1.60
Table 2 C -norm and L2 -norm of the error, θ = 1.
Nξ ,Nz
hξ = hz =
C -norm of the error
RC
L2 -norm of the error
RC
0.1 0.05 0.025 0.0125 0.00625
3.977 E−1 1.208 E−1 3.175 E−2 8.013 E−3 2.166 E−3
– 1.20 1.49 1.57 1.43
8.256 E−2 1.942 E−2 4.301 E−3 1.002 E−3 2.904 E−4
– 1.68 1.81 1.72 1.29
Nξ ,Nz
Here Can and Cnum are respectively the analytical solution and the numerical solution computed at the mesh with Nξ and Nz subintervals, ∥.∥ is the mesh C -norm or L2 -norm of the error, where N ,Nz
∥Canξ
N ,Nz
∥Canξ
N ,Nz
ξ − Cnum ∥C =
N ξ ,N z
Can max 1 ≤ i ≤ Nξ + 1 1 ≤ j ≤ Nz + 1
i ,j
N ,N ξ z − Cnum
, i,j
Nξ +1 N +1 N ,N 2 z Nξ ,Nz Nξ ,Nz ξ z − Cnum ∥L2 = mi lj Can . − Cnum i =1
j =1
i,j
i,j
Results from computational experiments, concerning the error and the rate of convergence with respect to space variables at fixed value of t = T , are presented in Table 1 for θ = 0.5 and in Table 2 for θ = 1 respectively. (2) For the second kind of numerical experiments we solve the real problem. We have used the same values of the coefficients a, ws , k, Π , kmd , T , but take z 0 = 0.0002, zp = 0.1. The exact solution at these data is shown in Fig. 3. In this case we use the Runge method for the practical estimation of the rate of convergence s of the scheme with respect to the space variables at a fixed value of t = T . In the case when the exact solution is not known the formula for s is
Ch − Ch/2 s = ln C −C h/2
h/4
ln 2,
where Ch , Ch/2 and Ch/4 are computed values of the numerical solution at grid nodes with step h, h/2 and h/4 respectively. We use three inserted grids with 40 × 20, 80 × 40 and 160 × 80 subintervals respectively by ξ and z in the domain (ξ , z ) ∈ [−1, 1] × [z 0 , 1] and τ = τk = 0.0001. A part of the results from the calculations for the rate of convergence is presented in Table 3. It can be seen from Table 3 that the rate of convergence is very high when the nodes are not near to the point (0, zp ) where the Dirac-delta functions are concentrated.
12
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
Fig. 3. Exact solution at t = T . Table 3 Runge method for the rate of convergence, θ = 0.5. z
0.0002 0.05 0.10 0.15 0.20 0.25 0.30
ξ −0.20
−0.15
−0.10
−0.05
0
7.12 6.91 0.48 1.23 3.44 5.43 7.38
6.91 4.73 0.16 1.37 3.63 5.56 7.34
6.71 3.16 −0.16 1.59 3.85 5.70 7.39
6.50 2.67 −0.46 1.89 4.05 5.82 7.43
6.30 2.92 −0.74 2.23 4.16 5.86 7.44
0.05
0.10
0.15
0.20
6.08 2.42 −0.46 1.89 4.05 5.82 7.41
5.86 2.02 −0.16 1.59 3.84 5.70 7.33
5.64 1.72 0.16 1.35 3.62 5.53 7.22
5.42 1.48 0.51 1.16 3.39 5.34 7.07
6. Conclusions Here we deal with a domain transformation method, often used in computational mechanics. We study the transformed degenerate parabolic problem modeling settling particles dispersion for non-negativity. We construct and discuss a fitted finite volume difference scheme for the problem transformed on the finite domain Ω = (−1, 1) × (z 0 , 1). We may use any interval (−l, l), l > 0 instead of (−1, 1). This leads to insignificant computational cost of the solution. We have shown that the numerical scheme results in a monotone numerical approximation, i.e., it conserves the NP of the continuous solution. Experiments demonstrate the efficiency of our scheme. Besides the full time-dependent model presented here, it is interesting to consider the analogous steady problem. We intend to do this in the future. Acknowledgments The authors are grateful to the referees for their comment and suggestions which helped improve the quality of the manuscript. The first author is supported by the Sofia University Foundation under Grant No. 106/2013. The second author is supported by the Bulgarian Fund of Sciences under Grant No. DID 02/37/09. References [1] B.S. Mazumder, S. Bandyopadhyay, On solute dispersion from an elevated line source in an open-channel flow, J. Eng. Math. 40 (2001) 197–209. [2] K.K. Mondal, B.S. Mazumder, On dispersion of settling particles from an elevated source in an open-channel flow, J. Comput. Appl. Math. 193 (2006) 22–37. [3] K.K. Mondal, B.S. Mazumder, Dispersion of fine settling particles from an elevated source in an oscillatory turbulent flow, Eur. J. Mech. B Fluids 27 (6) (2008) 707–725. [4] N. Madden, K.K. Mondal, Improved mathematical and numerical modelling of dispersion of a solute from a continuous source, in: Lect. Notes Comp. Sci. Enj., vol. 81, Springer, 2011, pp. 177–185. [5] I. Nezu, H. Nakagawa, Turbulence in open-chanel flows, in: IAHR Monograph, A.A. Balkema Publishers, Rotterdam, 1993, p. 281. [6] P.J. Sullivan, H. Yip, Near-field contaminant dispersion from an elevated line-source, Z. Angew. Math. Phys. 38 (1987) 409–423. [7] K.W. Morton, Numerical Solution of Convection–Diffusion Problems, Chapman and Hall, London, 1996. [8] T. Chernogorova, R. Valkov, Finite volume difference scheme for a degenerate parabolic equation in the zero-coupon bond pricing, Math. Comput. Modelling 54 (2011) 2659–2671. [9] J.J.H. Miller, S. Wang, An exponentially fitted element volume method for the numerical solution of 2D unsteady incompressible flow problems, J. Comput. Phys. 115 (1994) 56–64. [10] S. Wang, A novel finite volume method for Black–Scholes equation governing option pricing, IMA J. Numer. Anal. 24 (2004) 699–720.
T. Chernogorova, L. Vulkov / Computers and Mathematics with Applications (
)
–
13
[11] O.A. Oleinik, E.V. Radkevich, Second Order Equations with Nonnegative Characteristic Form, Plenum, Press, New York, 1973, MR0457908. [12] T. Apel, O. Benedix, D. Sirch, B. Vexler, A priori mesh grading for an elliptic problem with Dirac right hand side, SIAM J. Numer. Anal. 49 (3) (2011) 992–1005. [13] X. Yang, X. Zhang, Z. Li, Guo-Wei He, A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations, J. Comput. Phys. 228 (2009) 7821–7836. [14] A. Kufner, Weighted Sobolev Spaces, John Wiley and Sons, New York, 1985. [15] J. Haslinger, M. Miettinen, P.D. Panagiotopoulos, Finite Element Method for Hemivariational Inequalities, Kluwer, Dordrecht, The Netherlands, 1999. [16] I. Farago, I. Horvath, Discrete maximum principle and adequate discretizations of linear parabolic problems, SIAM Sci. Comput. 28 (2006) 2313–2336. [17] D. Gilbert, N. Trudinger, Elliptic Partial Differential Equations of Second Order, 32nd ed., Springer, 2011. [18] O.A. Ladyzenskaja, V.A. Solonnikov, N.N. Ural’tseva, Linear and Quasilinear Equations of Parabolic Type, in: Amer. Math. Soc. Transl. Monographs, vol. 23, 1968. [19] A.A. Samarskii, Theory of Finite Difference Schemes, Marcel Decker, New York, 2003.