(*  CR418 10/12/94						*)
(*  Homework #6							*)
(*	Multiple reactions, reversible		*)
(* 	reactions, and phase change			*)

(*  Reaction Network:                   *)
(*  CO + 2H2 <--> CH3OH                 *)
(*  CO + H2O <--> CO2 + H2O             *)
(*  CH3OH --> CH2O + H2                 *)


(* Always clear all variables first     *)

Clear [ca, cb, cc, cd, ce, cf, r, t, p,
		r1, r2, r3, a1, a2, a3, ea1, ea2, ea3,
		k1, k2, k3, k1eq, k2eq, fa, fb, fc, fd,
		fe, ff, yao, ybo, ydo, yinerto, fao, fbo, 
		fdo, finert, vdot, vdoto, v, ftot,
		far, fbr, fcr, fdr, fer, ffr, sumfi,
		faa, fba, fca, fda, fea, ffa]

(* Write the reaction rates and Arrhenius *)
(* expressions                            *)

r1 = k1(ca cb^2 - (1/k1eq) cc);
r2 = k2(ca cd - (1/k2eq) ce ca);
r3 = k3 cc;

k1 = a1 Exp[-ea1/(r t)];
k2 = a2 Exp[-ea2/(r t)];
k3 = a3 Exp[-ea3/(r t)];
k1eq = 8.642 10^-3 Exp[20.698 10^3/(r t)];
k2eq = 2.24 10^-2  Exp[8.874 10^3/(r t)];

(* Relate concentrations to the molar flow *)
(* rates - the dependent variables         *)

ca = fa[v]/vdot; cb = fb[v]/vdot; cc = fc[v]/vdot;
cd = fd[v]/vdot; ce = fe[v]/vdot; cf = ff[v]/vdot;
vdot = vdoto (sumfi/ftoto); 

(* Fill in the constants.  Time is varied *)
(* iteratively.  This would be more cleanly *)
(* done in a loop.                         *)

a1 = 2 10^11; ea1 = 1 10^4;
a2 =  1.2 10^9  ;  ea2 = 15 10^3;
a3 =   5 10^16 ; ea3 = 35 10^3;
r = 1.987;  p = 50; vfinal = 200;

t = 368;

ftoto = 1 10^-3;  (*  mol/s  *)
vdoto = ftoto/p (82.1 t);
yao = 0.15; ybo = 0.25; ydo = 0.25; yeo = 0.10; yinerto = 0.25;
fao = ftoto yao; fbo = ftoto ybo; fdo = ftoto ydo; 
feo = yeo ftoto; finert = ftoto yinerto;

sumfi = fa[v] + fb[v] + fc[v] + fd[v] + fe[v] + ff[v] + finert;

(* Solve the numerical ODEs                 *)

sol =	NDSolve[{	fa'[v] == -r1	-r2,
					fb'[v] == -2r1	+r2		+r3,
					fc'[v] == +r1			-r3,
					fd'[v] == 		-r2,
					fe'[v] ==		+r2,
					ff'[v] ==				+r3,
					fa[0] == fao, fb[0] == fbo,
					fc[0] == 0, fd[0] == fdo,
					fe[0] == feo, ff[0] == 0},
					{fa, fb, fc, fd, fe, ff},
					{v, 0, vfinal}]

(* Extract the solutions for plotting and further analysis *)
far = fa /. sol[[1]]; fbr = fb /. sol[[1]];  fcr = fc /. sol[[1]];
fdr = fd /. sol[[1]]; fer = fe /. sol[[1]];  ffr = ff /. sol[[1]];

(* Determine the mole fractions              *)
yc[v_] := fcr[v]/(far[v] + fbr[v] + fcr[v] + fdr[v] + fer[v] 
		+ ffr[v] + finert);

g1 = Plot[{far[v], fbr[v], fcr[v], ffr[v]}, {v, 0,vfinal},
		GridLines -> Automatic, Frame -> True, AspectRatio -> 1,
		FrameLabel -> {"Volume (cm^3)", "Molar Flow (mol/s)"}]

Plot[{fdr[v], fer[v], ffr[v]}, {v,0,vfinal},
		GridLines -> Automatic, Frame -> True, AspectRatio -> 1,
		FrameLabel -> {"Volume (cm^3)", "Molar Flow (mol/s)"}]

(* Find the vapor pressure to determine when methanol begins *)
(* to condense                                          	 *)

harlach = Log[pvp] - (hara + harb/t + harc Log[t] + hard pvp/t^2);
hara = 72.268; harb = -7064.2; harc = -7.68; hard = -1.86;

pmeohsat =  pvp /.  FindRoot[harlach == 0, {pvp, 100, 10000}]

ycsat = (pmeohsat/760)/p

Plot[{yc[v], ycsat},{v,0, 200},
		GridLines -> Automatic, Frame -> True, AspectRatio -> 1,
		FrameLabel -> {"Volume (cm^3)", "MeOH Mole Fr."}]

vcrit = v /. FindRoot[yc[v] - ycsat == 0, {v, {1, 50}}]

(*  Now re-integrate the differential equations from vcrit		*)
(*  to vfinal													*)
(*	We must, now, correct sumfi to account for the gas phase	*)
(*	components only.											*)

sumfi = (faa[v] + fba[v]  + fda[v] + fea[v] + ffa[v] + finert)/
				(1 - ycsat);
ca = faa[v]/vdot; cb = fba[v]/vdot; 
cc = ycsat p/(82.1 t);
cd = fda[v]/vdot; ce = fea[v]/vdot; cf = ffa[v]/vdot;
vdot = vdoto (sumfi/ftoto);

favcrit = far[vcrit]; fbvcrit = fbr[vcrit]; fcvcrit = fcr[vcrit];
fdvcrit = fdr[vcrit]; fevcrit = fer[vcrit]; ffvcrit = ffr[vcrit];

sol2 =	NDSolve[{	faa'[v] == -r1	-r2,
					fba'[v] == -2r1	+r2		+r3,
					fca'[v] == +r1			-r3,
					fda'[v] == 		-r2,
					fea'[v] ==		+r2,
					ffa'[v] ==				+r3,
					faa[vcrit] == favcrit, 	fba[vcrit] == fbvcrit,
					fca[vcrit] == fcvcrit, 	fda[vcrit] == fdvcrit,
					fea[vcrit] == fevcrit, 	ffa[vcrit] == ffvcrit},
					{faa, fba, fca, fda, fea, ffa},
					{v, vcrit, vfinal}]

faar = faa /. sol2[[1]]; fbar = fba /. sol2[[1]]; 
fcar = fca /. sol2[[1]]; fdar = fda /. sol2[[1]]; 
fear = fea /. sol2[[1]]; ffar = ffa /. sol2[[1]];
ScientificForm[fcr[vfinal]]
ScientificForm[fcar[vfinal]]

g2 = Plot[{faar[v], fbar[v], fcar[v], ffar[v]}, {v,vcrit, vfinal},
		GridLines -> Automatic, Frame -> True, AspectRatio -> 1,
		FrameLabel -> {"Volume (cm^3)", "Molar Flow (mol/s)"}]

(*  Plot the results, with and without condensation, together *)
Show[g1, g2]