diff --git a/advanced-reactor.js b/advanced-reactor.js index e2eca51..9ca6959 100644 --- a/advanced-reactor.js +++ b/advanced-reactor.js @@ -1,99 +1,12 @@ -module.exports = function(RED) { - function reactor(config) { +const nameOfNode = "advanced-reactor"; // name of the node, should match file name and node type in Node-RED +const nodeClass = require('./src/nodeClass.js'); // node class + +module.exports = function (RED) { + // Register the node type + RED.nodes.registerType(nameOfNode, function (config) { + // Initialize the Node-RED node first RED.nodes.createNode(this, config); - var node = this; - - let name = config.name; - - const { Reactor_CSTR, Reactor_PFR } = require('./dependencies/reactor_class'); - - 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; - - switch (msg.topic) { - case "clock": - toggleUpdate = true; - break; - case "Fluent": - reactor.setInfluent = msg; - if (msg.payload.inlet == 0) { - toggleUpdate = true; - } - break; - case "OTR": - reactor.setOTR = msg; - break; - case "Dispersion": - reactor.setDispersion = msg; - break; - default: - console.log("Unknown topic: " + msg.topic); - } - - if (toggleUpdate) { - reactor.updateState(msg.timestamp); - send(reactor.getEffluent); - } - - if (done) { - done(); - } - }); - } - RED.nodes.registerType("advanced-reactor", reactor); + // Then create your custom class and attach it + this.nodeClass = new nodeClass(config, RED, this, nameOfNode); + }); }; diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js deleted file mode 100644 index 76a6340..0000000 --- a/dependencies/reactor_class.js +++ /dev/null @@ -1,263 +0,0 @@ -const ASM3 = require('./asm3_class') -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; - this.asm = new ASM3(); - - this.Vl = volume; // fluid volume reactor [m3] - 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.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; - } - - 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; - } - - 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 - const r = this.asm.compute_dC(this.state); - const dC_in = math.multiply(math.divide([this.Fs], this.Vl), this.Cs_in)[0]; - const dC_out = math.multiply(-1*math.sum(this.Fs)/this.Vl, this.state); - const t_O = Array(13).fill(0.0); - t_O[0] = isNaN(this.kla) ? this.OTR : this.calcOTR(this.state[0]); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR - - const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step); - - // clip value element-wise to each subarray to avoid negative concentrations - this.state = math.add(this.state, dC_total).map(val => val < 0 ? 0 : val); - return this.state; - } -} - -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.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 = 60; - - this.D_op = this.makeDoperator(true, true); - 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; - console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A)); - console.log("Pe local " + this.d_x*math.sum(this.Fs)/(this.D*this.A)); - console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x)); - console.log("Co D " + this.D*this.timeStep/(this.d_x*this.d_x)); - } - - 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(-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)); - const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0)); - - if (dispersion.some(row => row.some(Number.isNaN))) { - throw new Error("NaN detected in dispersion!"); - } - if (advection.some(row => row.some(Number.isNaN))) { - throw new Error("NaN detected in advection!"); - } - if (reaction.some(row => row.some(Number.isNaN))) { - throw new Error("NaN detected in reaction!"); - } - - 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); - const new_state = math.add(this.state, dC_total); - if (new_state.some(row => row.some(Number.isNaN))) { - throw new Error("NaN detected in new_state after dC_total update!"); - } - - // 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_gradient = Array(this.n_x).fill(0.0); - BC_gradient[0] = -1; - BC_gradient[1] = 1; - let Pe = this.length*math.sum(this.Fs)/(this.D*this.A) - const BC_dispersion = math.multiply((1-(1+4*this.volume/math.sum(this.Fs)/Pe)^0.5)/Pe, [BC_gradient], new_state)[0]; - console.log(math.add(BC_C_in, BC_dispersion)); - new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val); - - } 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] - - if (new_state.some(row => row.some(Number.isNaN))) { - throw new Error("NaN detected in new_state after enforcing boundary conditions!"); - } - - this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state - return new_state; - } - - makeDoperator(central=false, higher_order=false) { // create gradient operator - if (higher_order) { - if (central) { - const I = math.resize(math.diag(Array(this.n_x).fill(1/12), -2), [this.n_x, this.n_x]); - const A = math.resize(math.diag(Array(this.n_x).fill(-2/3), -1), [this.n_x, this.n_x]); - const B = math.resize(math.diag(Array(this.n_x).fill(2/3), 1), [this.n_x, this.n_x]); - const C = math.resize(math.diag(Array(this.n_x).fill(-1/12), 2), [this.n_x, this.n_x]); - const D = math.add(I, A, B, C); - const NearBoundary = Array(this.n_x).fill(0.0); - NearBoundary[0] = -1/4; - NearBoundary[1] = -5/6; - NearBoundary[2] = 3/2; - NearBoundary[3] = -1/2; - NearBoundary[4] = 1/12; - D[1] = NearBoundary; - NearBoundary.reverse(); - D[this.n_x-2] = math.multiply(-1, NearBoundary); - D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere - D[this.n_x-1] = Array(this.n_x).fill(0); - return D; - } else { - throw new Error("Upwind higher order method not implemented! Use central scheme instead."); - } - } else { - const I = math.resize(math.diag(Array(this.n_x).fill(1/(1+central)), central), [this.n_x, this.n_x]); - const A = math.resize(math.diag(Array(this.n_x).fill(-1/(1+central)), -1), [this.n_x, this.n_x]); - const D = math.add(I, A); - D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere - 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]); - const D2 = math.add(I, A, B); - D2[0] = Array(this.n_x).fill(0); // set by BCs elsewhere - D2[this.n_x-1] = Array(this.n_x).fill(0); - return D2; - } -} - - -// 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_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; -// Reactor.D = 0.01; -// let N = 0; -// while (N < 5000) { -// console.log(Reactor.tick_fe(0.001)); -// N += 1; -// } - -module.exports = { Reactor_CSTR, Reactor_PFR }; \ No newline at end of file diff --git a/src/nodeClass.js b/src/nodeClass.js new file mode 100644 index 0000000..7919df2 --- /dev/null +++ b/src/nodeClass.js @@ -0,0 +1,114 @@ +const { Reactor_CSTR, Reactor_PFR } = require('./specificClass.js'); + + +class nodeClass { + /** + * Node-RED node class for advanced-reactor. + * @param {object} uiConfig - Node-RED node configuration + * @param {object} RED - Node-RED runtime API + * @param {object} nodeInstance - Node-RED node instance + * @param {string} nameOfNode - Name of the node + */ + constructor(uiConfig, RED, nodeInstance, nameOfNode) { + // Preserve RED reference for HTTP endpoints if needed + this.node = nodeInstance; + this.RED = RED; + this.name = nameOfNode; + + this._loadConfig(uiConfig) + + this._setupClass(); + + this._attachInputHandler(); + } + + /** + * Handle node-red input messages + */ + _attachInputHandler() { + this.node.on('input', (msg, send, done) => { + let toggleUpdate = false; + + switch (msg.topic) { + case "clock": + toggleUpdate = true; + break; + case "Fluent": + this.reactor.setInfluent = msg; + if (msg.payload.inlet == 0) { + toggleUpdate = true; + } + break; + case "OTR": + this.reactor.setOTR = msg; + break; + case "Dispersion": + this.reactor.setDispersion = msg; + break; + default: + console.log("Unknown topic: " + msg.topic); + } + + if (toggleUpdate) { + this.reactor.updateState(msg.timestamp); + send(this.reactor.getEffluent); + } + + if (done) { + done(); + } + }); + } + + /** + * Parse node configuration + * @param {object} uiConfig Config set in UI in node-red + */ + _loadConfig(uiConfig) { + this.config = { + reactor_type: uiConfig.reactor_type, + volume: parseFloat(uiConfig.volume), + length: parseFloat(uiConfig.length), + resolution_L: parseInt(uiConfig.resolution_L), + n_inlets: parseInt(uiConfig.n_inlets), + kla: parseFloat(uiConfig.kla), + initialState: [ + parseFloat(uiConfig.S_O_init), + parseFloat(uiConfig.S_I_init), + parseFloat(uiConfig.S_S_init), + parseFloat(uiConfig.S_NH_init), + parseFloat(uiConfig.S_N2_init), + parseFloat(uiConfig.S_NO_init), + parseFloat(uiConfig.S_HCO_init), + parseFloat(uiConfig.X_I_init), + parseFloat(uiConfig.X_S_init), + parseFloat(uiConfig.X_H_init), + parseFloat(uiConfig.X_STO_init), + parseFloat(uiConfig.X_A_init), + parseFloat(uiConfig.X_TS_init) + ] + } + } + + /** + * Setup reactor class based on config + */ + _setupClass() { + let new_reactor; + + switch (this.config.reactor_type) { + case "CSTR": + new_reactor = new Reactor_CSTR(this.config); + break; + case "PFR": + new_reactor = new Reactor_PFR(this.config); + break; + default: + console.warn("Unknown reactor type: " + uiConfig.reactor_type); + } + + this.reactor = new_reactor; // protect from reassignment + } +} + +module.exports = nodeClass; \ No newline at end of file diff --git a/dependencies/asm3_class.js b/src/reaction_modules/asm3_class.js similarity index 81% rename from dependencies/asm3_class.js rename to src/reaction_modules/asm3_class.js index f5ac1df..ab272dc 100644 --- a/dependencies/asm3_class.js +++ b/src/reaction_modules/asm3_class.js @@ -1,11 +1,16 @@ const math = require('mathjs') +/** + * ASM3 class for the Activated Sludge Model No. 3 (ASM3). + */ class ASM3 { constructor() { + /** + * Kinetic parameters for ASM3 at 20 C. + * @property {Object} kin_params - Kinetic parameters + */ 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] @@ -31,9 +36,13 @@ class ASM3 { 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 + /** + * Stoichiometric and composition parameters for ASM3. + * @property {Object} stoi_params - Stoichiometric parameters + */ + this.stoi_params = { + // Fractions 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 @@ -63,6 +72,10 @@ class ASM3 { this.stoi_matrix = this._initialise_stoi_matrix(); } + /** + * Initialises the stoichiometric matrix for ASM3. + * @returns {Array} - The stoichiometric matrix for ASM3. (2D array) + */ _initialise_stoi_matrix() { // initialise stoichiometric matrix const { f_SI, f_XI, Y_STO_O, Y_STO_NO, Y_H_O, Y_H_NO, Y_A, i_CODN, i_CODNO, i_NSI, i_NSS, i_NXI, i_NXS, i_NBM, i_TSXI, i_TSXS, i_TSBM, i_TSSTO, i_cNH, i_cNO } = this.stoi_params; @@ -84,15 +97,32 @@ class ASM3 { return stoi_matrix[0].map((col, i) => stoi_matrix.map(row => row[i])); // transpose matrix } + /** + * Computes the Monod equation rate value for a given concentration and half-saturation constant. + * @param {number} c - Concentration of reaction species. + * @param {number} K - Half-saturation constant for the reaction species. + * @returns {number} - Monod equation rate value for the given concentration and half-saturation constant. + */ _monod(c, K){ return c / (K + c); } + /** + * Computes the inverse Monod equation rate value for a given concentration and half-saturation constant. Used for inhibition. + * @param {number} c - Concentration of reaction species. + * @param {number} K - Half-saturation constant for the reaction species. + * @returns {number} - Inverse Monod equation rate value for the given concentration and half-saturation constant. + */ _inv_monod(c, K){ return K / (K + c); } - compute_rates(state) { // computes reaction rates. state is optional + /** + * Computes the reaction rates for each process reaction based on the current state. + * @param {Array} state - State vector containing concentrations of reaction species. + * @returns {Array} - Reaction rates for each process reaction. + */ + compute_rates(state) { // 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 const rates = Array(12); const [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] = state; @@ -119,6 +149,11 @@ class ASM3 { return rates; } + /** + * Computes the change in concentrations of reaction species based on the current state. + * @param {Array} state - State vector containing concentrations of reaction species. + * @returns {Array} - Change in reaction species concentrations. + */ compute_dC(state) { // compute changes in concentrations // 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 return math.multiply(this.stoi_matrix, this.compute_rates(state)); diff --git a/src/specificClass.js b/src/specificClass.js new file mode 100644 index 0000000..384aaf0 --- /dev/null +++ b/src/specificClass.js @@ -0,0 +1,320 @@ +const ASM3 = require('./reaction_modules/asm3_class.js'); +const { create, all } = require('mathjs'); +const { assertNoNaN } = require('./utils.js'); + +const config = { + matrix: 'Array' // use Array as the matrix type +}; + +const math = create(all, config); + +const S_O_INDEX = 0; +const NUM_SPECIES = 13; +const DEBUG = false; + +class Reactor { + /** + * Reactor base class. + * @param {object} config - Configuration object containing reactor parameters. + */ + constructor(config) { + this.asm = new ASM3(); + + this.volume = config.volume; // fluid volume reactor [m3] + + this.Fs = Array(config.n_inlets).fill(0); // fluid debits per inlet [m3 d-1] + this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0)); // composition influents + this.OTR = 0.0; // oxygen transfer rate [g O2 d-1] + + this.kla = config.kla; // if NaN, use externaly provided OTR [d-1] + + this.currentTime = Date.now(); // milliseconds since epoch [ms] + this.timeStep = 1 / (24*60*15); // time step [d] + this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second + } + + /** + * Setter for influent data. + * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations. + */ + set setInfluent(input) { + let index_in = input.payload.inlet; + this.Fs[index_in] = input.payload.F; + this.Cs_in[index_in] = input.payload.C; + } + + /** + * Setter for OTR (Oxygen Transfer Rate). + * @param {object} input - Input object (msg) containing payload with OTR value [g O2 d-1]. + */ + set setOTR(input) { + this.OTR = input.payload; + } + + /** + * Calculate the oxygen transfer rate (OTR) based on the dissolved oxygen concentration and temperature. + * @param {number} S_O - Dissolved oxygen concentration [g O2 m-3]. + * @param {number} T - Temperature in Celsius, default to 20 C. + * @returns {number} - Calculated OTR [g O2 d-1]. + */ + _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); + } + + /** + * Clip values in an array to zero. + * @param {Array} arr - Array of values to clip. + * @returns {Array} - New array with values clipped to zero. + */ + _arrayClip2Zero(arr) { + if (Array.isArray(arr)) { + return arr.map(x => this._arrayClip2Zero(x)); + } else { + return arr < 0 ? 0 : arr; + } + } + + /** + * Update the reactor state based on the new time. + * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch. + */ + updateState(newTime) { // expect update with timestamp + 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(this.timeStep); + n += 1; + } + this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor; + } + } + +} + +class Reactor_CSTR extends Reactor { + /** + * Reactor_CSTR class for Continuous Stirred Tank Reactor. + * @param {object} config - Configuration object containing reactor parameters. + */ + constructor(config) { + super(config); + this.state = config.initialState; + } + + /** + * Getter for effluent data. + * @returns {object} Effluent data object (msg), defaults to inlet 0. + */ + 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 }; + } + + /** + * Tick the reactor state using the forward Euler method. + * @param {number} time_step - Time step for the simulation [d]. + * @returns {Array} - New reactor state. + */ + tick(time_step) { // tick reactor state using forward Euler method + const inflow = math.multiply(math.divide([this.Fs], this.volume), this.Cs_in)[0]; + const outflow = math.multiply(-1 * math.sum(this.Fs) / this.volume, this.state); + const reaction = this.asm.compute_dC(this.state); + const transfer = Array(NUM_SPECIES).fill(0.0); + transfer[S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[S_O_INDEX]); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR + + const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step) + this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations + if(DEBUG){ + assertNoNaN(dC_total, "change in state"); + assertNoNaN(this.state, "new state"); + } + return this.state; + } +} + +class Reactor_PFR extends Reactor { + /** + * Reactor_PFR class for Plug Flow Reactor. + * @param {object} config - Configuration object containing reactor parameters. + */ + constructor(config) { + super(config); + + this.length = config.length; // reactor length [m] + this.n_x = config.resolution_L; // number of slices + + this.d_x = this.length / this.n_x; + this.A = this.volume / this.length; // crosssectional area [m2] + + this.state = Array.from(Array(this.n_x), () => config.initialState.slice()) + + // console.log("Initial State: ") + // console.log(this.state) + + this.D = 0.0; // axial dispersion [m2 d-1] + + this.D_op = this._makeDoperator(true, true); + assertNoNaN(this.D_op, "Derivative operator"); + + this.D2_op = this._makeD2operator(); + assertNoNaN(this.D2_op, "Second derivative operator"); + } + + /** + * Setter for axial dispersion. + * @param {object} input - Input object (msg) containing payload with dispersion value [m2 d-1]. + */ + set setDispersion(input) { + this.D = input.payload; + } + + /** + * Setter for influent data. + * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations. + */ + set setInfluent(input) { + super.setInfluent = input; + if(DEBUG) { + console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A)); + console.log("Pe local " + this.d_x*math.sum(this.Fs)/(this.D*this.A)); + console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x)); + console.log("Co D " + this.D*this.timeStep/(this.d_x*this.d_x)); + } + } + + /** + * Getter for effluent data. + * @returns {object} Effluent data object (msg), defaults to inlet 0. + */ + get getEffluent() { + return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime }; + } + + /** + * Apply boundary conditions to the reactor state. + * for inlet, apply generalised Danckwerts BC, if there is not flow, apply Neumann BC with no flux + * for outlet, apply regular Danckwerts BC (Neumann BC with no flux) + * @param {Array} state - Current reactor state without enforced BCs. + */ + _applyBoundaryConditions(state) { + 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_gradient = Array(this.n_x).fill(0); + BC_gradient[0] = -1; + BC_gradient[1] = 1; + let Pe = this.length * math.sum(this.Fs) / (this.D * this.A); + let residence_time = this.volume/math.sum(this.Fs); + const BC_dispersion = math.multiply((1 - (1 + 4*residence_time/Pe)^0.5) / (Pe*this.d_x), [BC_gradient], state)[0]; + state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val); + } else { // Neumann BC (no flux) + state[0] = state[1]; + } + // Neumann BC (no flux) + state[this.n_x-1] = state[this.n_x-2] + } + + /** + * Tick the reactor state using explicit finite difference method. + * @param {number} time_step - Time step for the simulation [d]. + * @returns {Array} - New reactor state. + */ + tick(time_step) { + 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)); + const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).fill(0)); + + if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR + transfer.forEach((x) => { x[S_O_INDEX] = this.OTR; }); + } else { + transfer.forEach((x, i) => { x[S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX]); }); + } + + const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); + + const stateNew = math.add(this.state, dC_total); + this._applyBoundaryConditions(stateNew); + + if (DEBUG) { + assertNoNaN(dispersion, "dispersion"); + assertNoNaN(advection, "advection"); + assertNoNaN(reaction, "reaction"); + assertNoNaN(dC_total, "change in state"); + assertNoNaN(stateNew, "new state post BC"); + } + + this.state = this._arrayClip2Zero(stateNew); + return stateNew; + } + + /** + * Create finite difference first derivative operator. + * @param {boolean} central - Use central difference scheme if true, otherwise use upwind scheme. + * @param {boolean} higher_order - Use higher order scheme if true, otherwise use first order scheme. + * @returns {Array} - First derivative operator matrix. + */ + _makeDoperator(central = false, higher_order = false) { // create gradient operator + if (higher_order) { + if (central) { + const I = math.resize(math.diag(Array(this.n_x).fill(1/12), -2), [this.n_x, this.n_x]); + const A = math.resize(math.diag(Array(this.n_x).fill(-2/3), -1), [this.n_x, this.n_x]); + const B = math.resize(math.diag(Array(this.n_x).fill(2/3), 1), [this.n_x, this.n_x]); + const C = math.resize(math.diag(Array(this.n_x).fill(-1/12), 2), [this.n_x, this.n_x]); + const D = math.add(I, A, B, C); + const NearBoundary = Array(this.n_x).fill(0.0); + NearBoundary[0] = -1/4; + NearBoundary[1] = -5/6; + NearBoundary[2] = 3/2; + NearBoundary[3] = -1/2; + NearBoundary[4] = 1/12; + D[1] = NearBoundary; + NearBoundary.reverse(); + D[this.n_x-2] = math.multiply(-1, NearBoundary); + D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere + D[this.n_x-1] = Array(this.n_x).fill(0); + return D; + } else { + throw new Error("Upwind higher order method not implemented! Use central scheme instead."); + } + } else { + const I = math.resize(math.diag(Array(this.n_x).fill(1 / (1+central)), central), [this.n_x, this.n_x]); + const A = math.resize(math.diag(Array(this.n_x).fill(-1 / (1+central)), -1), [this.n_x, this.n_x]); + const D = math.add(I, A); + D[0] = Array(this.n_x).fill(0); // set by BCs elsewhere + D[this.n_x-1] = Array(this.n_x).fill(0); + return D; + } + } + + /** + * Create central finite difference second derivative operator. + * @returns {Array} - Second derivative operator matrix. + */ + _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]); + const D2 = math.add(I, A, B); + D2[0] = Array(this.n_x).fill(0); // set by BCs elsewhere + D2[this.n_x - 1] = Array(this.n_x).fill(0); + return D2; + } +} + +module.exports = { Reactor_CSTR, Reactor_PFR }; + +// DEBUG +// 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_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; +// Reactor.D = 0.01; +// let N = 0; +// while (N < 5000) { +// console.log(Reactor.tick(0.001)); +// N += 1; +// } \ No newline at end of file diff --git a/src/utils.js b/src/utils.js new file mode 100644 index 0000000..2ca1896 --- /dev/null +++ b/src/utils.js @@ -0,0 +1,18 @@ +/** + * Assert that no NaN values are present in an array. + * @param {Array} arr + * @param {string} label + */ +function assertNoNaN(arr, label = "array") { + if (Array.isArray(arr)) { + for (const el of arr) { + assertNoNaN(el, label); + } + } else { + if (Number.isNaN(arr)) { + throw new Error(`NaN detected in ${label}!`); + } + } +} + +module.exports = { assertNoNaN }; \ No newline at end of file