diff --git a/advanced-reactor.js b/advanced-reactor.js index 8348137..e2eca51 100644 --- a/advanced-reactor.js +++ b/advanced-reactor.js @@ -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); } diff --git a/dependencies/asm3_class.js b/dependencies/asm3_class.js index 9b7762a..f5ac1df 100644 --- a/dependencies/asm3_class.js +++ b/dependencies/asm3_class.js @@ -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; diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index 4027690..76a6340 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -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}; \ No newline at end of file +module.exports = { Reactor_CSTR, Reactor_PFR }; \ No newline at end of file