3
$\begingroup$

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.

$\endgroup$
  • $\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 '14 at 17:25
4
$\begingroup$

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.

$\endgroup$
  • $\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 '14 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$ – Geoff Oxberry Jan 21 '14 at 19:45
  • $\begingroup$ Thanks again! This looks promising. I use C, Fortran and Python, so IDA looks just fine. $\endgroup$ – seb Jan 22 '14 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 '14 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$ – Geoff Oxberry Jan 31 '14 at 19:09

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.