2 A Test Problem, the Convection Diffusion Equation
As a test case, we consider the one dimensional, steady state convection-diffusion equation with fixed boundary conditions
\[ u_x = bu_{xx} + 1 , \quad u(0) = u(1) = 0. \tag{2.1}\]
Here \(b\) is some physical parameter. Moreover, \(u(x)\) is defined on the interval \([0,1]\). This equation has a solution that is given by
\[ u(x) = x - \frac{e^{-(1-x)/b} - e^{-1/b}}{1-e^{-1/b}}. \tag{2.2}\]
We are however interested in solving this numerically, with a finite difference approach. We partition the interval \([0,1]\) using \(n+2\) equidistant points \(x_i, i = 0, \dots n+1\). We denote the distance between each points as \(\Delta x = \frac{1}{n+1}\). The approximated value of \(u\) at the point \(x_i\) is denoted by \(u^i\) and we have \(u^0 = u(0) = 0\) and \(u^{n+1} = u(1) = 0\). We approximate, for \(i \geq 1\), the derivative as
\[ u_x^i = \frac{u^i - u^{i-1}}{\Delta x}, \] and the second order derivative is approximated by \[ u^i_{xx} = \frac{u^{i+1} - 2u^i + u^{i-1}}{\Delta x ^2}. \]
Note that the first derivative is approximated backward in space, which aligns with the convection being in the right direction. For \(i = 1 , \dots , n\), we thus have the approximation
\[ \frac{u^i - u^{i-1}}{\Delta x} = b \frac{u^{i+1} - 2u^i + u^{i-1}}{\Delta x ^2} + 1. \]
This can be given in a matrix-vector format, by letting \(\mathbfit{u} = (u^1,\dots, u^{n})^\intercal\),
\[ \mathbfit{Au} = \mathbfit{Bu}+\mathbfit{d}, \]
where \(\mathbfit{d} = (1,1,\dots , 1)^\intercal\),
\[ \mathbfit{A} = \frac{1}{\Delta x}\begin{pmatrix} 1 & &&&\\ -1 & 1 &&\\ & -1 & 1 &\\ &&\ddots & \ddots &\\ &&&-1& 1 \end{pmatrix}, \]
and
\[ \mathbfit{B} = \frac{b}{\Delta x ^2}\begin{pmatrix} -2 & 1 &&&\\ 1 & \ddots & \ddots &&\\ & \ddots & \ddots & \ddots&\\ && \ddots & \ddots & 1 &\\ &&&1 & -2 \end{pmatrix} . \]
Note that from now on, matrices and vector are denoted in bold italic. With \(\mathbfit{N} = \mathbfit{A}-\mathbfit{B}\), the approximate solution of Equation eq-convec_diff_steady is then the solution of the linear system
\[ \mathbfit{Nu} = \mathbfit{d}, \tag{2.3}\]
where \(\mathbfit{N}\) is a square matrix of dimension \(n \times n\) and \(\mathbfit{d}\) is the vector of ones of dimension \(n\).
Remark. \(\mathbfit{N}\) is diagonally dominant. Since all elements of the diagonal are positive, we can use Gershgorin circle theorem to prove that all the eigenvalues of \(\mathbfit{N}\) have a positive real part. We thus only need to assume \(\mathbfit{N}\) is non singular to prove that \(-\mathbfit{N}\) is stable.
We plot two examples of what the exact solution (Equation eq-convec_diff_th_solution) and the discretized solution (Equation eq-UnscaledSystem) look like for different values of \(b\) in Figure fig-th_vs_dis
To solve this linear system, we use the method highlighted before. To make it easier for later, we choose to scale \(N\) so that its diagonal elements are \(1\). This allows us to have all eigenvalues in the circle centered around \(1\) with radius \(1\) independently of the parametrization. Setting \(\eta = \frac{1}{\Delta x} + \frac{2b}{\Delta x^2}\), solving Equation eq-UnscaledSystem is equivalent to solving the system
\[ \mathbfit{M} \mathbfit{u} = \mathbfit{e}, \tag{2.4}\]
where with \(\mathbfit{M} = \frac{\mathbfit{N}}{\eta}\), \(\mathbfit{e} = \frac{\mathbfit{d}}{\eta}\). The eigenvalues of \(\mathbfit{N}\) are also scaled by \(\frac{1}{\eta}\), and therefore \(-\mathbfit{M}\) is stable, assuming it is non singular. We are now ready to solve the system iteratively using an ODE solver. To do that, we introduce a (pseudo) time variable \(t\) and we consider the ODE
\[ \mathbfit{u'(t)} = \mathbfit{e} - \mathbfit{Mu(t)}. \tag{2.5}\]
where \(\mathbfit{M}\) and \(\mathbfit{e}\) depends on both \(n\) and \(b\). From now on, we call \(b\) and \(n\) the problem parameters. We can use Theorem thm-steadyState with the non singularity assumption to guarantee that \(\mathbfit{u(t)}\) converges to a steady state independently of its chosen initial value. In the next chapter, we introduce a numerical method to solve this differential equation, which we will use in this thesis.
Remark. The convection diffusion equation is derived in [1] , chap. 3 from the continuity equation for a scalar quantity \(u\) and is
\[ \frac{\partial u}{\partial t} + \Delta.(\vec{v}u - \nabla(Du)) = R. \]
We will assume that the quantity \(u\), is the temperature in Kelvin, and has the S.I unit \(K\). The physical quantities are:
- \(\vec{v}\), which is the velocity of the medium the quantity is in, in \(\text{ms}^{-1}\). (the advection/convection).
- \(D\) is the diffusion coefficient, in \(\text{m}^2\text{s}^{-1}\).
- \(R\) is governing whether the quantity is created when \(R>0\), or destructed when \(R<0\). The unit is \(K\text{s}^{-1}\).
We can now simplify the equation by considering it in a single dimension \(x\)
\[ \frac{\partial u}{\partial t} + \frac{\partial}{\partial x}.(vu - \frac{\partial}{\partial x}(Du)) = R. \]
This further simplifies to \(u_t + vu_x - Du_{xx} = R\).
Then, in the steady state, \(u_t = 0\) so we get
\[ u_x = \frac{D}{v}u_{xx} + \frac{R}{v}, \]
and we recognize Equation eq-convec_diff_steady, with \(b = \frac{D}{v}\), and \(1 = \frac{R}{v}\). This also means that we “lose” two parameters in the studied test problem for simplification purposes. Nevertheless, this can be used to give some degree of intuition behind Figure fig-th_vs_dis. When the diffusion is high compared to the convection, the quantity is more centered, but on the other hand, when the convection speed is high compared to the diffusion, the quantity \(u\) is “flushed” to the right (assuming \(v>0\)).