Finite volume discretisation for the one-dimensional convection diffusion-dissipation equation (original) (raw)
Finite volume discretisation for the one-dimensional convection diffusion-dissipation equation
Bienvenu ONDAMI 1∗{ }^{1 *}
1{ }^{1} Université Marien NGOUABI, Brazzaville, CONGO
∗{ }^{*} Corresponding author E-mail:bondami@gmail.com
Abstract
This paper is devoted to analysis of a finite volume scheme for a one-dimensional convection-diffusion-dissipation equation having application in pollution of water table. We analyse a scheme corresponding to a semi-descretization, also called method of lines. Results of numerical experiments using this approach are reported.
Keywords: Convection, diffusion, dissipation, finite volume method, method of lines.
1. Introduction
In this paper, we analyze a finite volume method for the following convection-diffusion-dissipation problem:
{ut(x,t)−(αux)x+v(x)ux+c(x)u=f(x,t) in QTu(x,0)=u0(x)u(0,t)=g(t)u(A,t)=h(t)\left\{\begin{aligned} u_{t}(x, t)-\left(\alpha u_{x}\right)_{x}+v(x) u_{x}+c(x) u=f(x, t) \text { in } Q_{T} \\ u(x, 0)=u^{0}(x) \\ u(0, t)=g(t) \\ u(A, t)=h(t) \end{aligned}\right.
where QT=]0,A[×]0,T[,AQ_{T}=] 0, A[\times] 0, T[, A and TT are two nonnegative integers, u(x,t)u(x, t) is the concentration of the pollutant, α\alpha is a constant, v(x)v(x) is the speed of water, c(x)c(x) is the dissipation, f(x,t)f(x, t) is the source function, xx is a space variable, and tt stands for time. utu_{t} denotes the time derivative of uu and uxu_{x} the derivative with respect to the space variable xx.
−(αux)x-\left(\alpha u_{x}\right)_{x} is the diffusion term, v(x)uxv(x) u_{x} is the convection term, and c(x)uc(x) u is the dissipation term.
Solving this problem has been recently addressed in [1], where the authors were able to calculate the exact solution by using the Adomian decomposition method on specific cases of α,v\alpha, v and cc. In general, it is not always possible or easy to calculate the exact solutions in many cases, this requires the use of a numerical method. The aim, here is to find an approximation to the solution, u(x,t)u(x, t) by a semi-discrete finite volume scheme.
The outline of the remainder of this paper is as follows. In next section, we present the semi-discrete finite volume scheme (FVM), for the problem (1). Section 3 is devoted to the fully discrete approximtion. Finally in Section 4, numerical results using this approach are shown and compared with the exact solution (Exact) obtained in [1].
2. Semi-discrete approximation
By semi-discretization we mean discretization only in space, not in time. This approach is also called method of lines (see, e.g.
[3]). We discretize space into NN equal size grid cells of size h=h= A/NA / N, and define xi=h/2+ihx_{i}=h / 2+i h, so that xix_{i} is the center of cell Ii=I_{i}= (xi−1/2,xi+1/2)\left(x_{i-1 / 2}, x_{i+1 / 2}\right).
In finite volume method the unknowns approximate the average of the solution over a grid cell. More precisely, we let ui(t)u_{i}(t) be the approximation
ui(t):=1h∫xi−1/2xi+1/2u(x,t)dxu_{i}(t):=\frac{1}{h} \int_{x_{i-1 / 2}}^{x_{i+1 / 2}} u(x, t) d x
Integrating (1) over the cell IiI_{i} and dividing par hh we get
1h∫xi−1/2xi+1/2ut(x,t)dx=αh∫xi−1/2xi+1/2(ux)x(x,t)dx−1h∫xi−1/2xi+1/2v(x)ux(x,t)dx−1h∫xi−1/2xi+1/2c(x)u(x,t)dx+1h∫xi−1/2xi+1/2f(x,t)dx∫xi−1/2xi+1/2(ux)x(x,t)dx=ux(xi+1/2,t)−ux(xi−1/2,t)\begin{aligned} & \frac{1}{h} \int_{x_{i-1 / 2}}^{x_{i+1 / 2}} u_{t}(x, t) d x=\frac{\alpha}{h} \int_{x_{i-1 / 2}}^{x_{i+1 / 2}}\left(u_{x}\right)_{x}(x, t) d x-\frac{1}{h} \int_{x_{i-1 / 2}}^{x_{i+1 / 2}} v(x) u_{x}(x, t) d x \\ & \quad-\frac{1}{h} \int_{x_{i-1 / 2}}^{x_{i+1 / 2}} c(x) u(x, t) d x+\frac{1}{h} \int_{x_{i-1 / 2}}^{x_{i+1 / 2}} f(x, t) d x \\ & \int_{x_{i-1 / 2}}^{x_{i+1 / 2}}\left(u_{x}\right)_{x}(x, t) d x=u_{x}\left(x_{i+1 / 2}, t\right)-u_{x}\left(x_{i-1 / 2}, t\right) \end{aligned}
Since the value in the midpoint of the cell is a second order approximation of average, we have for smooth uu,
ux(xi−1/2,t)=1h[u(xi,t)−u(xi−1,t)+O(h2)]u_{x}\left(x_{i-1 / 2}, t\right)=\frac{1}{h}\left[u\left(x_{i}, t\right)-u\left(x_{i-1}, t\right)+O\left(h^{2}\right)\right]
and
ux(xi+1/2,t)=1h[u(xi+1,t)−u(xi,t)+O(h2)]u_{x}\left(x_{i+1 / 2}, t\right)=\frac{1}{h}\left[u\left(x_{i+1}, t\right)-u\left(x_{i}, t\right)+O\left(h^{2}\right)\right].
To approximate the terms
∫xi−1/2xi+1/2v(x)ux(x,t)dx\int_{x_{i-1 / 2}}^{x_{i+1 / 2}} v(x) u_{x}(x, t) d x and ∫xi−1/2xi+1/2c(x)u(x,t)dx\int_{x_{i-1 / 2}}^{x_{i+1 / 2}} c(x) u(x, t) d x
we use the values of functions vv and cc in the midpoint of the cell, v(xi)v\left(x_{i}\right) et c(xi)c\left(x_{i}\right) respectively. So we get
∫xi−1/2xi+1/2c(x)u(x,t)dx≈hc(xi)u(xi,t)\int_{x_{i-1 / 2}}^{x_{i+1 / 2}} c(x) u(x, t) d x \approx h c\left(x_{i}\right) u\left(x_{i}, t\right) and
∫xi−1/2xi+1/2v(x)ux(x,t)dx≈v(xi)(u(xi+1/2)−u(xi−1/2))\int_{x_{i-1 / 2}}^{x_{i+1 / 2}} v(x) u_{x}(x, t) d x \approx v\left(x_{i}\right)\left(u\left(x_{i+1 / 2}\right)-u\left(x_{i-1 / 2}\right)\right)
Now, as in [2], we approximate v(xi)u(xi+1/2)v\left(x_{i}\right) u\left(x_{i+1 / 2}\right) by v(xi)u(xi)v\left(x_{i}\right) u\left(x_{i}\right) and v(xi)u(xi−1/2)v\left(x_{i}\right) u\left(x_{i-1 / 2}\right) by v(xi)u(xi−1)v\left(x_{i}\right) u\left(x_{i-1}\right). Fanally we get
dui(t)dt=αh2[(ui+1(t)−ui(t))−(ui(t)−ui−1(t))]−v(xi)h(ui(t)−ui−1(t))−c(xi)ui(t)+fi(t)\begin{aligned} & \frac{d u_{i}(t)}{d t}=\frac{\alpha}{h^{2}}\left[\left(u_{i+1}(t)-u_{i}(t)\right)-\left(u_{i}(t)-u_{i-1}(t)\right)\right] \\ & \quad-\frac{v\left(x_{i}\right)}{h}\left(u_{i}(t)-u_{i-1}(t)\right)-c\left(x_{i}\right) u_{i}(t)+f_{i}(t) \end{aligned}
So
dui(t)dt=(v(xi)h+αh2)ui−1−(2αh2+v(xi)h+c(xi))ui(t)+αh2ui+1(t)+fi(t)\begin{aligned} \frac{d u_{i}(t)}{d t}= & \left(\frac{v\left(x_{i}\right)}{h}+\frac{\alpha}{h^{2}}\right) u_{i-1}-\left(\frac{2 \alpha}{h^{2}}+\frac{v\left(x_{i}\right)}{h}+c\left(x_{i}\right)\right) u_{i}(t) \\ & +\frac{\alpha}{h^{2}} u_{i+1}(t)+f_{i}(t) \end{aligned}
for i=1,…N−2i=1, \ldots N-2. To complete the scheme (2) we need update formula also for the boundary points i=0i=0 and i=N−1i=N-1. These must be derived by taking the boundary conditions into account. We introduce the ghost cells I−1I_{-1} and INI_{N} with located juste outside the domain. The boundary conditions are used to fill these cells with values u−1u_{-1} and uNu_{N}, based on the values uiu_{i} in the interior cells. The same update formula (2) as before can then be used also for i=0i=0 and i=N−1i=N-1. Let us consider our boundary conditions u(0,t)=g(t)u(0, t)=g(t) and u(A,t)=h(t)u(A, t)=h(t). So u0(t)=g(t)u_{0}(t)=g(t) and uN(t)=h(t)u_{N}(t)=h(t). Since the center of cells I0I_{0} and INI_{N} are not on the boundary, we take the average of two cells to approximate the value in between,
g(t)=u(0,t)=u(0,t)+u(−1,t)2=u0(t)+u−1(t)2+O(h2)g(t)=u(0, t)=\frac{u(0, t)+u(-1, t)}{2}=\frac{u_{0}(t)+u_{-1}(t)}{2}+O\left(h^{2}\right)
and
h(t)=u(A,t)=u(N,t)+u(N−1,t)2=uN(t)+uN−1(t)2+O(h2)h(t)=u(A, t)=\frac{u(N, t)+u(N-1, t)}{2}=\frac{u_{N}(t)+u_{N-1}(t)}{2}+O\left(h^{2}\right)
leading to the approxiations
u−1(t)=2g(t)−u0(t)u_{-1}(t)=2 g(t)-u_{0}(t),
and
uN(t)=2h(t)−uN−1(t)u_{N}(t)=2 h(t)-u_{N-1}(t).
We now insert this into the update (2) for i=0i=0 and i=N−1i=N-1, we get
du0dt=αh2u1−(2αh2+v(x0)h+c(x0))u0+(v(x0)h+αh2)u−1+f0\frac{d u_{0}}{d t}=\frac{\alpha}{h^{2}} u_{1}-\left(\frac{2 \alpha}{h^{2}}+\frac{v\left(x_{0}\right)}{h}+c\left(x_{0}\right)\right) u_{0}+\left(\frac{v\left(x_{0}\right)}{h}+\frac{\alpha}{h^{2}}\right) u_{-1}+f_{0},
du0dt=−(3αh2+2v(x0)h+c(x0))u0+αh2u1+f0+2(v(x0)h+αh2)g(t)\frac{d u_{0}}{d t}=-\left(\frac{3 \alpha}{h^{2}}+\frac{2 v\left(x_{0}\right)}{h}+c\left(x_{0}\right)\right) u_{0}+\frac{\alpha}{h^{2}} u_{1}+f_{0}+2\left(\frac{v\left(x_{0}\right)}{h}+\frac{\alpha}{h^{2}}\right) g(t)
and
duN−1dt=αh2uN−(2αh2+v(xN−1)h+c(xN−1))uN−1+(v(xN−1)h+αh2)uN−2+fN−1duN−1dt=(v(xN−1)h+αh2)uN−2−(3αh2+v(xN−1)h+c(xN−1))uN−1+2αh2h(t)h2+fN−1\begin{aligned} \frac{d u_{N-1}}{d t}= & \frac{\alpha}{h^{2}} u_{N}-\left(\frac{2 \alpha}{h^{2}}+\frac{v\left(x_{N-1}\right)}{h}+c\left(x_{N-1}\right)\right) u_{N-1} \\ & +\left(\frac{v\left(x_{N-1}\right)}{h}+\frac{\alpha}{h^{2}}\right) u_{N-2}+f_{N-1} \\ \frac{d u_{N-1}}{d t}= & \left(\frac{v\left(x_{N-1}\right)}{h}+\frac{\alpha}{h^{2}}\right) u_{N-2}-\left(\frac{3 \alpha}{h^{2}}+\frac{v\left(x_{N-1}\right)}{h}+c\left(x_{N-1}\right)\right) u_{N-1} \\ & +\frac{2 \alpha}{h^{2}} \frac{h(t)}{h^{2}}+f_{N-1} \end{aligned}
We put all the equations (2), (3), (4) together and write them in matrix form, then we get a linear ordinary diffenrential equation (ODE) system of the form:
dU(t)dt=AU(t)+S(t)\frac{d U(t)}{d t}=A U(t)+S(t).
Here, in this semi-discretization the time-dependent partial differential equation has been approximated by a system of ODEs.
Remark
In a heterogeneous medium, α\alpha may be discontinuous, since the conductivities of differents components of the medium may be quite different. In this case the edges of cells have to coincide with the discontinuities of α\alpha, and the approximation the diffusion term is done by using the harmonic means of α\alpha values between the neighboring cells, as done in [2].
3. Fully discrete approximation
The system of ODEs (5) can be solve by standard numerical methods for ODEs with a time step Δt\Delta t e.g. the forward Euler method
Un+1=Un+Δt[AUn+S(tn)],Un≈U(tn),tn=nΔtU^{n+1}=U^{n}+\Delta t\left[A U^{n}+S\left(t_{n}\right)\right], U^{n} \approx U\left(t_{n}\right), t_{n}=n \Delta t.
That is the method we used in this paper, others methods can be use. As for any ODE method, here arises the question of stability. The numerical experiments that we will present in Section 4 have shown that it takes Δt≤Ch2\Delta t \leq C h^{2} as a condition of the method.
4. Numerical simulations
In this section we present numerical results, comparing the approximate solution described in this paper and the example of exact solution obtained in [1]. The first test problem involved simulations with following data:
QT=]0,10[×]0,T[,α=v(x)=c(x)=1,u0(x)=cos(x)Q_{T}=] 0,10[\times] 0, T[, \alpha=v(x)=c(x)=1, u^{0}(x)=\cos (x),
g(t)=cos(t),h(t)=cos(10)cos(t)g(t)=\cos (t), h(t)=\cos (10) \cos (t), and
f(x,t)=−cos(x)sin(t)+2cos(x)cos(t)−sin(x)cos(t)f(x, t)=-\cos (x) \sin (t)+2 \cos (x) \cos (t)-\sin (x) \cos (t).
The second test problem involved simulations with following data:
QT=]0,10[×]0,T[,α=0.5,v(x)=x(5−x)cos(x),c(x)=3x2Q_{T}=] 0,10[\times] 0, T[, \alpha=0.5, v(x)=x(5-x) \cos (x), c(x)=3 x^{2},
u0(x)=cos(x),g(t)=cos(t),h(t)=cos(10)cos(t)u^{0}(x)=\cos (x), g(t)=\cos (t), h(t)=\cos (10) \cos (t), and
f(x,t)=−cos(x)sin(t)+(0.5+3x3)cos(x)cos(t)−x(5−x)sin(x)cos(t)f(x, t)=-\cos (x) \sin (t)+\left(0.5+3 x^{3}\right) \cos (x) \cos (t)-x(5-x) \sin (x) \cos (t).
In both test cases the exact solution is (see, [1]):
u(x,t)=cos(x)cos(t)u(x, t)=\cos (x) \cos (t).
Figure 1: Test problem 1: h=0.1,T=1,Δt=0.001h=0.1, T=1, \Delta t=0.001
Figure 2: Test problem 1: h=0.2,T=10,Δt=0.0025h=0.2, T=10, \Delta t=0.0025
Figure 3: Test problem 1: h=0.1,T=10,Δt=0.0025h=0.1, T=10, \Delta t=0.0025
Figure 4: Test problem 2: h=0.2,T=1,Δt=0.002h=0.2, T=1, \Delta t=0.002
Figure 5: Test problem 2: h=0.1,T=1,Δt=0.001h=0.1, T=1, \Delta t=0.001
Figure 6: Test problem 2: h=0.1,T=15,Δt=0.001h=0.1, T=15, \Delta t=0.001
Acknowledgement
The author would like to thank the reviewers for their work.
5. Conclusion
The purpose of this paper was to apply and develop a finite volume scheme corresponding to the method of lines for the one-dimensional convection diffusion-dissipation equation. To approximate the diffusion and convection terms we used the values of the functions v(x)v(x) and c(x)c(x) in the midpoint of grid cells. If vv and cc are constant, one finds the approximations used in [2]. The numerical results indicate that the method of lines is well adapted to the discretization of this problem. The extension of the present technique to the two-dimensional problem with uniform rectangular grids is straightforward.
References
[1] J. Bonazebi Yindoula, P. Youssouf, G. Bissanga, F. Bassono, B. Some, “Application of the Adomian decomposition method and Laplace transform method to solving the convection diffusion-dissipation equation”, International Journal of Applied Mathematical Research, Vol.3, No.1, (2014), pp.30-35.
[2] R. Eymard, T. Gallouët, R. Herbin “The finite volume methods”, P.G. Curlet and J.L. Lions, editors, Techniques of Scientific Computing, Handbook of Numerical Analysis, VII, III, (2000), pp.713-1020.
[3] S. Hamdi, W. E. Schiesser, G. W, Griffiths, “Method of lines”, Scholarpedia, 2(7):2859, (2009).