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

Enhance reactor class with NaN checks and refactor methods for clarity

Browse Source
This commit is contained in:
p.vanderwilt
2025-07-04 14:58:38 +02:00
parent 09e7072d16
commit 3f5b0eea32
1 changed files with 37 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

66
src/reactor_class.js
View File

@@ -1,4 +1,4 @@
const ASM3 = require('./asm3_class')
const ASM3 = require('./asm3_class.js')
const { create, all } = require('mathjs')
const config = {
@@ -7,6 +7,12 @@ const config = {
const math = create(all, config)
function assertNoNaN(arr, label="array") {
if (math.isNaN(arr)) {
throw new Error("NaN detected in ${label}!");
}
}
class Reactor {
/**
* Reactor base class.
@@ -57,11 +63,24 @@ 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
_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(clipToZero);
} 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.
@@ -110,12 +129,11 @@ class Reactor_CSTR extends Reactor {
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
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);
this.state = this_.arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations
return this.state;
}
}
@@ -141,8 +159,8 @@ class Reactor_PFR extends Reactor {
this.D = 0.0; // axial dispersion [m2 d-1]
this.D_op = this.makeDoperator(true, true);
this.D2_op = this.makeD2operator();
this.D_op = this._makeDoperator(true, true);
this.D2_op = this._makeD2operator();
}
/**
@@ -172,27 +190,20 @@ class Reactor_PFR extends Reactor {
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!");
}
assertNoNaN(dispersion, "dispersion");
assertNoNaN(advection, "advection");
assertNoNaN(reaction, "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]); });
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!");
}
assertNoNaN(new_state, "new state");
// apply boundary conditions
if (math.sum(this.Fs) > 0) { // Danckwerts BC
@@ -210,11 +221,9 @@ class Reactor_PFR extends Reactor {
// 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!");
}
assertNoNaN(new_state, "new state post BC");
this.state = new_state.map(row => row.map(val => val < 0 ? 0 : val)); // apply the new state
this.state = this._arrayClip2Zero(new_state);
return new_state;
}
@@ -224,7 +233,7 @@ 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]);
@@ -261,7 +270,7 @@ class Reactor_PFR extends Reactor {
* Create central finite difference second derivative operator.
* @returns {Array} - Second derivative operator matrix.
*/
makeD2operator() { // create the central second derivative operator
_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]);
@@ -272,6 +281,7 @@ class Reactor_PFR extends Reactor {
}
}
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
@@ -284,6 +294,4 @@ class Reactor_PFR extends Reactor {
// while (N < 5000) {
// console.log(Reactor.tick(0.001));
// N += 1;
// }
module.exports = { Reactor_CSTR, Reactor_PFR };
// }