Tank Reaction Equations

Special Mathematica Functions: flow2D[.] & vectorFlow2D[.]

Off[General::spell];
Off[General::spell1];
Off[InterpolatingFunction::dmval];
Off[ParametricPlot::pptr];
Off[Graphics::gprim];
Off[Graphics::gptn];
Off[NDSolve::ndsz];

Needs["Graphics`Colors`"]
	
colorQ[a_]:=(Head[a]===RGBColor);
grayQ[a_]:=(Head[a]===GrayLevel);
thickQ[a_]:=(Head[a]===Thickness);
absthickQ[a_]:=(Head[a]===AbsoluteThickness);
notruleQ[u_] := !(Head[u]===Rule);
ruleQ[u_] := (Head[u]===Rule);

Norm::usage = "Norm[expr] returns the Euclidean length of the
		vector, 'expr'.";
Norm[a_?VectorQ] := Sqrt[a.a];

tol = .5*10^-6;

Vector2D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics[Point[{0,0}]]];
	,						
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (.075 Norm[ntip-ntail]/Norm[ntip-ntail]) (ntip-ntail) ;
	up = 0.15 {-Punit[[2]],Punit[[1]]};
	arrow = Line[{ntail,ntip,ntip-Punit+up,ntip,ntip-Punit-up}];
	]];
Return[Graphics[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

		
uArrow[vec_,opts___]:= 
	Module[{a,b},
		If[Norm[N[vec]]<1,b = vec,b=vec/Norm[vec]];
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0}};
				Return[Vector2D[a,a+b,opts]]
				];

flow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
											
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

