Add support for multiple reactor types (CSTR and PFR) with corresponding properties (Dichelet BC for now)

advanced-reactor.html

@@ -4,7 +4,10 @@
         color: "#c4cce0",
         defaults: {
             name: { value: "" },
-            volume: { value: 0., required: true},
+            reactor_type: { value: "CSTR", required: true },
+            volume: { value: 0., required: true },
+            length: { value: 0.},
+            resolution_L: { value: 0.},
             n_inlets: { value: 1, required: true},
             kla: { value: null },
             S_O_init: { value: 0., required: true },
@@ -45,6 +48,32 @@
                 type:"num",
                 types:["num"]
             });
+            $("#node-input-reactor_type").typedInput({
+                types: [
+                    {
+                        value: "CSTR",
+                        options: [
+                            { value: "CSTR", label: "CSTR"},
+                            { value: "PFR", label: "PFR"}
+                        ]
+                    }
+                ]
+            })
+            $("#node-input-reactor_type").on("change", function() {
+                const type = $("#node-input-reactor_type").typedInput("value");
+                if (type === "CSTR") {
+                    $(".PFR").hide();
+                } else {
+                    $(".PFR").show();
+                }
+            });
+            // Set initial visibility on dialog open
+            const initialType = $("#node-input-reactor_type").typedInput("value");
+            if (initialType === "CSTR") {
+                $(".PFR").hide();
+            } else {
+                $(".PFR").show();
+            }
         },
         oneditsave: function() {
             let volume = parseFloat($("#node-input-volume").typedInput("value"));
@@ -65,10 +94,22 @@
         <input type="text" id="node-input-name" placeholder="Name">
     </div>
     <h2> Reactor properties </h2>
+    <div class="form-row">
+        <label for="node-input-reactor_type"><i class="fa fa-tag"></i> Reactor type</label>
+        <input type="text" id="node-input-reactor_type">
+    </div>
     <div class="form-row">
         <label for="node-input-volume"><i class="fa fa-tag"></i> Fluid volume [m3]</label>
         <input type="text" id="node-input-volume" placeholder="m3">
     </div>
+    <div class="form-row PFR">
+        <label for="node-input-length"><i class="fa fa-tag"></i> Reactor length [m]</label>
+        <input type="text" id="node-input-length" placeholder="m">
+    </div>
+    <div class="form-row PFR">
+        <label for="node-input-resolution_L"><i class="fa fa-tag"></i> Resolution</label>
+        <input type="text" id="node-input-resolution_L" placeholder="#">
+    </div>
     <div class="form-row">
         <label for="node-input-n_inlets"><i class="fa fa-tag"></i> Number of inlets</label>
         <input type="text" id="node-input-n_inlets" placeholder="#">

advanced-reactor.js

@@ -7,26 +7,60 @@ module.exports = function(RED) {
 
         const Reactor = require('./dependencies/reactor_class');
 
-        const reactor = new Reactor(
+        let new_reactor;
-            parseFloat(config.volume),
+
-            parseInt(config.n_inlets),
+        switch (config.reactor_type) {
-            parseFloat(config.kla),
+            case "CSTR":
-            [
+                new_reactor = new Reactor(
-                parseFloat(config.S_O_init),
+                    parseFloat(config.volume),
-                parseFloat(config.S_I_init),
+                    parseInt(config.n_inlets),
-                parseFloat(config.S_S_init),
+                    parseFloat(config.kla),
-                parseFloat(config.S_NH_init),
+                    [
-                parseFloat(config.S_N2_init),
+                        parseFloat(config.S_O_init),
-                parseFloat(config.S_NO_init),
+                        parseFloat(config.S_I_init),
-                parseFloat(config.S_HCO_init),
+                        parseFloat(config.S_S_init),
-                parseFloat(config.X_I_init),
+                        parseFloat(config.S_NH_init),
-                parseFloat(config.X_S_init),
+                        parseFloat(config.S_N2_init),
-                parseFloat(config.X_H_init),
+                        parseFloat(config.S_NO_init),
-                parseFloat(config.X_STO_init),
+                        parseFloat(config.S_HCO_init),
-                parseFloat(config.X_A_init),
+                        parseFloat(config.X_I_init),
-                parseFloat(config.X_TS_init)
+                        parseFloat(config.X_S_init),
-            ]
+                        parseFloat(config.X_H_init),
-        );
+                        parseFloat(config.X_STO_init),
+                        parseFloat(config.X_A_init),
+                        parseFloat(config.X_TS_init)
+                    ]
+                );
+                break;
+            case "PFR":
+                new_reactor = new Reactor(
+                    parseFloat(config.volume),
+                    parseFloat(config.L),
+                    parseInt(config.resolution_L),
+                    parseInt(config.n_inlets),
+                    parseFloat(config.kla),
+                    [
+                        parseFloat(config.S_O_init),
+                        parseFloat(config.S_I_init),
+                        parseFloat(config.S_S_init),
+                        parseFloat(config.S_NH_init),
+                        parseFloat(config.S_N2_init),
+                        parseFloat(config.S_NO_init),
+                        parseFloat(config.S_HCO_init),
+                        parseFloat(config.X_I_init),
+                        parseFloat(config.X_S_init),
+                        parseFloat(config.X_H_init),
+                        parseFloat(config.X_STO_init),
+                        parseFloat(config.X_A_init),
+                        parseFloat(config.X_TS_init)
+                    ]
+                );
+                break;
+            default:
+                console.warn("Unknown reactor type: " + config.reactor_type);
+        }
+
+        const reactor = new_reactor; // protect from reassignment
 
         node.on('input', function(msg, send, done) {
             let toggleUpdate = false;

dependencies/reactor_class.js

@@ -5,7 +5,6 @@ class Reactor_CSTR {
     
     constructor(volume, n_inlets, kla, initial_state) {
         this.state = initial_state;
-        console.log(this.state);
         this.asm = new ASM3();
 
         this.Vl = volume; // fluid volume reactor [m3]
@@ -17,7 +16,7 @@ class Reactor_CSTR {
 
         this.currentTime = Date.now(); // milliseconds since epoch [ms]
         this.timeStep = 1/(24*60*15); // time step [d]
-        this.speedUpFactor = 30;
+        this.speedUpFactor = 1;
     }
 
     set setInfluent(input) { // setter for C_in (WIP)
@@ -69,6 +68,113 @@ class Reactor_CSTR {
     }
 }
 
+class Reactor_PFR {
+    
+    constructor(volume, length, resolution_L, n_inlets, kla, initial_state) {
+        this.asm = new ASM3();
+
+        this.Vl = volume; // fluid volume reactor [m3]
+        this.length = length; // reactor length [m]
+        this.n_x = resolution_L; // number of slices
+        this.d_x = length / resolution_L;
+
+        this.A = volume / length; // crosssectional area [m2]
+
+        this.state = Array.from(Array(this.n_x), () => initial_state.slice())
+
+        this.Fs = Array(n_inlets).fill(0.0); // fluid debits per inlet [m3 d-1]
+        this.Cs_in = Array.from(Array(n_inlets), () => new Array(13).fill(0.0)); // composition influents
+        this.OTR = 0.0; // oxygen transfer rate [g O2 d-1]
+        this.D = 0.0; // axial dispersion [m2 d-1]
+        
+        this.kla = kla; // if NaN, use external OTR [d-1]
+
+        this.currentTime = Date.now(); // milliseconds since epoch [ms]
+        this.timeStep = 1/(24*60*15); // time step [d]
+        this.speedUpFactor = 1;
+
+        this.D_op = makeDoperator();
+        this.D2_op = makeD2operator();
+    }
+
+    set setInfluent(input) { // setter for C_in (WIP)
+        let index_in = input.payload.inlet;
+        this.Fs[index_in] = input.payload.F;
+        this.Cs_in[index_in] = input.payload.C;
+    }
+
+    set setOTR(input) { // setter for OTR (WIP) [g O2 d-1]
+        this.OTR = input.payload;
+    }
+
+    set setDispersion(input) { // setter for Axial dispersion [m2 d-1]
+        this.D = input.payload;
+    }
+
+    get getEffluent() { // getter for Effluent, defaults to inlet 0
+        return {topic: "Fluent", payload: {inlet: 0, F: math.sum(this.Fs), C:this.state}, timestamp: this.currentTime};
+    }
+
+    calcOTR(S_O, T=20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
+        let S_O_sat = 14.652 - 4.1022e-1*T + 7.9910e-3*T*T + 7.7774e-5*T*T*T;
+        return this.kla * (S_O_sat - S_O);
+    }
+
+    // expect update with timestamp
+    updateState(newTime) {
+
+        const day2ms = 1000 * 60 * 60 * 24; 
+
+        let n_iter = Math.floor(this.speedUpFactor*(newTime - this.currentTime) / (this.timeStep * day2ms));
+        if (n_iter) {
+            let n = 0;
+            while (n < n_iter) {
+                this.tick_fe(this.timeStep);
+                n += 1;
+            }
+            this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
+        }
+    }
+
+    tick_fe(time_step) { // tick reactor state using forward Euler method
+        if (math.sum(this.Fs) > 0) {
+            this.state[0] = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0] // Dichelet boundary condition
+        }
+
+        const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
+        const advection = math.multiply(math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
+        const reaction = this.state.map(this.asm.compute_dC);
+        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
+            transfer.forEach((x) => { x[0] = this.OTR; });
+        } else {
+            transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); });
+        }
+
+        const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
+
+        this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
+        return this.state;
+    }
+
+    makeDoperator() { // create the upwind scheme gradient operator
+        const I = math.identity(this.n_x);
+        const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
+        I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L
+        return math.add(I, A);
+    }
+
+    makeD2operator() { // create the upwind scheme second derivative operator
+        const I = math.diag(Array(this.n_x).fill(2), 0);
+        const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
+        const B = math.diag(Array(this.n_x).fill(-1), -1).resize([this.n_x, this.n_x]);
+        I[0, 0] = 1;
+        return math.add(I, A, B);
+    }
+}
+
+
 // testing stuff
 // state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
 // let initial_state = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1];