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

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?

• Do you have quantitative boundary conditions we could play around with? Commented Apr 3, 2014 at 19:47
• I'm sorry, but I don't understand your question.
– user7847
Commented Apr 5, 2014 at 4:17

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).

• To add to Bill's answer, it is well known that finite difference or finite element schemes become unstable for advection dominated problems like yours. Search the literature for the term "transport stabilization" and you will find an infinite supply of schemes to remedy this, as well as a great deal of analysis. Commented Apr 3, 2014 at 13:03
• I thought that only finite-volume and finite-element had a Peclet number restriction for stability? Commented Apr 4, 2014 at 1:56
• Given that linear FEM and FVM methods can directly reduced to central-difference methods, you should expect that they should have similar $Pe$ restrictions. Commented Apr 4, 2014 at 3:42