From 01380c309f757622c5ad194f42ee938cedde2bc7 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 7 Jul 2025 14:47:50 +0200 Subject: [PATCH 1/4] Add boundary condition input and update reactor configuration handling --- advanced-reactor.html | 17 +++++++++++++++++ src/nodeClass.js | 1 + src/specificClass.js | 27 ++++++++++++++++++++++----- 3 files changed, 40 insertions(+), 5 deletions(-) diff --git a/advanced-reactor.html b/advanced-reactor.html index 512f3dd..6594947 100644 --- a/advanced-reactor.html +++ b/advanced-reactor.html @@ -8,6 +8,7 @@ volume: { value: 0., required: true }, length: { value: 0.}, resolution_L: { value: 0.}, + boundary_condition: {value: "Generalised"}, n_inlets: { value: 1, required: true}, kla: { value: null }, S_O_init: { value: 0., required: true }, @@ -75,6 +76,18 @@ $(".PFR").show(); } }); + $("#node-input-boundary_condition").typedInput({ + types: [ + { + value: "Generalised", + options: [ + { value: "Generalised", label: "Generalised"}, + { value: "Diriclet", label: "Diriclet"}, + { value: "Danckwerts", label: "Danckwerts"} + ] + } + ] + }) // Set initial visibility on dialog open const initialType = $("#node-input-reactor_type").typedInput("value"); if (initialType === "CSTR") { @@ -117,6 +130,10 @@
+
+
+ +
diff --git a/src/nodeClass.js b/src/nodeClass.js index 7919df2..e3c3fd6 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -70,6 +70,7 @@ class nodeClass { volume: parseFloat(uiConfig.volume), length: parseFloat(uiConfig.length), resolution_L: parseInt(uiConfig.resolution_L), + boundary_condition: uiConfig.boundary_condition, n_inlets: parseInt(uiConfig.n_inlets), kla: parseFloat(uiConfig.kla), initialState: [ diff --git a/src/specificClass.js b/src/specificClass.js index 384aaf0..e670eb5 100644 --- a/src/specificClass.js +++ b/src/specificClass.js @@ -10,7 +10,7 @@ const math = create(all, config); const S_O_INDEX = 0; const NUM_SPECIES = 13; -const DEBUG = false; +const DEBUG = true; class Reactor { /** @@ -30,7 +30,7 @@ class Reactor { this.currentTime = Date.now(); // milliseconds since epoch [ms] this.timeStep = 1 / (24*60*15); // time step [d] - this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second + this.speedUpFactor = 1; // speed up factor for simulation, 60 means 1 minute per simulated second } /** @@ -149,6 +149,8 @@ class Reactor_PFR extends Reactor { this.d_x = this.length / this.n_x; this.A = this.volume / this.length; // crosssectional area [m2] + this.BC = config.boundary_condition; + this.state = Array.from(Array(this.n_x), () => config.initialState.slice()) // console.log("Initial State: ") @@ -178,6 +180,7 @@ class Reactor_PFR extends Reactor { set setInfluent(input) { super.setInfluent = input; if(DEBUG) { + console.log("Inlet state max " + math.max(this.state[0])) 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)); @@ -205,9 +208,23 @@ class Reactor_PFR extends Reactor { const BC_gradient = Array(this.n_x).fill(0); BC_gradient[0] = -1; BC_gradient[1] = 1; - let Pe = this.length * math.sum(this.Fs) / (this.D * this.A); - let residence_time = this.volume/math.sum(this.Fs); - const BC_dispersion = math.multiply((1 - (1 + 4*residence_time/Pe)^0.5) / (Pe*this.d_x), [BC_gradient], state)[0]; + let BC_term = 0; + switch(this.BC) { + case "Diriclet": + BC_term = 0; + break; + case "Danckwerts": + BC_term = this.D*this.A / math.sum(this.Fs); + break; + case "Generalised": + let Pe = this.length * math.sum(this.Fs) / (this.D * this.A); + let residence_time = this.volume/math.sum(this.Fs); + BC_term = (1 - (1 + 4*residence_time/Pe)^0.5) / Pe + break; + default: + console.warn("Unknown boundary condition: " + this.BC); + } + const BC_dispersion = math.multiply(BC_term/this.d_x, [BC_gradient], state)[0]; state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val); } else { // Neumann BC (no flux) state[0] = state[1]; From 01318a2d3b2d968f2efed840e6dfe2dd08d73971 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Mon, 7 Jul 2025 14:58:52 +0200 Subject: [PATCH 2/4] Fix spelling of "Dirichlet" in advanced-reactor.html and specificClass.js --- advanced-reactor.html | 2 +- src/specificClass.js | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/advanced-reactor.html b/advanced-reactor.html index 6594947..9ea7e23 100644 --- a/advanced-reactor.html +++ b/advanced-reactor.html @@ -82,7 +82,7 @@ value: "Generalised", options: [ { value: "Generalised", label: "Generalised"}, - { value: "Diriclet", label: "Diriclet"}, + { value: "Dirichlet", label: "Dirichlet"}, { value: "Danckwerts", label: "Danckwerts"} ] } diff --git a/src/specificClass.js b/src/specificClass.js index e670eb5..1124b0d 100644 --- a/src/specificClass.js +++ b/src/specificClass.js @@ -210,7 +210,7 @@ class Reactor_PFR extends Reactor { BC_gradient[1] = 1; let BC_term = 0; switch(this.BC) { - case "Diriclet": + case "Dirichlet": BC_term = 0; break; case "Danckwerts": From 5c03dddb7923066f705edbd490afd49388c511bf Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 8 Jul 2025 10:03:03 +0200 Subject: [PATCH 3/4] Refactor boundary condition handling to use adjustable parameter alpha in advanced-reactor and specificClass --- advanced-reactor.html | 25 ++++++++++--------------- src/nodeClass.js | 2 +- src/specificClass.js | 20 ++------------------ 3 files changed, 13 insertions(+), 34 deletions(-) diff --git a/advanced-reactor.html b/advanced-reactor.html index 9ea7e23..40ea4f6 100644 --- a/advanced-reactor.html +++ b/advanced-reactor.html @@ -8,7 +8,7 @@ volume: { value: 0., required: true }, length: { value: 0.}, resolution_L: { value: 0.}, - boundary_condition: {value: "Generalised"}, + alpha: {value: 0}, n_inlets: { value: 1, required: true}, kla: { value: null }, S_O_init: { value: 0., required: true }, @@ -76,17 +76,9 @@ $(".PFR").show(); } }); - $("#node-input-boundary_condition").typedInput({ - types: [ - { - value: "Generalised", - options: [ - { value: "Generalised", label: "Generalised"}, - { value: "Dirichlet", label: "Dirichlet"}, - { value: "Danckwerts", label: "Danckwerts"} - ] - } - ] + $("#node-input-alpha").typedInput({ + type:"num", + types:["num"] }) // Set initial visibility on dialog open const initialType = $("#node-input-reactor_type").typedInput("value"); @@ -131,9 +123,12 @@
-
- - +
+

Inlet boundary condition parameter α (α = 0: Danckwerts BC / α = 1: Dirichlet BC)

+
+ + +
diff --git a/src/nodeClass.js b/src/nodeClass.js index e3c3fd6..7ec81d3 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -70,7 +70,7 @@ class nodeClass { volume: parseFloat(uiConfig.volume), length: parseFloat(uiConfig.length), resolution_L: parseInt(uiConfig.resolution_L), - boundary_condition: uiConfig.boundary_condition, + alpha: parseFloat(uiConfig.alpha), n_inlets: parseInt(uiConfig.n_inlets), kla: parseFloat(uiConfig.kla), initialState: [ diff --git a/src/specificClass.js b/src/specificClass.js index 1124b0d..18af512 100644 --- a/src/specificClass.js +++ b/src/specificClass.js @@ -149,7 +149,7 @@ class Reactor_PFR extends Reactor { this.d_x = this.length / this.n_x; this.A = this.volume / this.length; // crosssectional area [m2] - this.BC = config.boundary_condition; + this.alpha = config.alpha; this.state = Array.from(Array(this.n_x), () => config.initialState.slice()) @@ -208,23 +208,7 @@ class Reactor_PFR extends Reactor { const BC_gradient = Array(this.n_x).fill(0); BC_gradient[0] = -1; BC_gradient[1] = 1; - let BC_term = 0; - switch(this.BC) { - case "Dirichlet": - BC_term = 0; - break; - case "Danckwerts": - BC_term = this.D*this.A / math.sum(this.Fs); - break; - case "Generalised": - let Pe = this.length * math.sum(this.Fs) / (this.D * this.A); - let residence_time = this.volume/math.sum(this.Fs); - BC_term = (1 - (1 + 4*residence_time/Pe)^0.5) / Pe - break; - default: - console.warn("Unknown boundary condition: " + this.BC); - } - const BC_dispersion = math.multiply(BC_term/this.d_x, [BC_gradient], state)[0]; + const BC_dispersion = math.multiply((1 - this.alpha) * this.D*this.A / (math.sum(this.Fs) * this.d_x), [BC_gradient], state)[0]; state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val); } else { // Neumann BC (no flux) state[0] = state[1]; From c1e331b5f0134070b51d6a3f2762442d1aa87543 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 8 Jul 2025 11:02:39 +0200 Subject: [PATCH 4/4] Add temperature theta parameters and adjust reaction rate calculations in ASM3 class --- src/reaction_modules/asm3_class.js | 88 +++++++++++++++++++++++------- 1 file changed, 68 insertions(+), 20 deletions(-) diff --git a/src/reaction_modules/asm3_class.js b/src/reaction_modules/asm3_class.js index ab272dc..d228619 100644 --- a/src/reaction_modules/asm3_class.js +++ b/src/reaction_modules/asm3_class.js @@ -69,6 +69,28 @@ class ASM3 { 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] }; + + /** + * Temperature theta parameters for ASM3. + * These parameters are used to adjust reaction rates based on temperature. + * @property {Object} temp_params - Temperature theta parameters + */ + this.temp_params = { + // Hydrolysis + theta_H: this._compute_theta(2, 3, 10, 20), + // Heterotrophs + theta_STO: this._compute_theta(2.5, 5, 10, 20), + theta_mu_H: this._compute_theta(1, 2, 10, 20), + theta_b_H_O: this._compute_theta(0.1, 0.2, 10, 20), + theta_b_H_NO: this._compute_theta(0.05, 0.1, 10, 20), + theta_b_STO_O: this._compute_theta(0.1, 0.2, 10, 20), + theta_b_STO_NO: this._compute_theta(0.05, 0.1, 10, 20), + // Autotrophs + theta_mu_A: this._compute_theta(0.35, 1, 10, 20), + theta_b_A_O: this._compute_theta(0.05, 0.15, 10, 20), + theta_b_A_NO: this._compute_theta(0.02, 0.05, 10, 20) + }; + this.stoi_matrix = this._initialise_stoi_matrix(); } @@ -103,7 +125,7 @@ class ASM3 { * @param {number} K - Half-saturation constant for the reaction species. * @returns {number} - Monod equation rate value for the given concentration and half-saturation constant. */ - _monod(c, K){ + _monod(c, K) { return c / (K + c); } @@ -113,50 +135,76 @@ class ASM3 { * @param {number} K - Half-saturation constant for the reaction species. * @returns {number} - Inverse Monod equation rate value for the given concentration and half-saturation constant. */ - _inv_monod(c, K){ + _inv_monod(c, K) { return K / (K + c); } /** - * Computes the reaction rates for each process reaction based on the current state. + * Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter. + * @param {number} k - Rate constant at 20 degrees Celcius. + * @param {number} theta - Theta parameter. + * @param {number} T - Temperature in Celcius. + * @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation. + */ + _arrhenius(k, theta, T) { + return k * Math.exp(theta*(T-20)); + } + + /** + * Computes the temperature theta parameter based on two rate constants and their corresponding temperatures. + * @param {number} k1 - Rate constant at temperature T1. + * @param {number} k2 - Rate constant at temperature T2. + * @param {number} T1 - Temperature T1 in Celcius. + * @param {number} T2 - Temperature T2 in Celcius. + * @returns {number} - Theta parameter. + */ + _compute_theta(k1, k2, T1, T2) { + return Math.log(k1/k2)/(T1-T2); + } + + /** + * Computes the reaction rates for each process reaction based on the current state and temperature. * @param {Array} state - State vector containing concentrations of reaction species. + * @param {number} [T=20] - Temperature in degrees Celsius (default is 20). * @returns {Array} - Reaction rates for each process reaction. */ - compute_rates(state) { + compute_rates(state, T = 20) { // 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; - + const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params; + // Hydrolysis - rates[0] = X_H == 0 ? 0 : k_H * this._monod(X_S / X_H, K_X) * X_H; + rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * 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; + rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H; + rates[2] = this._arrhenius(k_STO, theta_STO, T) * 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 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H; + rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H; + rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H; + rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * 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; + rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * 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] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A; + rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A; return rates; } /** - * Computes the change in concentrations of reaction species based on the current state. + * Computes the change in concentrations of reaction species based on the current state and temperature. * @param {Array} state - State vector containing concentrations of reaction species. + * @param {number} [T=20] - Temperature in degrees Celsius (default is 20). * @returns {Array} - Change in reaction species concentrations. */ - compute_dC(state) { // compute changes in concentrations + compute_dC(state, T = 20) { // 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)); + return math.multiply(this.stoi_matrix, this.compute_rates(state, T)); } }