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?

  • $\begingroup$ Do you have quantitative boundary conditions we could play around with? $\endgroup$ – dearN Apr 3 '14 at 19:47
  • $\begingroup$ I'm sorry, but I don't understand your question. $\endgroup$ – watanabenoburu Apr 5 '14 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).

  • 2
    $\begingroup$ 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. $\endgroup$ – Wolfgang Bangerth Apr 3 '14 at 13:03
  • $\begingroup$ I thought that only finite-volume and finite-element had a Peclet number restriction for stability? $\endgroup$ – boyfarrell Apr 4 '14 at 1:56
  • 1
    $\begingroup$ Given that linear FEM and FVM methods can directly reduced to central-difference methods, you should expect that they should have similar $Pe$ restrictions. $\endgroup$ – Bill Barth Apr 4 '14 at 3:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.