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

Merge pull request 'MVP for dispersion model' (#5) from experimental into main

Browse Source 
Reviewed-on: p.vanderwilt/asm3#5
This commit is contained in:
p.vanderwilt
2025-06-24 10:36:09 +00:00
parent 3d3304f23a 6b57a46aab
commit c279552d7e
4 changed files with 301 additions and 96 deletions
Download Patch File Download Diff File Expand all files Collapse all files

51
advanced-reactor.html
View File

@@ -4,7 +4,10 @@
color: "#c4cce0",
defaults: {
name: { value: "" },
volume: { value: 0., required: true},
reactor_type: { value: "CSTR", required: true },
volume: { value: 0., required: true },
length: { value: 0.},
resolution_L: { value: 0.},
n_inlets: { value: 1, required: true},
kla: { value: null },
S_O_init: { value: 0., required: true },
@@ -37,6 +40,14 @@
type:"num",
types:["num"]
});
$("#node-input-length").typedInput({
type:"num",
types:["num"]
});
$("#node-input-resolution_L").typedInput({
type:"num",
types:["num"]
});
$("#node-input-kla").typedInput({
type:"num",
types:["num"]
@@ -45,6 +56,32 @@
type:"num",
types:["num"]
});
$("#node-input-reactor_type").typedInput({
types: [
{
value: "CSTR",
options: [
{ value: "CSTR", label: "CSTR"},
{ value: "PFR", label: "PFR"}
]
}
]
})
$("#node-input-reactor_type").on("change", function() {
const type = $("#node-input-reactor_type").typedInput("value");
if (type === "CSTR") {
$(".PFR").hide();
} else {
$(".PFR").show();
}
});
// Set initial visibility on dialog open
const initialType = $("#node-input-reactor_type").typedInput("value");
if (initialType === "CSTR") {
$(".PFR").hide();
} else {
$(".PFR").show();
}
},
oneditsave: function() {
let volume = parseFloat($("#node-input-volume").typedInput("value"));
@@ -65,10 +102,22 @@
<input type="text" id="node-input-name" placeholder="Name">
</div>
<h2> Reactor properties </h2>
<div class="form-row">
<label for="node-input-reactor_type"><i class="fa fa-tag"></i> Reactor type</label>
<input type="text" id="node-input-reactor_type">
</div>
<div class="form-row">
<label for="node-input-volume"><i class="fa fa-tag"></i> Fluid volume [m3]</label>
<input type="text" id="node-input-volume" placeholder="m3">
</div>
<div class="form-row PFR">
<label for="node-input-length"><i class="fa fa-tag"></i> Reactor length [m]</label>
<input type="text" id="node-input-length" placeholder="m">
</div>
<div class="form-row PFR">
<label for="node-input-resolution_L"><i class="fa fa-tag"></i> Resolution</label>
<input type="text" id="node-input-resolution_L" placeholder="#">
</div>
<div class="form-row">
<label for="node-input-n_inlets"><i class="fa fa-tag"></i> Number of inlets</label>
<input type="text" id="node-input-n_inlets" placeholder="#">

76
advanced-reactor.js
View File

@@ -5,28 +5,62 @@ module.exports = function(RED) {
let name = config.name;
const Reactor = require('./dependencies/reactor_class');
const { Reactor_CSTR, Reactor_PFR } = require('./dependencies/reactor_class');
const reactor = new Reactor(
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)
]
);
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;

118
dependencies/asm3_class.js vendored
View File

@@ -2,67 +2,65 @@ const math = require('mathjs')
class ASM3 {
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]
}
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]
}
constructor() {
this.stoi_matrix = this._initialise_stoi_matrix()
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

152
dependencies/reactor_class.js vendored
View File

@@ -1,11 +1,16 @@
const ASM3 = require('./asm3_class')
const math = require('mathjs')
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;
console.log(this.state);
this.asm = new ASM3();
this.Vl = volume; // fluid volume reactor [m3]
@@ -16,7 +21,8 @@ class Reactor_CSTR {
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.timeStep = 1/(24*60*15); // time step [d]
this.speedUpFactor = 1;
}
set setInfluent(input) { // setter for C_in (WIP)
@@ -39,20 +45,18 @@ class Reactor_CSTR {
}
// expect update with timestamp
updateState(timestamp) {
let newTime = timestamp;
updateState(newTime) {
const day2ms = 1000 * 60 * 60 * 24;
let n_iter = Math.floor((newTime - this.currentTime) / (this.timeStep * day2ms));
if (n_iter > 0) {
let n_iter = Math.floor(this.speedUpFactor*(newTime - this.currentTime) / (this.timeStep * day2ms));
if (n_iter) {
let n = 0;
while (n < n_iter) {
console.log(this.tick_fe(this.timeStep));
this.tick_fe(this.timeStep);
n += 1;
}
this.currentTime += n_iter * this.timeStep * day2ms;
n_iter = 0;
this.currentTime += n_iter * this.timeStep * day2ms / this.speedUpFactor;
}
}
@@ -70,15 +74,135 @@ class Reactor_CSTR {
}
}
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.1; // 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*60); // time step [d]
this.speedUpFactor = 1;
this.D_op = this.makeDoperator();
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;
}
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(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));
reaction[0] = Array(13).fill(0.0);
const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
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]); });
}
if (math.sum(this.Fs) > 0) { // Danckwerts BC
const BC_influx = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0];
const BC_gradient = Array(this.n_x).fill(0.0);
BC_gradient[0] = 1;
BC_gradient[1] = -1;
const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], this.state)[0];
this.state[0] = math.add(BC_influx, BC_dispersion);
}
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
return this.state;
}
makeDoperator() { // create the upwind scheme gradient operator
const I = math.identity(this.n_x);
const A = math.resize(math.diag(Array(this.n_x).fill(-1), 1), [this.n_x, this.n_x]);
I[0][0] = 0;
I[0][1] = 1;
I[this.n_x-1][this.n_x-1] = 0; // Neumann boundary condition at x=L
return math.add(I, A);
}
makeD2operator() { // create the upwind scheme second derivative operator
const I = math.identity(this.n_x);
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]);
I[0][0] = 0;
I[0][1] = 1;
return math.add(I, A, B);
}
}
// 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_CSTR(initial_state);
// Reactor.C_in = [0.0, 30., 100., 16., 0., 0., 5., 25., 75., 30., 0., 0., 125.];
// N = 0;
// 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;
// let N = 0;
// while (N < 500) {
// console.log(Reactor.tick_fe(0.001));
// N += 1;
// }
module.exports = Reactor_CSTR;
module.exports = {Reactor_CSTR, Reactor_PFR};