diff --git a/advanced-reactor.html b/advanced-reactor.html
index 73d66c2..512f3dd 100644
--- a/advanced-reactor.html
+++ b/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 },
@@ -37,6 +40,14 @@
type:"num",
types:["num"]
});
+ $("#node-input-length").typedInput({
+ type:"num",
+ types:["num"]
+ });
+ $("#node-input-resolution_L").typedInput({
+ type:"num",
+ types:["num"]
+ });
$("#node-input-kla").typedInput({
type:"num",
types:["num"]
@@ -45,6 +56,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 +102,22 @@
Reactor properties
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/advanced-reactor.js b/advanced-reactor.js
index edfad14..8348137 100644
--- a/advanced-reactor.js
+++ b/advanced-reactor.js
@@ -5,28 +5,62 @@ module.exports = function(RED) {
let name = config.name;
- const Reactor = require('./dependencies/reactor_class');
+ const { Reactor_CSTR, Reactor_PFR } = require('./dependencies/reactor_class');
- const reactor = new Reactor(
- parseFloat(config.volume),
- 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)
- ]
- );
+ let new_reactor;
+
+ switch (config.reactor_type) {
+ case "CSTR":
+ new_reactor = new Reactor_CSTR(
+ parseFloat(config.volume),
+ 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;
+ case "PFR":
+ new_reactor = new Reactor_PFR(
+ parseFloat(config.volume),
+ parseFloat(config.length),
+ 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;
diff --git a/dependencies/asm3_class.js b/dependencies/asm3_class.js
index 9ad03d9..9b7762a 100644
--- a/dependencies/asm3_class.js
+++ b/dependencies/asm3_class.js
@@ -2,67 +2,65 @@ const math = require('mathjs')
class ASM3 {
- kin_params = {
- // Kinetic parameters (20 C for now)
-
- // Hydrolysis
- k_H: 3., // hydrolysis rate constant [g X_S g-1 X_H d-1]
- K_X: 1., // hydrolysis saturation constant [g X_S g-1 X_H]
- // Heterotrophs
- k_STO: 5., // storage rate constant [g S_S g-1 X_H d-1]
- nu_NO: 0.6, // anoxic reduction factor [-]
- K_O: 0.2, // saturation constant S_0 [g O2 m-3]
- K_NO: 0.5, // saturation constant S_NO [g NO3-N m-3]
- K_S: 2., // saturation constant S_s [g COD m-3]
- K_STO: 1., // saturation constant X_STO [g X_STO g-1 X_H]
- mu_H_max: 2., // maximum specific growth rate [d-1]
- K_NH: 0.01, // saturation constant S_NH3 [g NH3-N m-3]
- K_HCO: 0.1, // saturation constant S_HCO [mole HCO3 m-3]
- b_H_O: 0.2, // aerobic respiration rate [d-1]
- b_H_NO: 0.1, // anoxic respiration rate [d-1]
- b_STO_O: 0.2, // aerobic respitation rate X_STO [d-1]
- b_STO_NO: 0.1, // anoxic respitation rate X_STO [d-1]
- // Autotrophs
- mu_A_max: 1.0, // maximum specific growth rate [d-1]
- K_A_NH: 1., // saturation constant S_NH3 [g NH3-N m-3]
- K_A_O: 0.5, // saturation constant S_0 [g O2 m-3]
- K_A_HCO: 0.5, // saturation constant S_HCO [mole HCO3 m-3]
- b_A_O: 0.15, // aerobic respiration rate [d-1]
- b_A_NO: 0.05 // anoxic respiration rate [d-1]
- }
-
- stoi_params = {
- // Stoichiometric and composition parameters
-
- f_SI: 0., // fraction S_I from hydrolysis [g S_I g-1 X_S]
- f_XI: 0.2, // fraction X_I from decomp X_H [g X_I g-1 X_H]
- // Yields
- Y_STO_O: 0.85, // aerobic yield X_STO per S_S [g X_STO g-1 S_S]
- Y_STO_NO: 0.80, // anoxic yield X_STO per S_S [g X_STO g-1 S_S]
- Y_H_O: 0.63, // aerobic yield X_H per X_STO [g X_H g-1 X_STO]
- Y_H_NO: 0.54, // anoxic yield X_H per X_STO [g X_H g-1 X_STO]
- Y_A: 0.24, // anoxic yield X_A per S_NO [g X_A g-1 NO3-N]
- // Composition (COD via DoR)
- i_CODN: -1.71, // COD content (DoR) [g COD g-1 N2-N]
- i_CODNO: -4.57, // COD content (DoR) [g COD g-1 NO3-N]
- // Composition (nitrogen)
- i_NSI: 0.01, // nitrogen content S_I [g N g-1 S_I]
- i_NSS: 0.03, // nitrogen content S_S [g N g-1 S_S]
- i_NXI: 0.02, // nitrogen content X_I [g N g-1 X_I]
- i_NXS: 0.04, // nitrogen content X_S [g N g-1 X_S]
- i_NBM: 0.07, // nitrogen content X_H / X_A [g N g-1 X_H / X_A]
- // Composition (TSS)
- i_TSXI: 0.75, // TSS content X_I [g TS g-1 X_I]
- i_TSXS: 0.75, // TSS content X_S [g TS g-1 X_S]
- i_TSBM: 0.90, // TSS content X_H / X_A [g TS g-1 X_H / X_A]
- i_TSSTO: 0.60, // TSS content X_STO (PHB based) [g TS g-1 X_STO]
- // Composition (charge)
- i_cNH: 1/14, // charge per S_NH [mole H+ g-1 NH3-N]
- i_cNO: -1/14 // charge per S_NO [mole H+ g-1 NO3-N]
- }
-
constructor() {
- this.stoi_matrix = this._initialise_stoi_matrix()
+ this.kin_params = {
+ // Kinetic parameters (20 C for now)
+
+ // Hydrolysis
+ k_H: 3., // hydrolysis rate constant [g X_S g-1 X_H d-1]
+ K_X: 1., // hydrolysis saturation constant [g X_S g-1 X_H]
+ // Heterotrophs
+ k_STO: 5., // storage rate constant [g S_S g-1 X_H d-1]
+ nu_NO: 0.6, // anoxic reduction factor [-]
+ K_O: 0.2, // saturation constant S_0 [g O2 m-3]
+ K_NO: 0.5, // saturation constant S_NO [g NO3-N m-3]
+ K_S: 2., // saturation constant S_s [g COD m-3]
+ K_STO: 1., // saturation constant X_STO [g X_STO g-1 X_H]
+ mu_H_max: 2., // maximum specific growth rate [d-1]
+ K_NH: 0.01, // saturation constant S_NH3 [g NH3-N m-3]
+ K_HCO: 0.1, // saturation constant S_HCO [mole HCO3 m-3]
+ b_H_O: 0.2, // aerobic respiration rate [d-1]
+ b_H_NO: 0.1, // anoxic respiration rate [d-1]
+ b_STO_O: 0.2, // aerobic respitation rate X_STO [d-1]
+ b_STO_NO: 0.1, // anoxic respitation rate X_STO [d-1]
+ // Autotrophs
+ mu_A_max: 1.0, // maximum specific growth rate [d-1]
+ K_A_NH: 1., // saturation constant S_NH3 [g NH3-N m-3]
+ K_A_O: 0.5, // saturation constant S_0 [g O2 m-3]
+ K_A_HCO: 0.5, // saturation constant S_HCO [mole HCO3 m-3]
+ b_A_O: 0.15, // aerobic respiration rate [d-1]
+ b_A_NO: 0.05 // anoxic respiration rate [d-1]
+ };
+ this.stoi_params = {
+ // Stoichiometric and composition parameters
+
+ f_SI: 0., // fraction S_I from hydrolysis [g S_I g-1 X_S]
+ f_XI: 0.2, // fraction X_I from decomp X_H [g X_I g-1 X_H]
+ // Yields
+ Y_STO_O: 0.85, // aerobic yield X_STO per S_S [g X_STO g-1 S_S]
+ Y_STO_NO: 0.80, // anoxic yield X_STO per S_S [g X_STO g-1 S_S]
+ Y_H_O: 0.63, // aerobic yield X_H per X_STO [g X_H g-1 X_STO]
+ Y_H_NO: 0.54, // anoxic yield X_H per X_STO [g X_H g-1 X_STO]
+ Y_A: 0.24, // anoxic yield X_A per S_NO [g X_A g-1 NO3-N]
+ // Composition (COD via DoR)
+ i_CODN: -1.71, // COD content (DoR) [g COD g-1 N2-N]
+ i_CODNO: -4.57, // COD content (DoR) [g COD g-1 NO3-N]
+ // Composition (nitrogen)
+ i_NSI: 0.01, // nitrogen content S_I [g N g-1 S_I]
+ i_NSS: 0.03, // nitrogen content S_S [g N g-1 S_S]
+ i_NXI: 0.02, // nitrogen content X_I [g N g-1 X_I]
+ i_NXS: 0.04, // nitrogen content X_S [g N g-1 X_S]
+ i_NBM: 0.07, // nitrogen content X_H / X_A [g N g-1 X_H / X_A]
+ // Composition (TSS)
+ i_TSXI: 0.75, // TSS content X_I [g TS g-1 X_I]
+ i_TSXS: 0.75, // TSS content X_S [g TS g-1 X_S]
+ i_TSBM: 0.90, // TSS content X_H / X_A [g TS g-1 X_H / X_A]
+ i_TSSTO: 0.60, // TSS content X_STO (PHB based) [g TS g-1 X_STO]
+ // Composition (charge)
+ i_cNH: 1/14, // charge per S_NH [mole H+ g-1 NH3-N]
+ i_cNO: -1/14 // charge per S_NO [mole H+ g-1 NO3-N]
+ };
+ this.stoi_matrix = this._initialise_stoi_matrix();
}
_initialise_stoi_matrix() { // initialise stoichiometric matrix
diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index db8b3bb..4027690 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -1,11 +1,16 @@
const ASM3 = require('./asm3_class')
-const math = require('mathjs')
+const { create, all } = require('mathjs')
+
+const config = {
+ matrix: 'Array' // Choose 'Matrix' (default) or 'Array'
+}
+
+const math = create(all, config)
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]
@@ -16,7 +21,8 @@ class Reactor_CSTR {
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.timeStep = 1/(24*60*15); // time step [d]
+ this.speedUpFactor = 1;
}
set setInfluent(input) { // setter for C_in (WIP)
@@ -39,20 +45,18 @@ class Reactor_CSTR {
}
// expect update with timestamp
- updateState(timestamp) {
- let newTime = timestamp;
+ updateState(newTime) {
const day2ms = 1000 * 60 * 60 * 24;
- let n_iter = Math.floor((newTime - this.currentTime) / (this.timeStep * day2ms));
- if (n_iter > 0) {
+ let n_iter = Math.floor(this.speedUpFactor*(newTime - this.currentTime) / (this.timeStep * day2ms));
+ if (n_iter) {
let n = 0;
while (n < n_iter) {
- console.log(this.tick_fe(this.timeStep));
+ this.tick_fe(this.timeStep);
n += 1;
}
- this.currentTime += n_iter * this.timeStep * day2ms;
- n_iter = 0;
+ this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
}
}
@@ -70,15 +74,135 @@ 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())
+
+ // console.log("Initial State: ")
+ // console.log(this.state)
+
+ 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.1; // 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*60); // time step [d]
+ this.speedUpFactor = 1;
+
+ this.D_op = this.makeDoperator();
+ this.D2_op = this.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.at(-1)}, 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
+ 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((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
+ transfer.forEach((x) => { x[0] = this.OTR; });
+ } else {
+ transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); });
+ }
+
+ if (math.sum(this.Fs) > 0) { // Danckwerts BC
+ const BC_influx = math.multiply(math.divide([this.Fs], this.A), 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];
+ this.state[0] = math.add(BC_influx, BC_dispersion);
+ }
+
+ 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.resize(math.diag(Array(this.n_x).fill(-1), 1), [this.n_x, this.n_x]);
+ I[0][0] = 0;
+ I[0][1] = 1;
+ 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.identity(this.n_x);
+ 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;
+ I[0][1] = 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];
-// const Reactor = new Reactor_CSTR(initial_state);
-// Reactor.C_in = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
-// N = 0;
+// const Reactor = new Reactor_PFR(200, 10, 10, 1, 100, initial_state);
+// Reactor.Cs_in[0] = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
+// Reactor.Fs[0] = 10;
+// let N = 0;
// while (N < 500) {
// console.log(Reactor.tick_fe(0.001));
// N += 1;
// }
-module.exports = Reactor_CSTR;
\ No newline at end of file
+module.exports = {Reactor_CSTR, Reactor_PFR};
\ No newline at end of file