# 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? – dearN Apr 3 '14 at 19:47
• I'm sorry, but I don't understand your question. – watanabenoburu Apr 5 '14 at 4:17

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).
• Given that linear FEM and FVM methods can directly reduced to central-difference methods, you should expect that they should have similar $Pe$ restrictions. – Bill Barth Apr 4 '14 at 3:42