RnD/reactor
3
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request 'Refactor code to align with other projects in EVOLV' (#7) from refactor into main

Browse Source 
Reviewed-on: p.vanderwilt/asm3#7
This commit is contained in:
p.vanderwilt
2025-07-07 10:26:46 +00:00
parent 0d12192ccc 302780726a
commit c89d5b0024
6 changed files with 502 additions and 365 deletions
Download Patch File Download Diff File Expand all files Collapse all files

107
advanced-reactor.js
View File

@@ -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);
});
};

263
dependencies/reactor_class.js vendored
View File

@@ -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 };

114
src/nodeClass.js Normal file
View File

@@ -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;

45
dependencies/asm3_class.js → src/reaction_modules/asm3_class.js
View File

@@ -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));

320
src/specificClass.js Normal file
View File

@@ -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;
// }

18
src/utils.js Normal file
View File

@@ -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 };