vectorFlow2D[{f_,g_},{x_,y_,t_},initVals_,tFinal_,dt_,opns___] :=
 Module[{soln,arr,pt,crv,tt},
soln = Flatten[Map[
	NDSolve[{x[0]==#[[1]],y[0]==#[[2]],
	x'[t] == f[x[t],y[t]],y'[t] == g[x[t],y[t]]},
	{x[t],y[t]},{t,0,tFinal}]&,
										initVals],1];
								
	arr[tt_] := Map[uArrow[{f[#[[1]],#[[2]]],g[#[[1]],#[[2]]]},Tail->#,Red]&,
	{x[t],y[t]}/.soln/.t->tt];
	
pt[tt_] := Map[Graphics[{Blue,AbsolutePointSize[4],Point[{#[[1]],#[[2]]}]}]&,
{x[t],y[t]}/.soln/.t->tt];
										
crv[tt_] := ParametricPlot[Evaluate[{x[t],y[t]}/.soln],{t,0,tt},opns,
			AspectRatio->1,DisplayFunction->Identity];
			
Show[{crv[dt/100],pt[0],arr[dt/100]},opns,AspectRatio->1,
DisplayFunction->$DisplayFunction];
Do[Show[{crv[tt],arr[tt],pt[tt]},DisplayFunction->$DisplayFunction];,
{tt,dt,tFinal,dt}];
	Return[soln];	
];

BeginPackage["Graphics`flow2D`"]
flow2D::usage =
"flow2D[{f,g},{x,y},initialValues,tFinal,dt,(options)]
	shows the flow solution of dx = f[x,y] dt, dy = g[x,y] dt 
	starting with the initial (x,y)-values = initialValues = {{a,b},{c,d},...}
	flowing up to time tFinal with animation frames every dt time units."

Begin["`Private`"]
stuff goes here 
End[]
EndPackage[]

Variables & Change Functions in General

t - time in hours
c - concentration of input chemical in pound-moles
T - temperature in tank in degrees Farenheit
r = T/Tin

\(Clear[f, g, c, r, p, a, b, h]; 

    f[c_, r_]\  := \ p*\((cIN - c)\)\  - \ c\ *Exp[a - b/r]; 

    g[c_, r_]\  := \ p*\((1 - r)\)\  + \ c*h*Exp[a - b/r]; 
f[c, r]\)g[c, r]

Linearization

Matrix of partials of the change functions

\(m[cE_, rE_]\  = 
      \ {{D[f[c, r], c], D[f[c, r], r]}, 

          		{D[g[c, r], c], D[g[c, r], r]}} /. {c -> cE, r -> rE}; 

    TableForm[m[cE, rE]]\)

Characteristic equation

Det[ m[cE,rE] - x Identity] = 0

\(TraditionalForm
      \`x\^2\  + \ β[cE, rE]\ x\  + \ γ[cE, rE]\  = \ 0\)

γ[c_, r_] = Det[m[c, r]]

β[c_, r_] = 
      -m[c, r][[1, 1]] - m[c, r][[2, 2]]

Both characteristic roots have negative real parts when both  β  and  γ  are positive.

A case study

Measured constants

\(q = 100; 
cIN = 1; 
Tin\  = \ 300; 
    \  (*250, \ fire\ goes\ out, \ 350\ everything\ used\ up*) 
V = \ 100; 

    A\  = \ E^25; 
B\  = \ 9000; 
H = 200; 
sC\  = \ 1; \)

Computed parameters

\(p\  = \ q/V; 
a\  = \ Log[A]; 
b\  = \ B/Tin; 

    h\  = \ H\ /\((sC*Tin)\); \)

Change Functions in Specific

f[c, r]g[c, r]

Equilibria

Solve[{f[c, r] == 0, g[c, r] == 0}, {c, r}]

verticalField = 
      ContourPlot[f[c, r], {c, 0, 3}, {r, 1, 2}, Contours -> {0}, 
        ContourShading -> False, ContourStyle -> {RGBColor[1, 0, 0]}];

horizontalField = 
      ContourPlot[g[c, r], {c, 0, 3}, {r, 1, 2}, Contours -> {0}, 
        ContourShading -> False, ContourStyle -> {RGBColor[0, 0, 1]}, 
        PlotPoints -> 60];

Show[verticalField, horizontalField];

\(verticalField = 
      ContourPlot[f[c, r], {c, 0, 0.01}, {r, 1.6, 1.7}, Contours -> {0}, 
        ContourShading -> False, ContourStyle -> {RGBColor[1, 0, 0]}]; 
    horizontalField = 
      ContourPlot[g[c, r], {c, 0,  .01}, {r, 1.6, 1.7}, Contours -> {0}, 
        ContourShading -> False, ContourStyle -> {RGBColor[0, 0, 1]}, 
        PlotPoints -> 50]; 
Show[verticalField, horizontalField]; \)

Accurate Root 1

roots = FindRoot[{f[c, r] == 0, g[c, r] == 0}, {c, 1}, {r, 1}]

Stability @ (cE,rE)

β[c, r] /. rootsγ[c, r] /. roots

Stable: β and γ  both positive.
Unstable: β and γ NOT both positive.

Accurate Root 2

roots = 
      FindRoot[{f[c, r] == 0, g[c, r] == 0}, {c, 0.001}, {r, 1.7}]

Stability @ (cE,rE)

β[c, r] /. rootsγ[c, r] /. roots

Stable: β and γ  both positive.
Unstable: β and γ NOT both positive.

Accurate Root 3

roots = FindRoot[{f[c, r] == 0, g[c, r] == 0}, {c, 0.8}, {r, 1.2}]

Stability @ (cE,rE)

β[c, r] /. rootsγ[c, r] /. roots

Stable: β and γ  both positive.
Unstable: β and γ NOT both positive.

The Flow of The Reaction

Be patient, this computation takes a long time, because there are a lot of initial conditions and a lot of time steps.

\(initconds = Table[{x, 1.13}, {x, 1, 1.1,  .1}]; 

    initconds\  = \ Union[initconds, Table[{x, 1.1}, {x,  .9, 2,  .1}]]; 

    initconds\  = \ Union[initconds, Table[{x, 1.17}, {x, 0.5,  .7,  .1}]]; 

    initconds\  = \ Union[initconds, Table[{x, 1.2}, {x, 0.1,  .3,  .1}]]; 

    initconds\  = \ Union[initconds, Table[{x, 1.1}, {x, 0.7,  .8,  .1}]]; 

    initconds\  = \ Union[initconds, Table[{0.05, x}, {x, 1.1, 1.3,  .05}]]; 
    
tfinal = 10; 
\(dt = 1; \)\)Print["Flow solution of the differential equations:"]; \(Print["dc = ", f[c, r], " dt"]; 

    Print["dr = ", g[c, r], " dt"]; 

    flow2D[{f, g}, {c, r, t}, initconds, tfinal, dt, 
      AxesLabel -> {"c", "r"}, 
      PlotRange -> {{0, 2}, {1, 2
			}}]; \)

Here's a closer look at the stable point:

\(f[c_, r_]\  := \ \((p*\((cIN - c)\)\  - \ c\ *Exp[a - b/r])\)/20; 

    g[c_, r_]\  := \((\ p*\((1 - r)\)\  + \ c*h*Exp[a - b/r])\); \)

\(initconds = Table[{x, 1.6}, {x,  .002,  .005,  .001}]; 

    initconds\  = \ 
      Union[initconds, Table[{x, 1.8}, {x,  .002,  .005,  .001}]]; 

    initconds\  = \ 
      Union[initconds, Table[{ .0001, x}, {x, 1.6, 1.8,  .05}]]; 

    initconds\  = \ Union[initconds, Table[{ .005, x}, {x, 1.6, 1.8,  .05}]]; 
    
tfinal = 10; 
dt = 1; 

    \(Print["Flow solution of the differential equations:"]; \)\)\(Print["dc = ", f[c, r], " dt"]; 

    Print["dr = ", g[c, r], " dt"]; 

    flow2D[{f, g}, {c, r, t}, initconds, tfinal, dt, 
      AxesLabel -> {"c", "r"}, 
      PlotRange -> {{0,  .005}, {1.6, 1.8
			}}]; \)