(*	CR418		HW 7.1	Fogler Problem 4-32		        *)
(*	Unsteady Flow Problem - Startup of 3 Series CSTRs  	*)
(*  Rxn: A + B --> C, elementary.  Equimolar feed rate  *)
(*  Find the time required for the exit of the third 	*)
(*  reactor to reach 99% of its steady-state value		*)

(*  Important!  Start every problem by clearing			*)
(*	previous values										*)
Clear[dna1dt, dna2dt, dna3dt, fao, vdot, ca1, ca2, ca3,
		na1, na2, na3, na1r, na2r, na3r, v, fao, vdot,
		tfinal, k]

(* Mole balances for species A in each reactor. Writing	*)
(* these equations in the form:							*)
(* "accum. = in - out + generation" aids debugging		*)
dna1dt = fao - vdot ca1 - k ca1^2 v;
dna2dt = vdot ca1 - vdot ca2 - k ca2^2 v;
dna3dt = vdot ca2 - vdot ca3 - k ca3^2 v;

(* Auxiliary equations.  Relate all the variables in 	*)
(* problem to the dependent variables of the ODEs. Also	*)
(* assign values to the constants.  "ca3ss" is the		*)
(* steady-state value solved in a separate notbook.		*)
ca1 = na1[t]/v; ca2 = na2[t]/v; ca3 = na3[t]/v;
fao = 20;  v = 200; vdot = 20; tfinal = 65; k = 0.025;
ca3ss = 0.609;

(* NDSolve numerically evaluates the coupled ODEs given	*)
(* the initial conditions.								*)
sol = NDSolve[{	na1'[t] == dna1dt,
				na2'[t] == dna2dt,
				na3'[t] == dna3dt,
				na1[0] == 0,
				na2[0] == 0,
				na3[0] == 0},
				{na1, na2, na3},
				{t, 0, tfinal}];

(* The direct output of NDSolve is a List.  The set		*)
(* replace operator (/.) extracts the desired function	*)
(* from the List.									    *)
na1r = na1 /. sol[[1]];  na2r = na2 /. sol[[1]];
na3r = na3 /. sol[[1]];

Plot[{na1r[t]/v, na2r[t]/v, na3r[t]/v}, {t, 0, tfinal},
		GridLines -> Automatic, Frame -> True,
		FrameLabel -> {"Time (min)", "Conc. (mol/L)"}];

(* Find the point at which the exit conc. of Reactor 3 is	*)
(* 99% of its steady-state value.							*)
FindRoot[na3r[t]/v == 0.99 ca3ss, {t, 0}]

(*	Part d.	Split the feed of B		*)

Clear[dna1dt, dna2dt, dna3dt,ca1, ca2, ca3,
		na1, na2, na3, na1r, na2r, na3r,
		dnb1dt, dnb2dt, dnb3dt, nb1, nb2, nb3]

dna1dt = fao - vdot1 ca1 - k ca1 cb1 v;
dna2dt = vdot1 ca1 - vdot2 ca2 - k ca2 cb2 v;
dna3dt = vdot2 ca2 - vdot3 ca3 - k ca3 cb3 v;
dnb1dt = fbo1 - vdot1 cb1 - k ca1 cb1 v;
dnb2dt = vdot1 cb1 + fbo2 - vdot2 cb2 - k ca2 cb2 v;
dnb3dt = vdot2 cb2 - vdot3 cb3 - k ca2 cb3 v;

vdot1 = 15; vdot2 = vdot3 = 20; fbo1 = fbo2 = 10;
ca1 = na1[t]/v; ca2 = na2[t]/v; ca3 = na3[t]/v;
cb1 = nb1[t]/v; cb2 = nb2[t]/v; cb3 = nb3[t]/v;

sol2 = NDSolve[{	na1'[t] == dna1dt,
					na2'[t] == dna2dt,
					na3'[t] == dna3dt,
					nb1'[t] == dnb1dt,
					nb2'[t] == dnb2dt,
					nb3'[t] == dnb3dt,
					na1[0] == 0, na2[0] == 0, na3[0] == 0,
					nb1[0] == 0, nb2[0] == 0, nb3[0] == 0},
					{na1, na2, na3, nb1, nb2, nb3},
					{t, 0, tfinal}]

na1r = na1 /. sol2[[1]];  na2r = na2 /. sol2[[1]];
na3r = na3 /. sol2[[1]];  nb1r = nb1 /. sol2[[1]];
nb2r = nb2 /. sol2[[1]];  nb3r = nb3 /. sol2[[1]];

Plot[{na1r[t]/v, na2r[t]/v, na3r[t]/v}, {t, 0, tfinal},
		GridLines -> Automatic, Frame -> True, FrameLabel -> 
		{"Time (min)", "Conc. (mol/L)", "Species A"," "}]
Plot[{nb1r[t]/v, nb2r[t]/v, nb3r[t]/v}, {t, 0, tfinal},
		GridLines -> Automatic, Frame -> True, FrameLabel -> 
		{"Time (min)", "Conc. (mol/L)", "Species B", " "}]

FindRoot[na3r[t]/v == 0.99 ca3ss, {t, 0}]