- Email: [email protected]
A numerical computation of a non-dimensional form of stream water quality model with hydrodynamic advection–dispersion–reaction equations
Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
Contents lists available at ScienceDirect
Nonlinear Analysis: Hybrid Systems journal homepage: www.elsevier.com/locate/nahs
A numerical computation of a non-dimensional form of stream water quality model with hydrodynamic advection–dispersion–reaction equations Nopparat Pochai ∗ Department of Mathematics and Computer Science, Faculty of Science, King Mongkut’s Institute of Technology Ladkrabang, Bangkok, 10520, Thailand
article
info
Article history: Received 21 May 2009 Accepted 2 June 2009 Keywords: Finite differences Crank–Nicolson scheme Backward time central space scheme One-dimensional Hydrodynamic model Advection-dispersion-reaction equations
abstract Mathematical models of water quality assessment problems often arise in environmental science. The modelling often involves numerical methods to solve the equations. In this research, two mathematical models are used to simulate pollution due to sewage effluent in the nonuniform flow of water in a stream with varied current velocity. The first is a hydrodynamic model that provides the velocity field and elevation of the water flow. The second is a dispersion model, where the commonly used governing factor is the one-dimensional advection–dispersion–reaction equation that gives the pollutant concentration fields. In the simulation processes, we used the Crank–Nicolson method system of a hydrodynamic model and the backward time central space scheme for the dispersion model. Finally, we present a numerical simulation that confirms the results of the techniques. © 2009 Elsevier Ltd. All rights reserved.
Methods for detecting the amount of pollutants both in the air and water mostly are conducted by field measurement and a mathematical simulation. For the shallow water mass transport problems that are presented in [1], the method of characteristics has been reported as being applied with success, but it presents in real cases some difficulties. In [2], the finite element method for solving steady water pollution control to achieve the minimum cost is presented. The numerical techniques for solving the uniform flow of stream water quality model, the one-dimensional advection-diffusion-reaction equation in presented in [3–6]. Most nonuniform flow models require data concerned with the velocity of the current at any point and any time in the domain. The hydrodynamics model provides the velocity field and tidal elevation of the water. In [7,8], the hydrodynamics model and convection-diffusion equation were used to approximate the velocity of the water current in a bay and a channel, respectively. In this paper, the results from hydrodynamic model are data for the nonuniform flow of the advection-diffusion-reaction equation, which provides the pollutant concentration field. The term of the friction forces due to the drag of sides of the stream is considered. We also find the theoretical solution of the model at the end point of the domain. We then use the analytical solution to check the accuracy of our approximate solution. The stream has a simple one space dimension as shown in Fig. 1. Averaging the equation over the depth, and discarding the term due to the Coriolis force, it follows that the one-dimensional shallow water and advection-diffusion-reaction equations are applicable. We use the Crank–Nicolson method and the backward time central space method to approximate the velocity, the tidal elevation and the concentration from the foresaid models.
∗
Tel.: +66 2 326 4341; fax: +66 2 326 4341. E-mail address: [email protected].
1751-570X/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.nahs.2009.06.003
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
667
u(x,t)
z=ζ
z
z=0 x h(x)
Fig. 1. The shallow water system.
1. Mathematical model 1.1. The hydrodynamic model The continuity and momentum equations govern the hydrodynamic behavior of the stream [11,8]. Averaging the equations over the depth, discarding the term due to the Coriolis parameter, shearing stresses and surface wind. We now introduce the one-dimensional shallow water equations:
∂ζ ∂ + [(h + ζ )u] = 0, ∂t ∂x ∂u ∂ζ +g =0 ∂t ∂x
(1) (2)
where x is longitudinal distance along the stream (m), t is the time (s), h(x) is the depth measured from the mean water level to the stream bed (m), ζ (x, t ) is the elevation from the mean water level to the temporary water surface or the tidal elevation (m/s), and u(x, t ) are the velocity components (m/s), for all x ∈ [0, l]. We assume that h is a constant and ζ h. Then the Eqs. (1) and (2) leads to
∂u . ∂ζ +h = 0, ∂t ∂x ∂u ∂ζ +g = 0. ∂t ∂x
(3) (4)
√
We will consider the equation in the dimensionless problem by letting U = u/ gh, X = x/l, Z = ζ /h and T = t Substituting them into the Eqs. (3) and (4) leads to
∂Z ∂U + = 0, ∂T ∂X ∂U ∂Z + = 0. ∂T ∂X
√
gh/l.
(5) (6)
We now introduce a damping term into Eqs. (5)–(6) to represent frictional forces due to the drag of the sides of the stream, thus :
∂Z ∂U + = 0, ∂T ∂X ∂Z ∂U + = −U ∂T ∂X
(7) (8)
with the initial conditions at t = 0 and 0 ≤ X ≤ 1 being specified: Z = 0 and U = 0. The boundary conditions for t > 0 are specified: Z = eit at X = 0 and ∂∂ XZ = 0 at X = 1. The system of Eqs. (7)–(8) is called the damped equation. We solve the damped equation by using the finite difference method. In order to solve the Eqs. (7)–(8) in [0, 1] × [0, T ], for convenience using u, d for U and Z respectively,
∂u ∂d + = −u, ∂t ∂x ∂d ∂u + =0 ∂t ∂x with the initial conditions u = 0, d = 0 at t = 0 and the boundary conditions d(0, t ) = f (t ) and ∂∂ dx = 0 at x = 1.
(9) (10)
668
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
1.2. Dispersion model In a stream water quality model, the governing equations are the dynamic one-dimensional advection-diffusion-reaction equations (ADRE). A simplified representation by averaging the equation over the depth is shown in [3–6] as
∂C ∂C ∂ 2C +u = D 2 − KC (11) ∂t ∂x ∂x where C (x, t ) is the concentration averaged by depth at the point x and at time t, D is the diffusion coefficient, K is the mass decay rate, and u(x, t ) is the velocity component, for all x ∈ [0, 1]. We will consider the model with following conditions. The initial condition C (x, 0) = 0 at t = 0 for all x > 0. The boundary conditions C (0, t ) = C0 at x = 0 and ∂∂Cx = 0 at x = 1 where C0 is a constant.
2. Numerical experiment The hydrodynamic model provides the velocity field and elevation of the water. The calculated results of the model are input to the dispersion model, which provides the pollutant concentration field. 2.1. Numerical solution of the hydrodynamic model We will follow the numerical techniques of [8]. To find the water velocity and water elevation from Eqs. (9)–(10), we make the following change of variable, v = et u and substituting them into Eqs. (9)–(10), we have
∂v ∂d + e− t = 0, ∂t ∂x ∂d ∂v + e− t = 0. ∂t ∂x
(12) (13)
The Eqs. (12)–(13) can be written in matrix form
v d
et 0
0 + −t e
t
v d
= x
0 . 0
(14)
That is Ut + AUx = 0¯ ,
(15)
where 0 e− t
et , 0
v
and
A=
U =
d
(16)
v d
= t
∂v/∂ t ∂ d/∂ t
(17)
with the initial condition d = v = 0 at t = 0. The left boundary conditions for x = 0, t > 0 are specified: d(0, t ) = sin t and ∂v = −et cos t, and the right boundary conditions for x = 1, t > 0 are specified: ∂∂ dx = 0 and v(0, t ) = 0. ∂x We now discretize Eq. (15) by dividing the interval [0, 1] into M subintervals such that M ∆x = 1 and the interval [0, T ] into N subintervals such that N ∆t = T . We can then approximate d(xi , tn ) by dni , the value of the difference approximation of d(x, t ) at point x = i∆x and t = n∆t, where 0 ≤ i ≤ M and 0 ≤ n ≤ N, and similarly defined for vin and Uin . The grid points (xn , tn ) are defined by xi = i∆x for all i = 0, 1, 2, . . . , M and tn = n∆t for all n = 0, 1, 2, . . . , N in which M and N are positive integers. Using the Crank–Nicolson method [9] to Eq. (15), we get the following finite difference equation:
I−
1 4
λA[(∆x + ∇x )]
Uin+1
1
= I + λA(∆x + ∇x ) Uin 4
(18)
where Uin =
n vi dni
,
∆x Uin = Uin+1 − Uin and ∇x Uin = Uin − Uin−1 ,
(19)
I is the unit matrix of order 2 and λ = ∆t /∆x. Applying the initial and boundary conditions given for Eqs. (12)–(13) we obtained the general form Gn+1 U¯ n+1 = E n U¯ n + F n
(20)
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
669
where
Gn+1
1
λ n+1 a2 4 0 λ = a2n+1 4 . ..
λ 4
0
0
λ − an1+1
0
0
1
λ − an2+1
0
0
0
4
an1+1
1
0
0
λ − a1n+1
0
0
1
λ − a2n+1
0
..
..
..
..
.
0
1
λ n − a2 4 0 λ n E = − an 4 2 . ..
.
0
0
4
.
λ
0
λ
0
4
4
0
an2
0
0
0
1
0
0
0
0
1
..
..
..
1
λ − an1
.
λ 4
4
.
.
λ 4
an2
..
.
λ
..
a1n+1
1
0
0
4
an1
0
..
.
0
0
0
λ − an1
1
0
0
λ − an2
0
0
4
.
an2+1 0
0
4
4
λ − an1
0
4
4
.
, λ n+1 − a1 4
(21)
1
, λ n a1 4
U1n+1 U n +1 2
. , ..
U¯ n =
(22)
n+1 UM −1
1
λ λ − an1+1 sin(tn+1 ) − an1 sin(tn ) 4 4 λ λ n +1 −tn+1 − a2 ∆xe cos(tn+1 ) − an2 ∆e−tn cos(tn ) 4 4 0 Fn = 0 .. .
(23)
0 0
where an1 = etn , an2 = e−tn and tn = n∆t for all n = 0, 1, 2, . . . , N.
2.2. Numerical solution of the dynamic advection-diffusion-reaction equation We can then approximate C (xi , tn ) by Cin , the value of the difference approximation of C (x, t ) at point x = i∆x and t = n∆t, where 0 ≤ i ≤ M and 0 ≤ n ≤ N. The grid points (xn , tn ) are defined by xi = i∆x for all i = 0, 1, 2, . . . , M and tn = n∆t for all n = 0, 1, 2, . . . , N in which M and N are positive integers. Taking the backward time and central space technique [9] into Eq. (11), we get the following discretization: C ∼ = Cin+1
∂C ∼ = ∂t
(24)
Cin+1
− ∆t
Cin
∂ C ∼ Cin++11 − Cin−+11 = ∂x 2∆x ∂ 2 C ∼ Cin++11 − 2Cin+1 + Cin−+11 = ∂ x2 (∆x)2 U ∼ U n+1 . =b i
(25) (26) (27) (28)
670
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
Substitute Eqs. (24)–(28) into Eq. (11), we get Cin+1 − Cin
∆t
Uin+1
+
Cin++11 − Cin−+11
!
Cin++11 − 2Cin+1 + Cin−+11
=D
2∆x
! − KCin+1
(∆x)2
for 1 ≤ i ≤ M and 0 ≤ n ≤ N − 1. Let λ = (D∆∆x)t2 and γin+1 =
∆t bn+1 U , ∆x i
(29)
Eq. (29) becomes
1 n +1 1 n +1 n +1 n+1 Ci−1 + (1 + 2λ + K ∆t )Ci − λ − γi Cin++11 = Cin . − λ + γi 2
(30)
2
For i = 1, move the known value of the left boundary C0n+1 in Eq. (30) to the right-hand side and rearranging, we obtain
(1 + 2λ + K ∆t )
C1n+1
1
− λ− γ 2
n +1 1
C2n+1
=
C1n
1
+ λ+ γ 2
n +1 1
C0n+1 .
(31)
n +1 n+1 n +1 For i = M, substitute the approximate unknown value of the right boundary by CM +1 = 2CM − CM −1 and rearranging, we obtain
− γMn+1 CMn+−11 + (1 + K ∆t + γMn+1 )CMn+1 = CMn .
(32)
We then have Eq. (30) in the general form
− → n +1 C
1 − →n − → n +1 = Hn−1 C + λ + γ1n+1 Hn−1 C 0
(33)
2
for n = 1, 2, . . . , N − 1, where
(1 + 2λ + K ∆t ) − λ + 1 γ n+1 2 i Hn =
1 − λ − γin+1
2
( 1 + 2 λ + K ∆t ) ..
− λ− γ 2
..
.
1
− λ+ γ 2
1
n+1 i
n +1 i
.
(1 + 2λ + K ∆t ) −γMn+1
C1n n C2
− →n C
= .. . n CM
M ×1
C0n+1 0
,
− → n +1 C
0
.. . 1 − λ − γin+1 2
(1 + K ∆t + γMn+1 )
M ×M
= .. . 0
.
(34)
M ×1
Both Crank–Nicolson and backward time central space methods are unconditionally stable [9,10]. 3. The accuracy of the hydrodynamic approximation It is not hard to find the analytical solution d(x, t ) in the Eqs. (9)–(10) with f (t ) = sin t. By changing variables, d(x, t ) = eit D (x) and u(x, t ) = eit U(x) for some D (x), U(x) ∈ C02 [0, 1] and substitute into Eqs. (9)–(10). Using a separable variables technique, we can obtain d(1, t ) a solution d(1, t ) =
sin t cos β cosh α − cos t sin β sinh α cos2 β cosh2 α + sin2 β sinh2 α
(35)
where α = 21/4 cos(3π /8) and β = 21/4 sin(3π /8). Anyhow, it is not easy to find the analytical solution u(x, t ) of Eqs. (9)–(10). We use the solution d(1, t ) obtained in Eq. (35) to verify to our approximate solution obtained by the Crank–Nicolson method Eq. (20). Actually when using the Crank–Nicolson method, we get the approximate solution both d(x, t ) and u(x, t ). We assume that when we get a good approximation for d(x, t ) this implies that the method gives a good approximation for u(x, t ). The verification of the approximate solution d¯ (1, t ) is shown in Table 1 and Fig. 2. Fig. 2 shows the comparison between the analytical solutions d(1, t ) and the approximate solutions d¯ (1, t ) only at the end of the domain. Table 1 shows that an estimate of the maximum error is less than 7.0%.
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
671
2 Height of water elevation d(1,t)
Analytical solution Numerical solution
1.5 1 0.5 0 -0.5 -1 -1.5 -2
0
100
200
300
400
500
600
700
800
900
1000
Number of time step with Δt = 0.05
Fig. 2. Comparison of analytical solution for height of water elevation with results obtained by numerical technique at the end point of the domain. Table 1 The error defined by error (T ) = max|d(1, t ) − d¯ (1, t )| for all T − π ≤ t ≤ T + π at fixed ∆t /∆x = 0.25. T
∆x
∆t
d(1, t )
d¯ (1, t )
error (T )
20
0.200 0.100 0.050 0.025 0.200 0.100 0.050 0.025 0.200 0.100 0.050 0.025
0.05000 0.02500 0.01250 0.00625 0.05000 0.02500 0.01250 0.00625 0.05000 0.02500 0.01250 0.00625
1.50158 0.09392 −0.25831 −0.33323 1.49666 −0.09358 0.25798 0.33290 1.50146 0.09324 −0.25765 −0.33257
1.45659 0.38856 −0.40244 −0.40433 1.45165 −0.38821 0.40198 0.40389 1.45644 0.38787 −0.40165 −0.40356
0.04499 0.29465 0.14413 0.07110 0.04501 0.29463 0.14401 0.07099 0.04502 0.29463 0.14401 0.07099
30
40
Table 2 The velocity of water flow u(x, t ). t
10 20 30 40 50
x 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
−0.5478
−0.5695
−0.5666
−0.5199
−0.4711
−0.4061
−0.3352
−0.2571
−0.1747
−0.0878
1.3101 −0.4468 −1.0361 1.0939
1.2213 −0.3731 −0.9898 0.9918
1.1232 −0.3085 −0.9258 0.8867
1.0160 −0.2527 −0.8459 0.7791
0.8931 −0.2057 −0.7513 0.6700
0.7504 −0.1651 −0.6439 0.5594
0.6097 −0.1252 −0.5258 0.4479
0.4595 −0.0875 −0.4004 0.3356
0.3099 −0.0573 −0.2711 0.2233
0.1554 −0.0276 −0.1367 0.1114
0.0000 0.0000 0.0000 0.0000 0.0000
4. Application to the stream water quality assessment problem Suppose we want to measure the pollutant concentration C of a uniform stream. A stream is aligned with longitudinal distance, 1.0 km total length and 1.0 m depth. There is a plant which discharges waste water into the stream and the pollutant concentration at discharge point is C (0, t ) = C0 = 1 mg/L at x = 0 for all t > 0 and C (x, 0) = 0 at t = 0. The elevation of water at the discharge point can be described as a function d(0, t ) = f (t ) = sin t m for all t > 0, and the elevation does not change at x = 1.0 km. The stream system characteristics are: the diffusion coefficient D is 0.05 m2 /s, and the first-order reaction rate is 10−5 s−1 . In the analysis conducted in this study, a mesh with 80 elements with ∆x = 0.0125 and time increment is 1.6 s with ∆t = 0.005, characterizing a one-dimensional flow. Using the Eq. (20), we get the water velocity u(x, t ) on Table 2 and Fig. 3. Next, we can plug the velocity of the water into Eq. (33)- (34), we then have the pollutant concentration C in Table 3 and Fig. 4. 5. Discussion and conclusions In this paper, we can combine the hydrodynamic model and the convection-diffusion-reaction equation for an approximation of the pollutant concentration in a stream when the velocity of the current is not uniform. With the technique developed in this paper the response of the stream to the different external inputs: the elevation of water and the pollutant concentration at the discharge point can be obtained. It is also possible to find tentative better locations and periods of time of different discharge points to the stream. This model can be applied to the real cases for sewage effluent in a canal that connected to another stream or reservoir.
672
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
Table 3 The pollutant concentration C (x, t ). t
x 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
10 20 30 40 50
1.0000 1.0000 1.0000 1.0000 1.0000
0.3211 0.9996 0.9960 0.2842 0.9922
0.1128 0.9972 0.9899 0.1808 0.9639
0.0491 0.9876 0.9772 0.1481 0.8950
0.0267 0.9589 0.9514 0.1244 0.7719
0.0167 0.8909 0.9032 0.1024 0.6027
0.0111 0.7655 0.8209 0.0814 0.4184
0.0074 0.5847 0.6933 0.0610 0.2560
0.0046 0.3790 0.5143 0.0411 0.1366
0.0023 0.1869 0.2883 0.0216 0.0588
0.0002 0.0201 0.0327 0.0024 0.0060
Velocity of water (m/s)
1.5 1 0.5 0 -0.5 -1 -1.5 0 10
20 30
40
50 60 0
70
80
90 0
2500 3000 1500 2000 500 1000 0
3500
4000
Fig. 3. The water velocity u(x, t ) at t = 1 h 46 min.
Fig. 4. The pollutant concentration C (x, t ) at t = 1 h 46 min.
Acknowledgements The author would like to thanks: S. Tangmanee, L.J. Crane and J.J.H. Miller for their valuable comments. References [1] A. Garzon, L. D’Alpaos, A modified method of the characteristic technique combined with Gelerkin finite element method to solve shallow water mass transport problems, in: Proceedings 23rd International Conference in Coastal Engineering, vol. 3, 1992, pp. 3068–3080. [2] N. Pochai, S. Tangmanee, L.J. Crane, J.J.H. Miller, A mathematical model of water pollution control using the finite element method, Proceedings in Applied Mathematics and Mechanics 6 (1) (2006) 755–756. [3] J.Y. Chen, C.H. Ko, S. Bhattacharjee, M. Elimelech, Role of spatial distribution of porous medium surface charge heterogeneity in colloid transport, Colloids and Surfaces A: Physicochemical and Engineering Aspects 191 (12) (2001) 3–15. [4] G. Li, C.R. Jackson, Simple, accurate and efficient revisions to MacCormack and Saulyev schemes: High Peclet numbers, Applied Mathematics and Computation 186 (2007) 610–622. [5] E.M. O’Loughlin, K.H. Bowmer, Dilution and decay of aquatic herbicides in flowing channels, Journal of Hydrology 26 (34) (1975) 217–235. [6] A.I. Stamou, Improving the numerical modeling of river water quality by using high order difference schemes, Water Research 26 (12) (1992) 1563–1570. [7] P. Tabuenca, J. Vila, J. Cardona, A. Samartin, Finite element simulation of dispersion in the bay of Santander, Advanced in Engineering Software 28 (1997) 313–332.
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
673
[8] N. Pochai, S. Tangmanee, L.J. Crane, J.J.H. Miller, A water quality computation in the uniform channel, Journal of Interdisciplinary Mathematics 11 (6) (2008) 803–814. [9] A.R. Mitchell, Computational Methods in Partial Differential Equations, John-Wiley and Sons, 1969. [10] S.C. Chapra, Surface Water-Quality Modeling, McGraw-Hill, NewYork, 1997. [11] H. Ninomiya, K. Onishi, Flow analysis using a PC, in: Computational Mechanics Publications, CRC Press, Boca Raton, 1991.
Contents lists available at ScienceDirect
Nonlinear Analysis: Hybrid Systems journal homepage: www.elsevier.com/locate/nahs
A numerical computation of a non-dimensional form of stream water quality model with hydrodynamic advection–dispersion–reaction equations Nopparat Pochai ∗ Department of Mathematics and Computer Science, Faculty of Science, King Mongkut’s Institute of Technology Ladkrabang, Bangkok, 10520, Thailand
article
info
Article history: Received 21 May 2009 Accepted 2 June 2009 Keywords: Finite differences Crank–Nicolson scheme Backward time central space scheme One-dimensional Hydrodynamic model Advection-dispersion-reaction equations
abstract Mathematical models of water quality assessment problems often arise in environmental science. The modelling often involves numerical methods to solve the equations. In this research, two mathematical models are used to simulate pollution due to sewage effluent in the nonuniform flow of water in a stream with varied current velocity. The first is a hydrodynamic model that provides the velocity field and elevation of the water flow. The second is a dispersion model, where the commonly used governing factor is the one-dimensional advection–dispersion–reaction equation that gives the pollutant concentration fields. In the simulation processes, we used the Crank–Nicolson method system of a hydrodynamic model and the backward time central space scheme for the dispersion model. Finally, we present a numerical simulation that confirms the results of the techniques. © 2009 Elsevier Ltd. All rights reserved.
Methods for detecting the amount of pollutants both in the air and water mostly are conducted by field measurement and a mathematical simulation. For the shallow water mass transport problems that are presented in [1], the method of characteristics has been reported as being applied with success, but it presents in real cases some difficulties. In [2], the finite element method for solving steady water pollution control to achieve the minimum cost is presented. The numerical techniques for solving the uniform flow of stream water quality model, the one-dimensional advection-diffusion-reaction equation in presented in [3–6]. Most nonuniform flow models require data concerned with the velocity of the current at any point and any time in the domain. The hydrodynamics model provides the velocity field and tidal elevation of the water. In [7,8], the hydrodynamics model and convection-diffusion equation were used to approximate the velocity of the water current in a bay and a channel, respectively. In this paper, the results from hydrodynamic model are data for the nonuniform flow of the advection-diffusion-reaction equation, which provides the pollutant concentration field. The term of the friction forces due to the drag of sides of the stream is considered. We also find the theoretical solution of the model at the end point of the domain. We then use the analytical solution to check the accuracy of our approximate solution. The stream has a simple one space dimension as shown in Fig. 1. Averaging the equation over the depth, and discarding the term due to the Coriolis force, it follows that the one-dimensional shallow water and advection-diffusion-reaction equations are applicable. We use the Crank–Nicolson method and the backward time central space method to approximate the velocity, the tidal elevation and the concentration from the foresaid models.
∗
Tel.: +66 2 326 4341; fax: +66 2 326 4341. E-mail address: [email protected].
1751-570X/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.nahs.2009.06.003
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
667
u(x,t)
z=ζ
z
z=0 x h(x)
Fig. 1. The shallow water system.
1. Mathematical model 1.1. The hydrodynamic model The continuity and momentum equations govern the hydrodynamic behavior of the stream [11,8]. Averaging the equations over the depth, discarding the term due to the Coriolis parameter, shearing stresses and surface wind. We now introduce the one-dimensional shallow water equations:
∂ζ ∂ + [(h + ζ )u] = 0, ∂t ∂x ∂u ∂ζ +g =0 ∂t ∂x
(1) (2)
where x is longitudinal distance along the stream (m), t is the time (s), h(x) is the depth measured from the mean water level to the stream bed (m), ζ (x, t ) is the elevation from the mean water level to the temporary water surface or the tidal elevation (m/s), and u(x, t ) are the velocity components (m/s), for all x ∈ [0, l]. We assume that h is a constant and ζ h. Then the Eqs. (1) and (2) leads to
∂u . ∂ζ +h = 0, ∂t ∂x ∂u ∂ζ +g = 0. ∂t ∂x
(3) (4)
√
We will consider the equation in the dimensionless problem by letting U = u/ gh, X = x/l, Z = ζ /h and T = t Substituting them into the Eqs. (3) and (4) leads to
∂Z ∂U + = 0, ∂T ∂X ∂U ∂Z + = 0. ∂T ∂X
√
gh/l.
(5) (6)
We now introduce a damping term into Eqs. (5)–(6) to represent frictional forces due to the drag of the sides of the stream, thus :
∂Z ∂U + = 0, ∂T ∂X ∂Z ∂U + = −U ∂T ∂X
(7) (8)
with the initial conditions at t = 0 and 0 ≤ X ≤ 1 being specified: Z = 0 and U = 0. The boundary conditions for t > 0 are specified: Z = eit at X = 0 and ∂∂ XZ = 0 at X = 1. The system of Eqs. (7)–(8) is called the damped equation. We solve the damped equation by using the finite difference method. In order to solve the Eqs. (7)–(8) in [0, 1] × [0, T ], for convenience using u, d for U and Z respectively,
∂u ∂d + = −u, ∂t ∂x ∂d ∂u + =0 ∂t ∂x with the initial conditions u = 0, d = 0 at t = 0 and the boundary conditions d(0, t ) = f (t ) and ∂∂ dx = 0 at x = 1.
(9) (10)
668
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
1.2. Dispersion model In a stream water quality model, the governing equations are the dynamic one-dimensional advection-diffusion-reaction equations (ADRE). A simplified representation by averaging the equation over the depth is shown in [3–6] as
∂C ∂C ∂ 2C +u = D 2 − KC (11) ∂t ∂x ∂x where C (x, t ) is the concentration averaged by depth at the point x and at time t, D is the diffusion coefficient, K is the mass decay rate, and u(x, t ) is the velocity component, for all x ∈ [0, 1]. We will consider the model with following conditions. The initial condition C (x, 0) = 0 at t = 0 for all x > 0. The boundary conditions C (0, t ) = C0 at x = 0 and ∂∂Cx = 0 at x = 1 where C0 is a constant.
2. Numerical experiment The hydrodynamic model provides the velocity field and elevation of the water. The calculated results of the model are input to the dispersion model, which provides the pollutant concentration field. 2.1. Numerical solution of the hydrodynamic model We will follow the numerical techniques of [8]. To find the water velocity and water elevation from Eqs. (9)–(10), we make the following change of variable, v = et u and substituting them into Eqs. (9)–(10), we have
∂v ∂d + e− t = 0, ∂t ∂x ∂d ∂v + e− t = 0. ∂t ∂x
(12) (13)
The Eqs. (12)–(13) can be written in matrix form
v d
et 0
0 + −t e
t
v d
= x
0 . 0
(14)
That is Ut + AUx = 0¯ ,
(15)
where 0 e− t
et , 0
v
and
A=
U =
d
(16)
v d
= t
∂v/∂ t ∂ d/∂ t
(17)
with the initial condition d = v = 0 at t = 0. The left boundary conditions for x = 0, t > 0 are specified: d(0, t ) = sin t and ∂v = −et cos t, and the right boundary conditions for x = 1, t > 0 are specified: ∂∂ dx = 0 and v(0, t ) = 0. ∂x We now discretize Eq. (15) by dividing the interval [0, 1] into M subintervals such that M ∆x = 1 and the interval [0, T ] into N subintervals such that N ∆t = T . We can then approximate d(xi , tn ) by dni , the value of the difference approximation of d(x, t ) at point x = i∆x and t = n∆t, where 0 ≤ i ≤ M and 0 ≤ n ≤ N, and similarly defined for vin and Uin . The grid points (xn , tn ) are defined by xi = i∆x for all i = 0, 1, 2, . . . , M and tn = n∆t for all n = 0, 1, 2, . . . , N in which M and N are positive integers. Using the Crank–Nicolson method [9] to Eq. (15), we get the following finite difference equation:
I−
1 4
λA[(∆x + ∇x )]
Uin+1
1
= I + λA(∆x + ∇x ) Uin 4
(18)
where Uin =
n vi dni
,
∆x Uin = Uin+1 − Uin and ∇x Uin = Uin − Uin−1 ,
(19)
I is the unit matrix of order 2 and λ = ∆t /∆x. Applying the initial and boundary conditions given for Eqs. (12)–(13) we obtained the general form Gn+1 U¯ n+1 = E n U¯ n + F n
(20)
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
669
where
Gn+1
1
λ n+1 a2 4 0 λ = a2n+1 4 . ..
λ 4
0
0
λ − an1+1
0
0
1
λ − an2+1
0
0
0
4
an1+1
1
0
0
λ − a1n+1
0
0
1
λ − a2n+1
0
..
..
..
..
.
0
1
λ n − a2 4 0 λ n E = − an 4 2 . ..
.
0
0
4
.
λ
0
λ
0
4
4
0
an2
0
0
0
1
0
0
0
0
1
..
..
..
1
λ − an1
.
λ 4
4
.
.
λ 4
an2
..
.
λ
..
a1n+1
1
0
0
4
an1
0
..
.
0
0
0
λ − an1
1
0
0
λ − an2
0
0
4
.
an2+1 0
0
4
4
λ − an1
0
4
4
.
, λ n+1 − a1 4
(21)
1
, λ n a1 4
U1n+1 U n +1 2
. , ..
U¯ n =
(22)
n+1 UM −1
1
λ λ − an1+1 sin(tn+1 ) − an1 sin(tn ) 4 4 λ λ n +1 −tn+1 − a2 ∆xe cos(tn+1 ) − an2 ∆e−tn cos(tn ) 4 4 0 Fn = 0 .. .
(23)
0 0
where an1 = etn , an2 = e−tn and tn = n∆t for all n = 0, 1, 2, . . . , N.
2.2. Numerical solution of the dynamic advection-diffusion-reaction equation We can then approximate C (xi , tn ) by Cin , the value of the difference approximation of C (x, t ) at point x = i∆x and t = n∆t, where 0 ≤ i ≤ M and 0 ≤ n ≤ N. The grid points (xn , tn ) are defined by xi = i∆x for all i = 0, 1, 2, . . . , M and tn = n∆t for all n = 0, 1, 2, . . . , N in which M and N are positive integers. Taking the backward time and central space technique [9] into Eq. (11), we get the following discretization: C ∼ = Cin+1
∂C ∼ = ∂t
(24)
Cin+1
− ∆t
Cin
∂ C ∼ Cin++11 − Cin−+11 = ∂x 2∆x ∂ 2 C ∼ Cin++11 − 2Cin+1 + Cin−+11 = ∂ x2 (∆x)2 U ∼ U n+1 . =b i
(25) (26) (27) (28)
670
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
Substitute Eqs. (24)–(28) into Eq. (11), we get Cin+1 − Cin
∆t
Uin+1
+
Cin++11 − Cin−+11
!
Cin++11 − 2Cin+1 + Cin−+11
=D
2∆x
! − KCin+1
(∆x)2
for 1 ≤ i ≤ M and 0 ≤ n ≤ N − 1. Let λ = (D∆∆x)t2 and γin+1 =
∆t bn+1 U , ∆x i
(29)
Eq. (29) becomes
1 n +1 1 n +1 n +1 n+1 Ci−1 + (1 + 2λ + K ∆t )Ci − λ − γi Cin++11 = Cin . − λ + γi 2
(30)
2
For i = 1, move the known value of the left boundary C0n+1 in Eq. (30) to the right-hand side and rearranging, we obtain
(1 + 2λ + K ∆t )
C1n+1
1
− λ− γ 2
n +1 1
C2n+1
=
C1n
1
+ λ+ γ 2
n +1 1
C0n+1 .
(31)
n +1 n+1 n +1 For i = M, substitute the approximate unknown value of the right boundary by CM +1 = 2CM − CM −1 and rearranging, we obtain
− γMn+1 CMn+−11 + (1 + K ∆t + γMn+1 )CMn+1 = CMn .
(32)
We then have Eq. (30) in the general form
− → n +1 C
1 − →n − → n +1 = Hn−1 C + λ + γ1n+1 Hn−1 C 0
(33)
2
for n = 1, 2, . . . , N − 1, where
(1 + 2λ + K ∆t ) − λ + 1 γ n+1 2 i Hn =
1 − λ − γin+1
2
( 1 + 2 λ + K ∆t ) ..
− λ− γ 2
..
.
1
− λ+ γ 2
1
n+1 i
n +1 i
.
(1 + 2λ + K ∆t ) −γMn+1
C1n n C2
− →n C
= .. . n CM
M ×1
C0n+1 0
,
− → n +1 C
0
.. . 1 − λ − γin+1 2
(1 + K ∆t + γMn+1 )
M ×M
= .. . 0
.
(34)
M ×1
Both Crank–Nicolson and backward time central space methods are unconditionally stable [9,10]. 3. The accuracy of the hydrodynamic approximation It is not hard to find the analytical solution d(x, t ) in the Eqs. (9)–(10) with f (t ) = sin t. By changing variables, d(x, t ) = eit D (x) and u(x, t ) = eit U(x) for some D (x), U(x) ∈ C02 [0, 1] and substitute into Eqs. (9)–(10). Using a separable variables technique, we can obtain d(1, t ) a solution d(1, t ) =
sin t cos β cosh α − cos t sin β sinh α cos2 β cosh2 α + sin2 β sinh2 α
(35)
where α = 21/4 cos(3π /8) and β = 21/4 sin(3π /8). Anyhow, it is not easy to find the analytical solution u(x, t ) of Eqs. (9)–(10). We use the solution d(1, t ) obtained in Eq. (35) to verify to our approximate solution obtained by the Crank–Nicolson method Eq. (20). Actually when using the Crank–Nicolson method, we get the approximate solution both d(x, t ) and u(x, t ). We assume that when we get a good approximation for d(x, t ) this implies that the method gives a good approximation for u(x, t ). The verification of the approximate solution d¯ (1, t ) is shown in Table 1 and Fig. 2. Fig. 2 shows the comparison between the analytical solutions d(1, t ) and the approximate solutions d¯ (1, t ) only at the end of the domain. Table 1 shows that an estimate of the maximum error is less than 7.0%.
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
671
2 Height of water elevation d(1,t)
Analytical solution Numerical solution
1.5 1 0.5 0 -0.5 -1 -1.5 -2
0
100
200
300
400
500
600
700
800
900
1000
Number of time step with Δt = 0.05
Fig. 2. Comparison of analytical solution for height of water elevation with results obtained by numerical technique at the end point of the domain. Table 1 The error defined by error (T ) = max|d(1, t ) − d¯ (1, t )| for all T − π ≤ t ≤ T + π at fixed ∆t /∆x = 0.25. T
∆x
∆t
d(1, t )
d¯ (1, t )
error (T )
20
0.200 0.100 0.050 0.025 0.200 0.100 0.050 0.025 0.200 0.100 0.050 0.025
0.05000 0.02500 0.01250 0.00625 0.05000 0.02500 0.01250 0.00625 0.05000 0.02500 0.01250 0.00625
1.50158 0.09392 −0.25831 −0.33323 1.49666 −0.09358 0.25798 0.33290 1.50146 0.09324 −0.25765 −0.33257
1.45659 0.38856 −0.40244 −0.40433 1.45165 −0.38821 0.40198 0.40389 1.45644 0.38787 −0.40165 −0.40356
0.04499 0.29465 0.14413 0.07110 0.04501 0.29463 0.14401 0.07099 0.04502 0.29463 0.14401 0.07099
30
40
Table 2 The velocity of water flow u(x, t ). t
10 20 30 40 50
x 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
−0.5478
−0.5695
−0.5666
−0.5199
−0.4711
−0.4061
−0.3352
−0.2571
−0.1747
−0.0878
1.3101 −0.4468 −1.0361 1.0939
1.2213 −0.3731 −0.9898 0.9918
1.1232 −0.3085 −0.9258 0.8867
1.0160 −0.2527 −0.8459 0.7791
0.8931 −0.2057 −0.7513 0.6700
0.7504 −0.1651 −0.6439 0.5594
0.6097 −0.1252 −0.5258 0.4479
0.4595 −0.0875 −0.4004 0.3356
0.3099 −0.0573 −0.2711 0.2233
0.1554 −0.0276 −0.1367 0.1114
0.0000 0.0000 0.0000 0.0000 0.0000
4. Application to the stream water quality assessment problem Suppose we want to measure the pollutant concentration C of a uniform stream. A stream is aligned with longitudinal distance, 1.0 km total length and 1.0 m depth. There is a plant which discharges waste water into the stream and the pollutant concentration at discharge point is C (0, t ) = C0 = 1 mg/L at x = 0 for all t > 0 and C (x, 0) = 0 at t = 0. The elevation of water at the discharge point can be described as a function d(0, t ) = f (t ) = sin t m for all t > 0, and the elevation does not change at x = 1.0 km. The stream system characteristics are: the diffusion coefficient D is 0.05 m2 /s, and the first-order reaction rate is 10−5 s−1 . In the analysis conducted in this study, a mesh with 80 elements with ∆x = 0.0125 and time increment is 1.6 s with ∆t = 0.005, characterizing a one-dimensional flow. Using the Eq. (20), we get the water velocity u(x, t ) on Table 2 and Fig. 3. Next, we can plug the velocity of the water into Eq. (33)- (34), we then have the pollutant concentration C in Table 3 and Fig. 4. 5. Discussion and conclusions In this paper, we can combine the hydrodynamic model and the convection-diffusion-reaction equation for an approximation of the pollutant concentration in a stream when the velocity of the current is not uniform. With the technique developed in this paper the response of the stream to the different external inputs: the elevation of water and the pollutant concentration at the discharge point can be obtained. It is also possible to find tentative better locations and periods of time of different discharge points to the stream. This model can be applied to the real cases for sewage effluent in a canal that connected to another stream or reservoir.
672
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
Table 3 The pollutant concentration C (x, t ). t
x 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
10 20 30 40 50
1.0000 1.0000 1.0000 1.0000 1.0000
0.3211 0.9996 0.9960 0.2842 0.9922
0.1128 0.9972 0.9899 0.1808 0.9639
0.0491 0.9876 0.9772 0.1481 0.8950
0.0267 0.9589 0.9514 0.1244 0.7719
0.0167 0.8909 0.9032 0.1024 0.6027
0.0111 0.7655 0.8209 0.0814 0.4184
0.0074 0.5847 0.6933 0.0610 0.2560
0.0046 0.3790 0.5143 0.0411 0.1366
0.0023 0.1869 0.2883 0.0216 0.0588
0.0002 0.0201 0.0327 0.0024 0.0060
Velocity of water (m/s)
1.5 1 0.5 0 -0.5 -1 -1.5 0 10
20 30
40
50 60 0
70
80
90 0
2500 3000 1500 2000 500 1000 0
3500
4000
Fig. 3. The water velocity u(x, t ) at t = 1 h 46 min.
Fig. 4. The pollutant concentration C (x, t ) at t = 1 h 46 min.
Acknowledgements The author would like to thanks: S. Tangmanee, L.J. Crane and J.J.H. Miller for their valuable comments. References [1] A. Garzon, L. D’Alpaos, A modified method of the characteristic technique combined with Gelerkin finite element method to solve shallow water mass transport problems, in: Proceedings 23rd International Conference in Coastal Engineering, vol. 3, 1992, pp. 3068–3080. [2] N. Pochai, S. Tangmanee, L.J. Crane, J.J.H. Miller, A mathematical model of water pollution control using the finite element method, Proceedings in Applied Mathematics and Mechanics 6 (1) (2006) 755–756. [3] J.Y. Chen, C.H. Ko, S. Bhattacharjee, M. Elimelech, Role of spatial distribution of porous medium surface charge heterogeneity in colloid transport, Colloids and Surfaces A: Physicochemical and Engineering Aspects 191 (12) (2001) 3–15. [4] G. Li, C.R. Jackson, Simple, accurate and efficient revisions to MacCormack and Saulyev schemes: High Peclet numbers, Applied Mathematics and Computation 186 (2007) 610–622. [5] E.M. O’Loughlin, K.H. Bowmer, Dilution and decay of aquatic herbicides in flowing channels, Journal of Hydrology 26 (34) (1975) 217–235. [6] A.I. Stamou, Improving the numerical modeling of river water quality by using high order difference schemes, Water Research 26 (12) (1992) 1563–1570. [7] P. Tabuenca, J. Vila, J. Cardona, A. Samartin, Finite element simulation of dispersion in the bay of Santander, Advanced in Engineering Software 28 (1997) 313–332.
N. Pochai / Nonlinear Analysis: Hybrid Systems 3 (2009) 666–673
673
[8] N. Pochai, S. Tangmanee, L.J. Crane, J.J.H. Miller, A water quality computation in the uniform channel, Journal of Interdisciplinary Mathematics 11 (6) (2008) 803–814. [9] A.R. Mitchell, Computational Methods in Partial Differential Equations, John-Wiley and Sons, 1969. [10] S.C. Chapra, Surface Water-Quality Modeling, McGraw-Hill, NewYork, 1997. [11] H. Ninomiya, K. Onishi, Flow analysis using a PC, in: Computational Mechanics Publications, CRC Press, Boca Raton, 1991.