Solving the convection-diffusion equation using finite differences at high Peclet numbers

Question

I am trying to solve the 2D convection-diffusion equation, which in non-dimensional variables can be expressed for my problem as:

$\frac{\partial c}{\partial t} = \lambda^2\frac{\partial^2c}{\partial x^2} + \frac{\partial^2c}{\partial y^2} - 6\, Pe\, y (1 - y)\frac{\partial c}{\partial x}$

Here, $\lambda$ is a geometrical parameter (an aspect ratio of sorts) and $Pe$ is the Peclet number. The $y\,(1 - y)$ appears due to a parabolic velocity profile at the inlet (left wall). The boundary conditions are:

• Left wall: Dirichlet type
• Top and Right walls: No flux, so $y$ and $x$ gradients are zero respectively
• Bottom wall: Robin type boundary conditions. Basically a surface reaction that follows Langmuir kinetics.

I have discretized the equation using central differences for space and am using the Crank-Nicholson scheme to implicitly march in time. My code blows up (extremely large values of c, ~ $10^{134}$) for $Pe > 100$. It works fine when $Pe$ is below that.

Why does this happen? What do I do?

Answer 1

To quote Gresho: "The wiggles are telling you something." Namely, you need more mesh.

Without seeing your solution, or knowing $\lambda$, that's my best guess. I suspect that if you look at your first time step for $Pe$ greater than but near 100, the solution hasn't completely blown up, but it does appear oscillatory. If that's the case, you either need to stabilize, or you need more mesh in the boundary layer(s).