Difference between revisions of "Isothermal CSTR"

From Sutherland_wiki
Jump to: navigation, search
(beta1)
(Replaced content with "<h3>CSTR Reactor</h3> <h5>By Taylor Geisler</h5>")
Line 1: Line 1:
 
+
<h3>CSTR Reactor</h3>
        <script language="JavaScript" src="numeric-1.2.6.js"></script>
+
<h5>By Taylor Geisler</h5>
        <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 + " &ndash; ");
 
              });
 
        }
 
 
 
       
 
       
 
////////////////////////////////////////////////////////////////////////
 
 
 
            //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>
 
        <h3>CSTR Reactor</h3>
 
        <h5>By Taylor Geisler</h5>
 
        <img src="batch.png" width="100" height="166">
 
        <div id="placeholder" style="width:900px;height:300px"></div>
 
        <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<br>
 
            <INPUT type="text" name= "Kcform" value = "1.2">Kc<br>
 
            <INPUT type="text" name= "Tform" value="289">Temperature (K)<br>
 
            <INPUT type="text" name= "tstepform" value = "0.1">Time Step (min)<br>
 
            <INPUT type="text" name= "tfinalform" value = "1000">Final Time (min)<br>
 
            <INPUT type="text" name= "V0form" value="100">Volume of Reactor (m^3)<br>
 
               
 
            <INPUT type="text" name= "Ca0form" value="4">Inital Concentration of A in Reactor<br>
 
            <INPUT type="text" name= "Cb0form" value = "3">Inital Concentration of B in Reactor<br>
 
            <INPUT type="text" name= "Cc0form" value="2">Inital Concentration of C in Reactor<br>
 
            <INPUT type="text" name= "Cd0form" value = "1">Inital Concentration of D in Reactor<br>
 
               
 
            <INPUT type="text" name= "Cainform" value = "1">Inlet Concentration of A<br>
 
            <INPUT type="text" name= "Cbinform" value="1">Inlet Concentration of B<br>
 
            <INPUT type="text" name= "Ccinform" value="1">Inlet Concentration of C<br>
 
            <INPUT type="text" name= "Cdinform" value = "1">Inlet Concentration of D<br>
 
               
 
            <INPUT type="text" name= "finform" value="1">Flow In<br>
 
            <INPUT type="text" name= "foutform" value = "1">Flow Out<br>
 
        </FORM>
 
       
 
    </body>
 

Revision as of 14:55, 16 November 2013

CSTR Reactor

By Taylor Geisler