Impose Neumann Boundary Condition in advection-diffusion equation 1D

Question

when solving the advection equation in 1D that is:

$$ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0 $$ with $ u'(t,0) = 0$ and $u(t,L) = 0$ , $u(0,x) = u_{0} $

one numerical scheme is the FTCS (Forward time-centered space), but this numerical scheme is unstable.

$$\frac{u_{j}^{n+1}-u_{j}^{n}}{ h_{t}} = c \frac{u_{j+1}^{n}-u_{j-1}^{n}}{ 2h_{x}} $$

But when solving

$$ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = \alpha \frac{\partial^2 u}{\partial x^2} $$ the advection-diffusion equation in 1D with $ u'(t,0) = 0$ and $u(t,L) = 0$ , $u(0,x) = u_{0} $

Since the advection-diffusion equation is a second order equation I'd like to use a second order approximation.

if we define $u_{k}^{n} := u(t_{n},x_{k}) $; $\ \ x_{k} = kh $ and $ \ \ k = 0,1,2,...,N$. $h$ is known as the mesh size or step size.

For the second derivative:

$$ \frac{\partial^2u_{k}^{n}}{\partial x^2} \approx \frac{u_{k+1}^{n}-2u_{k}^{n}+u_{k-1}^{n}}{h^2} = \frac{u_{k-1}^{n}-2u_{k}^{n}+u_{k+1}^{n}}{h^2} $$ for $k=0,1,...,N-1$

Since $u'(t,0) = 0$ and $ u(t_{n},L) = u(t_{n},x_{N}) = u_{N}^{n} = 0 $ we get the following matrix representation of the second derivative operator

\begin{equation} \frac{\partial^2}{\partial x^2} \approx L_{2} = \frac{1}{h^2}\left(\begin{matrix} -2 & 1 & & 0\\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ 0 & & 1 & -2 \end{matrix} \right) \end{equation}

for $k=0$ , we get

$$ \frac{\partial u_{0}^{n}}{\partial x} = \frac{u_{0+1}^{n}-u_{0-1}^{n}}{ 2h} = 0 $$ this implies that $ u_{1}^{n}=u_{-1}^{n} $ and

$$ \frac{\partial^2u_{0}^{n}}{\partial x^2} \approx \frac{u_{0+1}^{n}-2u_{0}^{n}+u_{0-1}^{n}}{h^2} = \frac{u_{0-1}^{n}-2u_{0}^{n}+u_{0+1}^{n}}{h^2} = \frac{-2u_{0}^{n}+2u_{1}^{n}}{h^2} $$

thus we have to modify the entry $1,2$ of $L_{2}$

\begin{equation} L_{2} = \frac{1}{h^2}\left(\begin{matrix} -2 & 2 & & 0\\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ 0 & & 1 & -2 \end{matrix} \right) \end{equation}

What I have done, is $\mathbf{impose \ the \ Neumann \ boundary \ condition}$ in $L_{2}$ .

I want to approximate the first derivative using central difference(Second order approximation):

$$ \frac{\partial u_{k}^{n}}{\partial x} = \frac{u_{k+1}^{n}-u_{k-1}^{n}}{ 2h} $$

The matrix representation is: \begin{equation} \frac{\partial}{\partial x} \approx L_{1} = \frac{1}{2h}\left(\begin{matrix} 0 & 1 & & 0\\ -1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ 0 & & -1 & 0 \end{matrix} \right) \end{equation}

for $k=0$ , we get

$$ \frac{\partial u_{0}^{n}}{\partial x} = \frac{u_{0+1}^{n}-u_{0-1}^{n}}{ 2h} = 0 $$ this implies that $ u_{1}^{n}=u_{-1}^{n} $

But I'm stuck when I try to $\mathbf{impose \ the \ neumann \ boundary \ condition \ in}$ $L_{1}$. I don't know how to do that.

If we solve that problem, we can solve the differential equation

$$ \frac{ \partial u }{\partial t} = \Big(-cL_{1} +\alpha L_{2} \Big)u $$ integrating in time( C-N, Back-Euler, RK4 )

$\mathbf{please \ help!\ \ How \ do \ you \ impose \ the \ Neumann \ Boundary \ Condition \ in \ L_{1}?}$

Answer 1

First, write out the semi-discrete equation for $u_0$, assuming it's an interior node:

$ \frac{d}{dt}(u_0) = - c \frac{\left( u_1 - u_{-1} \right)}{2h} + \alpha \frac{\left(u_1 - 2u_0 + u_{-1}\right)}{h^2}$

Then, eliminate the ghost node $u_{-1}$ by using the equation you obtained from the Neumann BC: $u_{-1} = u_1$. Substituting this condition into the above equation gives:

$ \frac{d}{dt}(u_0) = \alpha \frac{\left(- 2u_0 + 2u_1 \right)}{h^2}$

Note that the advection term vanishes, so the first row of your $L_1$ matrix should be all zeros. The first row of the $L_2$ matrix should have coefficients of $-2$ and $2$ in the first two columns, as you have already derived.

Also something to keep in mind, using a central difference approximation for the first derivative in the advection term will cause stability issues if your problem is advection dominant. In particular, this central difference scheme may develop oscillations if the grid Peclet number is greater than $2$, where the Peclet number is defined as:

$Pe = \frac{c h}{\alpha}$.

You can avoid these oscillations by using an upwinded (biased) approximation, such as this first-order upwinded scheme:

$\frac{\partial u_k}{\partial x} \approx \frac{u_k - u_{k-1}}{h}$, if $c > 0$.

Higher-order upwinded approximations also do exist, and typically involve wider stencils.