Solution to a System of Ordinary Differential Equations with 
Mixed, Algebraic, Boundary Conditions

by:  Richard C. Miller, Graduate Student
Department of Chemical Engineering
Michigan Technological University, Houghton, MI  49931
e-mail:  rcmiller@mtu.edu
voice: (906) 487-3145

Abstract

This notebook illustrates the solution of four (4) simultaneous ordinary differential equations with three (3) constant initial boundary conditions and one (1) algebraic final boundary condition. Forward shooting was used to solve the linear system.  The example provided illustrates the temperature profile of an infinitely long cylindrical reactor where the internal reaction is temperature dependent.

Summary

Physical Meaning

The differential equation shown below describes the temperature profile of an infinitely long cylindrical reactor where the internal chemical reaction is temperature dependent.  Heat is transferred away from the reactants through the verticle wall at the outer boundary.  Figure 1 shows the resulting solution to the steady state temperature profile at an arbitrary Z-axis cross section (i.e. infinite length).  The solution predicts that the temperature will be the greatest in the center and will decrease at radial distances closer to the wall.

Mathematical Technique

The solution to the above differential equation can be reduced to solving

Four (4) Simultaneous Ordinary Differential Equations with
•    Three (3) Initial Boundary Conditions
•     One (1) Final Boundary Condition

by  substitution with

Yielding:

These ODE's can be solved using a Forward Shooting Method and standard Mathematica  Built-in Objects.  

The procedure is as follows:

1.    Estimate the Single Unknown Initial Boundary Condition (IBguess).
             
             2.    Obtain a Solution Around this Initial BC and Inspect the
                    Resulting Final Boundary Condition.
             
             3.    If Necessary, Iterate on steps 1. and 2. until the Final BC is Satisfied.

Implementation

Simultaneous Ordinary Differential Equations to Solve:

ODE1 = (y[1]'[r] == y[2][r]);

ODE2 = (y[2]'[r] == If[r == 0, 		-G/2*y[1][r]^n, 
	     			If[0 < r <= 1,	-1/r*y[2][r] - G*y[1][r]^n]]);
		 		 		   
ODE3 = (y[3]'[r] == y[4][r]);
		 				  
ODE4 = (y[4]'[r] == If[r == 0, 		-n*G/2*y[1][r]^(n-1)*y[3][r],
		 			If[0 < r <= 1,	-1/r*y[4][r] - n*G*y[1][r]^(n-1)*y[3][r]]]);

Initial Boundary Conditions:

IBC1[IBguess_] := (y[1][0] == IBguess);

IBC2 = (y[2][0] == 0);

IBC3 = (y[3][0] == 1);

IBC4 = (y[4][0] == 0);

Final Boundary Condition to Satisfy:

FBC1 = (bign/2*y[1][1][IBguess] + y[2][1][IBguess] - bign/2 == 0);

Numerical Assignments for Physical Situation:

G     = 0.05;
n     = 1.5;
bign  = 1;
rmin = 0;
rmax  = 1;

Defining the Solution to the System of ODE's at IBguess:

Clear[sol];
sol[IBguess_] := NDSolve[{ODE1,ODE2,ODE3,ODE4,

							IBC1[IBguess],IBC2,IBC3,IBC4},
		 
		 {y[1][r],y[2][r],y[3][r],y[4][r]},
			  
		 {r,rmin,rmax}]

Clear[y];
y[n_][m_][IBguess_] := First[(y[n][r] /. Evaluate[sol[IBguess]]) /. r->m]

Iterating to Find the Necessary IBguess that will Satisfy the Final Boundary Condition:

FindRoot[ (bign/2*y[1][1][IBguess] + y[2][1][IBguess] - bign/2 == 0),

		 {IBguess, {1.5,1.3}}]

Plotting the Results:

answer = Evaluate[sol[1.06842]]

Clear[temperature];
temperature[radius_] := First[(y[1][r] /. answer) /. r -> radius]

Clear[tempcartesian];
tempcartesian[x_,y_] := Which[	N[Sqrt[x^2+y^2]] <= 1,
								    temperature[N[Sqrt[x^2+y^2]]],
							    N[Sqrt[x^2+y^2]] >= 1,
							    	temperature[N[1]]]

tempprofile = ContourPlot[tempcartesian[x,y],{x,0,1},{y,0,1},
						AxesLabel -> {"r/Radius",None},
						PlotLabel -> FontForm["Temp. Profile",{"Helvetial-Bold",12}],
						Ticks -> {{0.0,0.5,1.0},{0.0,0.5,1.0}},
						Frame -> False,
						Axes -> True]

Acknowledgments

Author would like to acknowledge A. Pintar, Professor at Michigan Technological University, for introducing me to the forward shooting method via IMSL.

References

1).  Davis, Mark, E., Numerical Methods and Modeling for Chemical Engineers,  John Wiley and Sons, New York, 1984.

Appendix; Cells to Create the Graphics in the Notebook Summary

Reactor

reactor = ParametricPlot3D[{Sin[t], Cos[t], z},
									
												{t,0,2Pi},{z,-2,2},
												
												Boxed -> False,
												AxesEdge -> {None,None,Automatic},
												AxesLabel -> {None,None,"Z-Axis"},
												PlotLabel -> FontForm["Reactor",{"Helvetial-Bold",12}],
												Prolog -> {Text["To - Infinity", {0,0},{-1,0}],
																					  Text["To + Infinity", {0,1},{-1,0}]}]

Combining Results

Show[

	Graphics[Rectangle[{0,0},{.7,1},reactor]], 
	Graphics[Rectangle[{.6,.6},{1.0,1.2},tempprofile]],
	
	PlotLabel -> 
				FontForm["Figure 1: Temperature Profile in an Infinitely Long Cylindrical Reactor",{"Helvetica-Bold",12}]]