# Problem in Discretizing Convection-Diffusion-Reaction equation

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.

• 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? – seb Jan 20 '14 at 17:25

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}

solve

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

• 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). – seb Jan 21 '14 at 11:44
• 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? – Geoff Oxberry Jan 21 '14 at 19:45
• Thanks again! This looks promising. I use C, Fortran and Python, so IDA looks just fine. – seb Jan 22 '14 at 8:36
• @GeoffOxberry, by what definition of DAE would you describe your 2nd equation as a DAE but not the 1st? – Aurelius Jan 31 '14 at 18:39
• @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. – Geoff Oxberry Jan 31 '14 at 19:09