From e5c9010093e2cf8df67cead7d4ded041b1e758d1 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 24 Jun 2025 11:20:28 +0200 Subject: [PATCH] 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;