Isothermal CSTR
From Sutherland_wiki
<!DOCTYPE html> <html>
<head>
<script language="JavaScript" src="numeric-1.2.6.js"></script>
<script language="javascript" type="text/javascript" src="jquery.js"></script>
<script language="javascript" type="text/javascript" src="jquery.flot.js"></script>
<script language="JavaScript">
//Rate Equations
var ratelaw = function(CaR, CbR, CcR, CdR, kR) {
var rA = kR * CaR * CbR;
return -rA;
}
var f1 = function(Ca1F, Cb1F, Cc1F, Cd1F, kF, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin) {
var f1x = -Ca1F + Ca[i] + tstep*(ratelaw(Ca1F, Cb1F, Cc1F, Cd1F, kF)-Ca1F*fout[i]/V[i] + Cain*fin[i]/V[i]);
return f1x;
}
var f2 = function(Ca1F, Cb1F, Cc1F, Cd1F, kF, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin) {
var f2x = -Cb1F + Cb[i] + tstep*(ratelaw(Ca1F, Cb1F, Cc1F, Cd1F, kF)-Cb1F*fout[i]/V[i] + Cbin*fin[i]/V[i]);
return f2x;
}
var f3 = function(Ca1F, Cb1F, Cc1F, Cd1F, kF, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin) {
var f3x = -Cc1F + Cc[i] + tstep*(-ratelaw(Ca1F, Cb1F, Cc1F, Cd1F, kF)-Cc1F*fout[i]/V[i] + Ccin*fin[i]/V[i]);
return f3x;
}
var f4 = function(Ca1F, Cb1F, Cc1F, Cd1F, kF, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin) {
var f4x = -Cd1F + Cd[i] + tstep*(-ratelaw(Ca1F, Cb1F, Cc1F, Cd1F, kF)-Cd1F*fout[i]/V[i] + Cdin*fin[i]/V[i]);
return f4x;
}
var derriv = function(fname,Ca1, Cb1, Cc1, Cd1, deltaCa, deltaCb, deltaCc, deltaCd, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin) {
var A = fname(Ca1,Cb1,Cc1,Cd1, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
var B = fname(Ca1+deltaCa, Cb1+deltaCb, Cc1+deltaCc, Cd1 + deltaCd, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
var der = (B - A) / stepsize;
return der;
}
var jacob = function(CaguessJ, CbguessJ, CcguessJ, CdguessJ, stepsizeJ, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin) {
//Jacobian
JJ = numeric.rep([4,4],0);
JJ[0][0] = derriv(f1, CaguessJ, CbguessJ, CcguessJ, CdguessJ, stepsizeJ, 0, 0, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[0][1] = derriv(f1, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, stepsizeJ, 0, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[0][2] = derriv(f1, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, 0, stepsizeJ, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[0][3] = derriv(f1, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, 0, 0, stepsizeJ, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[1][0] = derriv(f2, CaguessJ, CbguessJ, CcguessJ, CdguessJ, stepsizeJ, 0, 0, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[1][1] = derriv(f2, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, stepsizeJ, 0, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[1][2] = derriv(f2, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, 0, stepsizeJ, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[1][3] = derriv(f2, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, 0, 0, stepsizeJ, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[2][0] = derriv(f3, CaguessJ, CbguessJ, CcguessJ, CdguessJ, stepsizeJ, 0, 0, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[2][1] = derriv(f3, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, stepsizeJ, 0, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[2][2] = derriv(f3, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, 0, stepsizeJ, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[2][3] = derriv(f3, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, 0, 0, stepsizeJ, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[3][0] = derriv(f4, CaguessJ, CbguessJ, CcguessJ, CdguessJ, stepsizeJ, 0, 0, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[3][1] = derriv(f4, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, stepsizeJ, 0, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[3][2] = derriv(f4, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, 0, stepsizeJ, 0, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
JJ[3][3] = derriv(f4, CaguessJ, CbguessJ, CcguessJ, CdguessJ, 0, 0, 0, stepsizeJ, kD, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
return JJ;
}
var RHSfun = function(CaguessR, CbguessR, CcguessR, CdguessR, kR, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin) {
RHSR = [f1(CaguessR, CbguessR, CcguessR, CdguessR, kR, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin), f2(CaguessR, CbguessR, CcguessR, CdguessR, kR, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin), f3(CaguessR, CbguessR, CcguessR, CdguessR, kR, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin), f4(CaguessR, CbguessR, CcguessR, CdguessR, kR, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin)];
return RHSR;
}
var NRHSfun = function(RHSN) {
nRHSR = [-1 * RHSN[0], -1 * RHSN[1], -1* RHSN[2], -1 * RHSN[3]];
return nRHSR;
}
var mainFunction = function(tstep, tfinal, V0, Ca0, Cb0, Cc0, Cd0, Cain, Cbin, Ccin, Cdin, fin0, fout0, T, k, Kc) {
//Inside Function
var numsteps = (tfinal / tstep) + 1;
var t = new Array(numsteps);
t[0] = 0;
var V = new Array(numsteps);
V[0] = V0;
var Ca = new Array(numsteps);
Ca[0] = Ca0;
var Cb = new Array(numsteps);
Cb[0] = Cb0;
var Cc = new Array(numsteps);
Cc[0] = Cc0;
var Cd = new Array(numsteps);
Cd[0] = Cd0;
var fin = new Array(numsteps);
fin[0] = fin0;
var fout = new Array(numsteps);
fout[0] = fout0;
var r;
var compare = new Array(3);
var norm;
var Caguess = 1;
var Cbguess = 1;
var Ccguess = 1;
var Cdguess = 1;
//////////////////////////////////////////////
for (var i=0; i<numsteps; i++) {
V[i+1] = V[i] + tstep*(fin0-fout0);
fin[i+1] = fin[i];
fout[i+1] = fout[i];
t[i+1] = t[i]+tstep;
//Make initial guesses with forward euler
r = ratelaw(Ca[i], Cb[i], Cc[i], Cd[i], k);
var Caguess = 1 + Ca[i] + tstep*(r - Ca[i]*fout[i]/V[i] + Cain*fin[i]/V[i]);
var Cbguess = 3 + Cb[i] + tstep*(r - Cb[i]*fout[i]/V[i] + Cbin*fin[i]/V[i]);
var Ccguess = 9 + Cc[i] + tstep*(-r - Cc[i]*fout[i]/V[i] + Ccin*fin[i]/V[i]);
var Cdguess = -0.1 + Cd[i] + tstep*(-r - Cd[i]*fout[i]/V[i] + Cdin*fin[i]/V[i]);
RHS = RHSfun(Caguess, Cbguess, Ccguess, Cdguess, k, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
nRHS = NRHSfun(RHS);
//Jacobian
J = jacob(Caguess, Cbguess, Ccguess, Cdguess, stepsize, k, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
// console.log(numeric.prettyPrint(RHS)); // console.log(numeric.prettyPrint(J));
//J*deltax = RHS
dlta = numeric.solve(J,nRHS);
// console.log(numeric.prettyPrint(dlta));
//x = x + deltax
Caguess = Caguess + dlta[0];
Cbguess = Cbguess + dlta[1];
Ccguess = Ccguess + dlta[2];
Cdguess = Cdguess + dlta[3];
RHS = RHSfun(Caguess, Cbguess, Ccguess, Cdguess, k, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
nRHS = NRHSfun(RHS);
norm = numeric.norm2(numeric.sub(RHS));
// console.log(norm);
for (var j = 0; j<5; j++) {
//Jacobian
J = jacob(Caguess, Cbguess, Ccguess, Cdguess, stepsize, k, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
//J*deltax = RHS
dlta = numeric.solve(J,nRHS);
//x = x + deltax
Caguess = Caguess + dlta[0];
Cbguess = Cbguess + dlta[1];
Ccguess = Ccguess + dlta[2];
Cdguess = Cdguess + dlta[3];
RHS = RHSfun(Caguess, Cbguess, Ccguess, Cdguess, k, Ca, Cb, Cc, Cd, i, fout, fin, V, tstep, Cain, Cbin, Ccin, Cdin);
nRHS = NRHSfun(RHS);
norm = numeric.norm2(numeric.sub(RHS));
// if (j === 0) {
console.log(norm);
// }
}
Ca[i+1] = Caguess;
Cb[i+1] = Cbguess;
Cc[i+1] = Ccguess;
Cd[i+1] = Cdguess;
}
console.log('Ca final = ' + Ca[numsteps]);
console.log('Cb final = ' + Cb[numsteps]);
console.log('Cc final = ' + Cc[numsteps]);
console.log('Cd final = ' + Cd[numsteps]);
console.log('The total concentration is ' + (Ca[numsteps]+Cb[numsteps]+Cc[numsteps]+Cd[numsteps]));
//Plot The thing
$(function() {
var d1 = [];
for (var i = 0; i < t.length; i ++) {
d1.push([t[i], Ca[i]]);
}
var d2 = [];
for (var i = 0; i < t.length; i ++) {
d2.push([t[i], Cb[i]]);
}
var d3 = [];
for (var i = 0; i < t.length; i ++) {
d3.push([t[i], Cc[i]]);
}
var d4 = [];
for (var i = 0; i < t.length; i ++) {
d4.push([t[i], Cd[i]]);
}
// $.plot("#placeholder", [d1, d2, d3]);
$.plot("#placeholder", [
{ label: "Ca", data: d1 },
{ label: "Cb", data: d2 },
{ label: "Cc", data: d3 },
{ label: "Cd", data: d4 }
]);
$("#footer").prepend("Flot " + $.plot.version + " – ");
});
}
////////////////////////////////////////////////////////////////////////
//Initialization
//Outside Function
// (tstep, tfinal, V0, Ca0, Cb0, Cc0, Cd0, Cain, Cbin, Ccin, Cdin, fin0, fout0, T, k, Kc)
// var tstep = 0.1; // var tfinal = 1000; // // var V0 = 100; // var Ca0 = 4; // var Cb0 = 3; // var Cc0 = 2; // var Cd0 = 1;
// var Cain = 1; // var Cbin = 1; // var Ccin = 1; // var Cdin = 1; // // var fin0 = 1; // var fout0 = 1;
// var T = 298;
// var k = 0.01; // var Kc = 1.2; //
var stepsize = 0.0001;
var tol = 0.0000000001;
</script> </head> <body>
CSTR Reactor
By Taylor Geisler
<img src="batch.png" width="100" height="166">
<button onclick="mainFunction(parseFloat(document.param.tstepform.value), parseFloat(document.param.tfinalform.value), parseFloat(document.param.V0form.value), parseFloat(document.param.Ca0form.value), parseFloat(document.param.Cb0form.value), parseFloat(document.param.Cc0form.value), parseFloat(document.param.Cd0form.value), parseFloat(document.param.Cainform.value), parseFloat(document.param.Cbinform.value), parseFloat(document.param.Ccinform.value), parseFloat(document.param.Cdinform.value), parseFloat(document.param.finform.value), parseFloat(document.param.foutform.value), parseFloat(document.param.Tform.value), parseFloat(document.param.kform.value), parseFloat(document.param.Kcform.value))">Update Plot</button>
<FORM name = "param">
<INPUT type="text" name= "kform" value="0.01">k
<INPUT type="text" name= "Kcform" value = "1.2">Kc
<INPUT type="text" name= "Tform" value="289">Temperature (K)
<INPUT type="text" name= "tstepform" value = "0.1">Time Step (min)
<INPUT type="text" name= "tfinalform" value = "1000">Final Time (min)
<INPUT type="text" name= "V0form" value="100">Volume of Reactor (m^3)
<INPUT type="text" name= "Ca0form" value="4">Inital Concentration of A in Reactor
<INPUT type="text" name= "Cb0form" value = "3">Inital Concentration of B in Reactor
<INPUT type="text" name= "Cc0form" value="2">Inital Concentration of C in Reactor
<INPUT type="text" name= "Cd0form" value = "1">Inital Concentration of D in Reactor
<INPUT type="text" name= "Cainform" value = "1">Inlet Concentration of A
<INPUT type="text" name= "Cbinform" value="1">Inlet Concentration of B
<INPUT type="text" name= "Ccinform" value="1">Inlet Concentration of C
<INPUT type="text" name= "Cdinform" value = "1">Inlet Concentration of D
<INPUT type="text" name= "finform" value="1">Flow In
<INPUT type="text" name= "foutform" value = "1">Flow Out
</FORM>
</body>
</html>
