Difference between revisions of "Isothermal CSTR"

Revision as of 14:45, 16 November 2013

       <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>

   
   
   <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>