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

Merge pull request 'New axial dispersion model with Generalised boundary conditions' (#6) from experimental into main

Browse Source 
Reviewed-on: p.vanderwilt/asm3#6
This commit is contained in:
p.vanderwilt
2025-07-03 20:30:34 +00:00
parent c279552d7e f517b7764d
commit 0d12192ccc
3 changed files with 96 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
advanced-reactor.js
View File

@@ -78,6 +78,9 @@ module.exports = function(RED) {
case "OTR":
reactor.setOTR = msg;
break;
case "Dispersion":
reactor.setDispersion = msg;
break;
default:
console.log("Unknown topic: " + msg.topic);
}

6
dependencies/asm3_class.js vendored
View File

@@ -99,13 +99,13 @@ class ASM3 {
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] = k_H * this._monod(X_S / X_H, K_X) * X_H;
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] = 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] = 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[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;

125
dependencies/reactor_class.js vendored
View File

@@ -2,7 +2,7 @@ const ASM3 = require('./asm3_class')
const { create, all } = require('mathjs')
const config = {
matrix: 'Array' // Choose 'Matrix' (default) or 'Array'
matrix: 'Array' // choose 'Matrix' (default) or 'Array'
}
const math = create(all, config)
@@ -69,7 +69,8 @@ class Reactor_CSTR {
const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step);
this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
// 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;
}
}
@@ -94,15 +95,15 @@ class Reactor_PFR {
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.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*60); // time step [d]
this.speedUpFactor = 1;
this.timeStep = 1/(24*60*15); // time step [d]
this.speedUpFactor = 60;
this.D_op = this.makeDoperator();
this.D_op = this.makeDoperator(true, true);
this.D2_op = this.makeD2operator();
}
@@ -110,6 +111,10 @@ class Reactor_PFR {
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]
@@ -131,7 +136,6 @@ class Reactor_PFR {
// 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));
@@ -147,10 +151,19 @@ class Reactor_PFR {
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 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));
reaction[0] = Array(13).fill(0.0);
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; });
@@ -158,37 +171,78 @@ class Reactor_PFR {
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);
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!");
}
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
// 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]
this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
return this.state;
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() { // 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);
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 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);
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;
}
}
@@ -199,10 +253,11 @@ class Reactor_PFR {
// 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 < 500) {
// while (N < 5000) {
// console.log(Reactor.tick_fe(0.001));
// N += 1;
// }
module.exports = {Reactor_CSTR, Reactor_PFR};
module.exports = { Reactor_CSTR, Reactor_PFR };