Quasilinear Solution Sequence for Nonlinear
Tubular Reactor with Dispersion
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :

Governing equations can be found in:

This problem was presented in the excellent book by W. F. Ramirez, Computational Methods for Process Simulation, Second Edition, 1997, Butterworth Heinemann, Oxford.

Off[NDSolve::"bvluc"]

Attempting to solve the nonlinear ODE with NDSolve fails because it is a nonlinear ODE and a split bounadary problem.

Pe = 10; R = 0.5; C0 = 5;

NDSolve[{D[c[z], {z, 2}] ==  -Pe D[c[z], {z, 1}] + 
            Pe R c[z]^2,  
        C0 ==  
          c[1] + 1/Pe (D[c[z], {z, 1}] /. z -> 1), (D[
                c[z], {z, 1}] ==  0 /. z ->  0)}, {c[z]}, {z, 
        0, 1}]

We assume a starting solution as constant 
concentration = 5. We then proceed with the quasilinear sequence until we get convergence to the solution of the original problem. Here four iterations are sufficient.

C0 = 5;

sol1 = 
      NDSolve[{D[c[z], {z, 2}] ==  -Pe D[c[z], {z, 1}] + 
              2 Pe R c[z] C0 - Pe R C0^2,  
          C0 ==  
            c[1] + 1/Pe c'[1], (D[c[z], {z, 1}] ==  0 /. 
              z ->  0)}, {c[z]}, {z, 0, 1}]

C0 = First[c[z] /. sol1]

sol2 = 
      NDSolve[{D[c[z], {z, 2}] ==  -Pe D[c[z], {z, 1}] + 
              2 Pe R c[z] C0 - Pe R C0^2,  
          5 ==  
            c[1] + 1/Pe c'[1], (D[c[z], {z, 1}] ==  0 /. 
              z ->  0)}, {c[z]}, {z, 0, 1}]

C0 = First[c[z] /. sol2]

sol3 = 
      NDSolve[{D[c[z], {z, 2}] ==  -Pe D[c[z], {z, 1}] + 
              2 Pe R c[z] C0 - Pe R C0^2,  
          5 ==  
            c[1] + 1/Pe c'[1], (D[c[z], {z, 1}] ==  0 /. 
              z ->  0)}, {c[z]}, {z, 0, 1}]

C0 = First[c[z] /. sol3]

sol4 = 
      NDSolve[{D[c[z], {z, 2}] ==  -Pe D[c[z], {z, 1}] + 
              2 Pe R c[z] C0 - Pe R C0^2,  
          5 ==  
            c[1] + 1/Pe c'[1], (D[c[z], {z, 1}] ==  0 /. 
              z ->  0)}, {c[z]}, {z, 0, 1}]

Plot[{c[z] /. sol1, c[z] /. sol2, c[z] /. sol3, c[z] /. sol4}, {z, 0, 
        1}, PlotStyle ->  {RGBColor[1, 0, 1], RGBColor[1, 0, 0], 
          RGBColor[0, 1, 0], RGBColor[0, 0, 1]}]