Refactor dispersion and boundary condition handling in Reactor_PFR

dependencies/reactor_class.js

@@ -103,7 +103,7 @@ class Reactor_PFR {
         this.timeStep = 1/(24*60*15); // time step [d]
         this.speedUpFactor = 60;
 
-        this.D_op = this.makeDoperator();
+        this.D_op = this.makeDoperator(false);
         this.D2_op = this.makeD2operator();
     }
 
@@ -111,6 +111,8 @@ class Reactor_PFR {
         let index_in = input.payload.inlet;
         this.Fs[index_in] = input.payload.F;
         this.Cs_in[index_in] = input.payload.C;
+        // console.log("Pe " + this.d_x*math.sum(this.Fs)/(this.D*this.A));
+        // console.log("Co " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
     }
 
     set setOTR(input) { // setter for OTR (WIP) [g O2 d-1]
@@ -150,7 +152,6 @@ class Reactor_PFR {
         const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
         const advection = math.multiply(-1*math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
         const reaction = this.state.map((state_slice) => this.asm.compute_dC(state_slice));
-        reaction[0] = Array(13).fill(0.0);
         const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
         
         if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
@@ -159,41 +160,45 @@ class Reactor_PFR {
             transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); });
         }
 
-        transfer[0][0] = 0;
+        const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
+        const new_state = math.add(this.state, dC_total);
 
+        // apply boundary conditions
         if (math.sum(this.Fs) > 0) { // Danckwerts BC
-            const BC_C_in = math.multiply(1/math.sum(this.Fs),[this.Fs], this.Cs_in)[0];
+            const BC_C_in = math.multiply(1/math.sum(this.Fs), [this.Fs], this.Cs_in)[0];
             const BC_gradient = Array(this.n_x).fill(0.0);
             BC_gradient[0] = -1;
             BC_gradient[1] = 1;
-            const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], this.state)[0];
+            const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], new_state)[0];
-            this.state[0] = math.add(BC_C_in, BC_dispersion);
+            new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
-            console.log(BC_dispersion);
+            console.log(new_state[0])
+        } else { // Neumann BC (no flux)
+            new_state[0] = new_state[1];
         }
+        // Neumann BC (no flux)
+        new_state[this.n_x-1] = new_state[this.n_x-2]
 
-        const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
+        this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state
-
+        return new_state;
-        // clip value element-wise to each subarray to avoid negative concentrations
-        this.state = math.add(this.state, dC_total).map(row => row.map(val => val < 0 ? 0 : val));
-        return this.state;
     }
 
-    makeDoperator() { // create the upwind scheme gradient operator
+    makeDoperator(central=false) { // create the upwind scheme gradient operator
-        const I = math.diag(Array(this.n_x).fill(1), 0);
+        const I = math.resize(math.diag(Array(this.n_x).fill(1), central), [this.n_x, this.n_x]);
         const A = math.resize(math.diag(Array(this.n_x).fill(-1), -1), [this.n_x, this.n_x]);
-        I[0][0] = 0;
+        const D = math.add(I, A);
-        I[0][1] = -1;
+        D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere
-        return math.add(I, A);
+        D[this.n_x-1] = Array(this.n_x).fill(0);
+        return D;
     }
 
     makeD2operator() { // create the central second derivative operator
         const I = math.diag(Array(this.n_x).fill(-2), 0);
         const A = math.resize(math.diag(Array(this.n_x).fill(1), 1), [this.n_x, this.n_x]);
         const B = math.resize(math.diag(Array(this.n_x).fill(1), -1), [this.n_x, this.n_x]);
-        I[0][0] = 0;
+        const D2 = math.add(I, A, B);
-        I[0][1] = -1;
+        D2[0] = Array(this.n_x).fill(0); // set by BCs elsewhere
-        I[0][0] = -1; // Dichelet boundary condition at outlet
+        D2[this.n_x-1] = Array(this.n_x).fill(0);
-        return math.add(I, A, B);
+        return D2;
     }
 }
 
@@ -211,4 +216,4 @@ class Reactor_PFR {
 //     N += 1;
 // }
 
-module.exports = {Reactor_CSTR, Reactor_PFR};
+module.exports = { Reactor_CSTR, Reactor_PFR };