Refactor reactor node registration
This commit is contained in:
128
src/asm3_class.js
Normal file
128
src/asm3_class.js
Normal file
@@ -0,0 +1,128 @@
|
||||
const math = require('mathjs')
|
||||
|
||||
class ASM3 {
|
||||
|
||||
constructor() {
|
||||
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
|
||||
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;
|
||||
|
||||
const stoi_matrix = Array(12);
|
||||
// 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
|
||||
stoi_matrix[0] = [0., f_SI, 1.-f_SI, i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI, 0., 0., (i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI)*i_cNH, 0., -1., 0., 0., 0., -i_TSXS];
|
||||
stoi_matrix[1] = [-(1.-Y_STO_O), 0, -1., i_NSS, 0., 0., i_NSS*i_cNH, 0., 0., 0., Y_STO_O, 0., Y_STO_O*i_TSSTO];
|
||||
stoi_matrix[2] = [0., 0., -1., i_NSS, -(1.-Y_STO_NO)/(i_CODNO-i_CODN), (1.-Y_STO_NO)/(i_CODNO-i_CODN), i_NSS*i_cNH + (1.-Y_STO_NO)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., Y_STO_NO, 0., Y_STO_NO*i_TSSTO];
|
||||
stoi_matrix[3] = [-(1.-Y_H_O)/Y_H_O, 0., 0., -i_NBM, 0., 0., -i_NBM*i_cNH, 0., 0., 1., -1./Y_H_O, 0., i_TSBM-i_TSSTO/Y_H_O];
|
||||
stoi_matrix[4] = [0., 0., 0., -i_NBM, -(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), (1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), -i_NBM*i_cNH+(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN))*i_cNO, 0., 0., 1., -1./Y_H_NO, 0., i_TSBM-i_TSSTO/Y_H_NO];
|
||||
stoi_matrix[5] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
|
||||
stoi_matrix[6] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
|
||||
stoi_matrix[7] = [-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., -i_TSSTO];
|
||||
stoi_matrix[8] = [0., 0., 0., 0., -1./(i_CODNO-i_CODN), 1./(i_CODNO-i_CODN), i_cNO/(i_CODNO-i_CODN), 0., 0., 0., -1., 0., -i_TSSTO];
|
||||
stoi_matrix[9] = [1.+i_CODNO/Y_A, 0., 0., -1./Y_A-i_NBM, 0., 1./Y_A, (-1./Y_A-i_NBM)*i_cNH+i_cNO/Y_A, 0., 0., 0., 0., 1., i_TSBM];
|
||||
stoi_matrix[10] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
|
||||
stoi_matrix[11] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
|
||||
|
||||
return stoi_matrix[0].map((col, i) => stoi_matrix.map(row => row[i])); // transpose matrix
|
||||
}
|
||||
|
||||
_monod(c, K){
|
||||
return c / (K + c);
|
||||
}
|
||||
|
||||
_inv_monod(c, K){
|
||||
return K / (K + c);
|
||||
}
|
||||
|
||||
compute_rates(state) { // computes reaction rates. state is optional
|
||||
// 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;
|
||||
const { k_H, K_X, k_STO, nu_NO, K_O, K_NO, K_S, K_STO, mu_H_max, K_NH, K_HCO, b_H_O, b_H_NO, b_STO_O, b_STO_NO, mu_A_max, K_A_NH, K_A_O, K_A_HCO, b_A_O, b_A_NO } = this.kin_params;
|
||||
|
||||
// Hydrolysis
|
||||
rates[0] = X_H == 0 ? 0 : k_H * this._monod(X_S / X_H, K_X) * X_H;
|
||||
|
||||
// Heterotrophs
|
||||
rates[1] = k_STO * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
|
||||
rates[2] = k_STO * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_S, K_S) * X_H;
|
||||
rates[3] = X_H == 0 ? 0 : mu_H_max * this._monod(S_O, K_O) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
|
||||
rates[4] = X_H == 0 ? 0 : mu_H_max * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
|
||||
rates[5] = b_H_O * this._monod(S_O, K_O) * X_H;
|
||||
rates[6] = b_H_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
|
||||
rates[7] = b_STO_O * this._monod(S_O, K_O) * X_H;
|
||||
rates[8] = b_STO_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
|
||||
|
||||
// Autotrophs
|
||||
rates[9] = mu_A_max * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
|
||||
rates[10] = b_A_O * this._monod(S_O, K_O) * X_A;
|
||||
rates[11] = b_A_NO * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
|
||||
|
||||
return rates;
|
||||
}
|
||||
|
||||
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));
|
||||
}
|
||||
}
|
||||
|
||||
module.exports = ASM3;
|
||||
109
src/nodeClass.js
Normal file
109
src/nodeClass.js
Normal file
@@ -0,0 +1,109 @@
|
||||
const { Reactor_CSTR, Reactor_PFR } = require('./reactor_class.js');
|
||||
|
||||
|
||||
class nodeClass {
|
||||
/**
|
||||
* Create a ReactorNode.
|
||||
* @param {object} uiConfig - Node-RED node configuration.
|
||||
* @param {object} RED - Node-RED runtime API.
|
||||
* @param {object} nodeInstance - The Node-RED node instance.
|
||||
* @param {string} nameOfNode - The name of the node, used for
|
||||
*/
|
||||
constructor(uiConfig, RED, nodeInstance, nameOfNode) {
|
||||
// Preserve RED reference for HTTP endpoints if needed
|
||||
this.node = nodeInstance;
|
||||
this.RED = RED;
|
||||
this.name = nameOfNode;
|
||||
|
||||
let new_reactor;
|
||||
|
||||
switch (uiConfig.reactor_type) {
|
||||
case "CSTR":
|
||||
new_reactor = new Reactor_CSTR(
|
||||
parseFloat(uiConfig.volume),
|
||||
parseInt(uiConfig.n_inlets),
|
||||
parseFloat(uiConfig.kla),
|
||||
[
|
||||
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)
|
||||
]
|
||||
);
|
||||
break;
|
||||
case "PFR":
|
||||
new_reactor = new Reactor_PFR(
|
||||
parseFloat(uiConfig.volume),
|
||||
parseFloat(uiConfig.length),
|
||||
parseInt(uiConfig.resolution_L),
|
||||
parseInt(uiConfig.n_inlets),
|
||||
parseFloat(uiConfig.kla),
|
||||
[
|
||||
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)
|
||||
]
|
||||
);
|
||||
break;
|
||||
default:
|
||||
console.warn("Unknown reactor type: " + uiConfig.reactor_type);
|
||||
}
|
||||
|
||||
const reactor = new_reactor; // protect from reassignment
|
||||
|
||||
this.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();
|
||||
}
|
||||
});
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
module.exports = nodeClass;
|
||||
262
src/reactor_class.js
Normal file
262
src/reactor_class.js
Normal file
@@ -0,0 +1,262 @@
|
||||
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];
|
||||
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 };
|
||||
Reference in New Issue
Block a user