I'm trying to solve the Convection-Diffusion-Reaction (CDR) equation on a rectangular domain, using cylindrical coordinates and Finite Difference Methods (FDM) (this approximates a flow reactor).

Specifically, I want to solve

$D ( \frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} ) - v(r) \frac{\partial u}{\partial z} = S(r)$

Here, r and z are the coordinates in the respective directions, u is the property I want to solve for, v is the laminar velocity profile, S is a given source and D is the diffusion coefficient.

The Boundary conditions are:

$\frac{\partial c}{\partial r} = 0, $ at $r=0$ and at $r=R$ (where R is the length of the rectangle in the r direction).

My approach is to use the Method of Lines, i.e. to solve the system in the r direction for each step in the z direction.

When discretizing the above equation, I run into the practical problem that I divide by the velocity v, i.e. my discretized version of the above equation reads:

$\frac{du}{dz} = \frac{D}{v_i}(\frac{u_{i+1} - 2u_i + u{i-1}}{2h}+\frac{1}{r_i}\frac{u_{i+1}-u_{i-1}}{2h}) - S_i$

The problem with this is that $v_i=0$ at one edge of the computational domain, i.e. where $r=R$

So far, I have approximated my Boundary Conditions by using central differences. This means that I divide entries in my solution matrix by 0.

I'm thinking of fixing this problem with using either unsymmetric difference approximations.

I just wanna ask if that is a good idea or if there is a better fix to this problem.

  • $\begingroup$ upon googling this and looking a bit closer, I think the simplest solution is not to use the analytical laminar flow profile here. In most of the plots that I have seen from some CFD slides, the velocity is never exactly 0 on the boundary. I guess this is because you have to use Neumann Boundary Conditions to find the velocity profile with Navier-Stokes-Eqs as well, i.e. you are always interpolating through the boundary one way or the other. Is this impression going in the right direction? $\endgroup$
    – seb
    Jan 20, 2014 at 17:25

1 Answer 1


The simplest fix here is to treat the problem as a differential-algebraic equation, which is consistent with many treatments of PDE discretizations. Instead of solving

\begin{align} \frac{du}{dz}(r_{i}) = \frac{D}{v_i}(\frac{u_{i+1} - 2u_i + u{i-1}}{2h}+\frac{1}{r_i}\frac{u_{i+1}-u_{i-1}}{2h}) - \frac{S_i}{v_{i}} \end{align}


\begin{align} v_{i}\frac{du}{dz}(r_{i}) = D(\frac{u_{i+1} - 2u_i + u{i-1}}{2h}+\frac{1}{r_i}\frac{u_{i+1}-u_{i-1}}{2h}) - S_i, \end{align}

which can be done using standard DAE solvers.

  • $\begingroup$ Thank you! I haven't used DAE-solvers nor DAEs before. If possible, could you direct me to a solver and or a getting started place (google is not of much help). $\endgroup$
    – seb
    Jan 21, 2014 at 11:44
  • $\begingroup$ I tend to recommend IDA, which is part of the SUNDIALS suite of solvers; DASSL and DASPK are earlier DAE solvers related to IDA. There are also DAE solvers in PETSc, and several Fortran-based DAE solvers that have been developed over the years (DASSL and DASPK are two such solvers; there are others). Are you looking for solvers in a particular language? $\endgroup$ Jan 21, 2014 at 19:45
  • $\begingroup$ Thanks again! This looks promising. I use C, Fortran and Python, so IDA looks just fine. $\endgroup$
    – seb
    Jan 22, 2014 at 8:36
  • $\begingroup$ @GeoffOxberry, by what definition of DAE would you describe your 2nd equation as a DAE but not the 1st? $\endgroup$
    – Aurelius
    Jan 31, 2014 at 18:39
  • $\begingroup$ @Aurelius: $v_{i} = 0$ at the centerline, so at the centerline, the second equation is algebraic (its left-hand side is zero) and the first equation is undefined due to division by zero. $\endgroup$ Jan 31, 2014 at 19:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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