From 62b034fb76421acd37e0b89474bb28216b2f8455 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Thu, 19 Jun 2025 20:55:42 +0200 Subject: [PATCH 1/5] Added speed-up factor --- dependencies/reactor_class.js | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index db8b3bb..1b2e3d9 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -16,7 +16,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 = 30; } set setInfluent(input) { // setter for C_in (WIP) @@ -39,20 +40,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; } } From 70531a3a597b6d68514b62f79621c5bc2ebb0933 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 23 Jun 2025 16:58:02 +0200 Subject: [PATCH 2/5] Add support for multiple reactor types (CSTR and PFR) with corresponding properties (Dichelet BC for now) --- advanced-reactor.html | 43 ++++++++++++- advanced-reactor.js | 74 ++++++++++++++++------- dependencies/reactor_class.js | 110 +++++++++++++++++++++++++++++++++- 3 files changed, 204 insertions(+), 23 deletions(-) diff --git a/advanced-reactor.html b/advanced-reactor.html index 73d66c2..5b4d41b 100644 --- a/advanced-reactor.html +++ b/advanced-reactor.html @@ -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 }, @@ -45,6 +48,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 +94,22 @@

Reactor properties

+
+ + +
+
+ + +
+
+ + +
diff --git a/advanced-reactor.js b/advanced-reactor.js index edfad14..68e9674 100644 --- a/advanced-reactor.js +++ b/advanced-reactor.js @@ -7,26 +7,60 @@ module.exports = function(RED) { const Reactor = 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( + 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( + parseFloat(config.volume), + parseFloat(config.L), + 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; diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 1b2e3d9..aa7f9dd 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -5,7 +5,6 @@ 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] @@ -17,7 +16,7 @@ class Reactor_CSTR { this.currentTime = Date.now(); // milliseconds since epoch [ms] this.timeStep = 1/(24*60*15); // time step [d] - this.speedUpFactor = 30; + this.speedUpFactor = 1; } set setInfluent(input) { // setter for C_in (WIP) @@ -69,6 +68,113 @@ 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()) + + 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 = 1; + + this.D_op = makeDoperator(); + this.D2_op = 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}, 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 + if (math.sum(this.Fs) > 0) { + this.state[0] = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0] // Dichelet boundary condition + } + + 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(this.asm.compute_dC); + 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]); }); + } + + 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.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]); + 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.diag(Array(this.n_x).fill(2), 0); + const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]); + const B = math.diag(Array(this.n_x).fill(-1), -1).resize([this.n_x, this.n_x]); + I[0, 0] = 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]; From e6c1e21c166ac90751e8e320344bf741b9b4bfc8 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 23 Jun 2025 17:46:55 +0200 Subject: [PATCH 3/5] Implement Danckwerts boundary condition in tick_fe method for Reactor_PFR --- dependencies/reactor_class.js | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index aa7f9dd..bb5a96c 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -137,14 +137,11 @@ class Reactor_PFR { } tick_fe(time_step) { // tick reactor state using forward Euler method - if (math.sum(this.Fs) > 0) { - this.state[0] = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0] // Dichelet boundary condition - } - 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(this.asm.compute_dC); - const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0)) + 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; }); @@ -152,6 +149,15 @@ 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); + 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 @@ -161,6 +167,8 @@ class Reactor_PFR { makeDoperator() { // create the upwind scheme gradient operator const I = math.identity(this.n_x); const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]); + I[0, 0] = 0; + I[0, 1] = 0; I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L return math.add(I, A); } @@ -169,7 +177,8 @@ class Reactor_PFR { const I = math.diag(Array(this.n_x).fill(2), 0); const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]); const B = math.diag(Array(this.n_x).fill(-1), -1).resize([this.n_x, this.n_x]); - I[0, 0] = 1; + I[0, 0] = 0; + I[0, 1] = 0; return math.add(I, A, B); } } From e5c9010093e2cf8df67cead7d4ded041b1e758d1 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 24 Jun 2025 11:20:28 +0200 Subject: [PATCH 4/5] Fixed various bugs --- dependencies/asm3_class.js | 118 +++++++++++++++++----------------- dependencies/reactor_class.js | 41 +++++++----- 2 files changed, 82 insertions(+), 77 deletions(-) diff --git a/dependencies/asm3_class.js b/dependencies/asm3_class.js index 9ad03d9..9b7762a 100644 --- a/dependencies/asm3_class.js +++ b/dependencies/asm3_class.js @@ -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 diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index bb5a96c..af3289e 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -1,5 +1,11 @@ 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 { @@ -93,8 +99,8 @@ class Reactor_PFR { this.timeStep = 1/(24*60*15); // time step [d] this.speedUpFactor = 1; - this.D_op = makeDoperator(); - this.D2_op = makeD2operator(); + this.D_op = this.makeDoperator(); + this.D2_op = this.makeD2operator(); } set setInfluent(input) { // setter for C_in (WIP) @@ -139,7 +145,7 @@ 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 reaction = this.state.map(this.asm.compute_dC); + const reaction = this.state.map((row) => this.asm.compute_dC(row)); reaction[0] = Array(13).fill(0.0); const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0)); @@ -154,7 +160,7 @@ class Reactor_PFR { 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); + 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); } @@ -166,19 +172,19 @@ class Reactor_PFR { makeDoperator() { // create the upwind scheme gradient operator const I = math.identity(this.n_x); - const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]); - I[0, 0] = 0; - I[0, 1] = 0; - I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L + 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.diag(Array(this.n_x).fill(2), 0); - const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]); - const B = math.diag(Array(this.n_x).fill(-1), -1).resize([this.n_x, this.n_x]); - I[0, 0] = 0; - I[0, 1] = 0; + 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); } } @@ -187,9 +193,10 @@ class Reactor_PFR { // 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; From 6b57a46aabd54511c7f3db31b9d480062bb18272 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 24 Jun 2025 12:32:11 +0200 Subject: [PATCH 5/5] Add typed input fields for reactor length and resolution in advanced-reactor, fixed NaN bug in reactor length --- advanced-reactor.html | 8 ++++++++ advanced-reactor.js | 8 ++++---- dependencies/reactor_class.js | 15 +++++++++------ 3 files changed, 21 insertions(+), 10 deletions(-) diff --git a/advanced-reactor.html b/advanced-reactor.html index 5b4d41b..512f3dd 100644 --- a/advanced-reactor.html +++ b/advanced-reactor.html @@ -40,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"] diff --git a/advanced-reactor.js b/advanced-reactor.js index 68e9674..8348137 100644 --- a/advanced-reactor.js +++ b/advanced-reactor.js @@ -5,13 +5,13 @@ module.exports = function(RED) { let name = config.name; - const Reactor = require('./dependencies/reactor_class'); + const { Reactor_CSTR, Reactor_PFR } = require('./dependencies/reactor_class'); let new_reactor; switch (config.reactor_type) { case "CSTR": - new_reactor = new Reactor( + new_reactor = new Reactor_CSTR( parseFloat(config.volume), parseInt(config.n_inlets), parseFloat(config.kla), @@ -33,9 +33,9 @@ module.exports = function(RED) { ); break; case "PFR": - new_reactor = new Reactor( + new_reactor = new Reactor_PFR( parseFloat(config.volume), - parseFloat(config.L), + parseFloat(config.length), parseInt(config.resolution_L), parseInt(config.n_inlets), parseFloat(config.kla), diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index af3289e..4027690 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -87,16 +87,19 @@ class Reactor_PFR { 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.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*15); // time step [d] + this.timeStep = 1/(24*60*60); // time step [d] this.speedUpFactor = 1; this.D_op = this.makeDoperator(); @@ -118,7 +121,7 @@ class Reactor_PFR { } 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.at(-1)}, timestamp: this.currentTime}; } calcOTR(S_O, T=20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C @@ -145,10 +148,10 @@ 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 reaction = this.state.map((row) => this.asm.compute_dC(row)); + 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 { @@ -202,4 +205,4 @@ class Reactor_PFR { // N += 1; // } -module.exports = Reactor_CSTR; \ No newline at end of file +module.exports = {Reactor_CSTR, Reactor_PFR}; \ No newline at end of file