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

Refactor advanced-reactor and nodeClass for improved readability and consistency

Browse Source
This commit is contained in:
p.vanderwilt
2025-07-04 15:14:03 +02:00
parent 3f5b0eea32
commit c6b0cab067
3 changed files with 64 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
advanced-reactor.js
View File

@@ -1,9 +1,9 @@
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) {
module.exports = function (RED) {
// Register the node type
RED.nodes.registerType(nameOfNode, function(config) {
RED.nodes.registerType(nameOfNode, function (config) {
// Initialize the Node-RED node first
RED.nodes.createNode(this, config);
// Then create your custom class and attach it

1
src/nodeClass.js
View File

@@ -20,7 +20,6 @@ class nodeClass {
this._setupClass();
this._attachInputHandler();
}
/**

107
src/reactor_class.js
View File

@@ -2,12 +2,17 @@ const ASM3 = require('./asm3_class.js')
const { create, all } = require('mathjs')
const config = {
matrix: 'Array' // choose 'Matrix' (default) or 'Array'
matrix: 'Array' // use Array as the matrix type
}
const math = create(all, config)
function assertNoNaN(arr, label="array") {
/**
* Assert that no NaN values are present in an array.
* @param {Array} arr
* @param {string} label
*/
function assertNoNaN(arr, label = "array") {
if (math.isNaN(arr)) {
throw new Error("NaN detected in ${label}!");
}
@@ -18,7 +23,7 @@ class Reactor {
* Reactor base class.
* @param {object} config - Configuration object containing reactor parameters.
*/
constructor(config){
constructor(config) {
this.asm = new ASM3();
this.Vl = config.volume; // fluid volume reactor [m3]
@@ -30,7 +35,7 @@ class Reactor {
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.timeStep = 1 / (24 * 60 * 15); // time step [d]
this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second
}
@@ -63,8 +68,8 @@ class Reactor {
* @param {number} T - Temperature in Celsius, default to 20 C.
* @returns
*/
_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;
_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);
}
@@ -88,7 +93,7 @@ class Reactor {
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));
let n_iter = Math.floor(this.speedUpFactor * (newTime - this.currentTime) / (this.timeStep * day2ms));
if (n_iter) {
let n = 0;
while (n < n_iter) {
@@ -116,7 +121,7 @@ class Reactor_CSTR extends Reactor {
* @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};
return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state }, timestamp: this.currentTime };
}
/**
@@ -127,13 +132,13 @@ class Reactor_CSTR extends Reactor {
tick(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 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);
this.state = this_.arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations
this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations
return this.state;
}
}
@@ -176,7 +181,25 @@ class Reactor_PFR extends Reactor {
* @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};
return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime };
}
_applyBoundaryConditions(newState) {
// 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], stateNew)[0];
newState[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
} else { // Neumann BC (no flux)
newState[0] = newState[1];
}
// Neumann BC (no flux)
newState[this.n_x - 1] = newState[this.n_x - 2]
return newState
}
/**
@@ -185,8 +208,8 @@ class Reactor_PFR extends Reactor {
* @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 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));
@@ -201,30 +224,16 @@ class Reactor_PFR extends Reactor {
}
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
const new_state = math.add(this.state, dC_total);
let stateNew = math.add(this.state, dC_total);
assertNoNaN(new_state, "new state");
assertNoNaN(stateNew, "new state");
// 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];
new_state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
stateNew = this._applyBoundaryConditions(stateNew);
} 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]
assertNoNaN(stateNew, "new state post BC");
assertNoNaN(new_state, "new state post BC");
this.state = this._arrayClip2Zero(new_state);
return new_state;
this.state = this._arrayClip2Zero(stateNew);
return stateNew;
}
/**
@@ -233,35 +242,35 @@ class Reactor_PFR extends Reactor {
* @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
_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 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;
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[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);
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 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);
D[this.n_x - 1] = Array(this.n_x).fill(0);
return D;
}
}
@@ -276,7 +285,7 @@ class Reactor_PFR extends Reactor {
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);
D2[this.n_x - 1] = Array(this.n_x).fill(0);
return D2;
}
}