From 3aea0e55c47183acec02d18be7ba58f19570a206 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Wed, 1 Oct 2025 16:50:48 +0200 Subject: [PATCH] Rewrite for improved boundary condition --- src/specificClass.js | 101 ++++++++++++++++++++----------------------- 1 file changed, 46 insertions(+), 55 deletions(-) diff --git a/src/specificClass.js b/src/specificClass.js index ac562d2..9a7e41d 100644 --- a/src/specificClass.js +++ b/src/specificClass.js @@ -12,6 +12,7 @@ const math = create(all, mathConfig); const S_O_INDEX = 0; const NUM_SPECIES = 13; +const BC_PADDING = 2; const DEBUG = false; class Reactor { @@ -250,11 +251,15 @@ class Reactor_PFR extends Reactor { this.alpha = config.alpha; - this.state = Array.from(Array(this.n_x), () => config.initialState.slice()) + this.state = Array.from(Array(this.n_x), () => config.initialState.slice()); + this.extendedState = Array.from(Array(this.n_x + 2*BC_PADDING), () => new Array(NUM_SPECIES).fill(0)); + + // initialise extended state + this.state.forEach((row, i) => this.extendedState[i+BC_PADDING] = row); this.D = 0.0; // axial dispersion [m2 d-1] - this.D_op = this._makeDoperator(true, true); + this.D_op = this._makeDoperator(); assertNoNaN(this.D_op, "Derivative operator"); this.D2_op = this._makeD2operator(); @@ -292,25 +297,26 @@ class Reactor_PFR extends Reactor { * @returns {Array} - New reactor state. */ tick(time_step) { - const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_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, this.temperature)); - const transfer = Array.from(Array(this.n_x), () => new Array(NUM_SPECIES).fill(0)); + this._applyBoundaryConditions(); + + const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.extendedState); + const advection = math.multiply(-1 * math.sum(this.Fs) / (this.A*this.d_x), this.D_op, this.extendedState); + const reaction = this.extendedState.map((state_slice) => this.asm.compute_dC(state_slice, this.temperrature)); + const transfer = Array.from(Array(this.n_x+2*BC_PADDING), () => new Array(NUM_SPECIES).fill(0)); if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR - for (let i = 1; i < this.n_x - 1; i++) { + for (let i = BC_PADDING+1; i < BC_PADDING+this.n_x - 1; i++) { transfer[i][S_O_INDEX] = this.OTR * this.n_x/(this.n_x-2); } } else { - for (let i = 1; i < this.n_x - 1; i++) { + for (let i = BC_PADDING+1; i < BC_PADDING+this.n_x - 1; i++) { transfer[i][S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX], this.temperature) * this.n_x/(this.n_x-2); } } - const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); + const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer).slice(BC_PADDING, this.n_x+BC_PADDING), time_step); const stateNew = math.add(this.state, dC_total); - this._applyBoundaryConditions(stateNew); if (DEBUG) { assertNoNaN(dispersion, "dispersion"); @@ -339,66 +345,50 @@ class Reactor_PFR extends Reactor { * Apply boundary conditions to the reactor state. * for inlet, apply generalised Danckwerts BC, if there is not flow, apply Neumann BC with no flux * for outlet, apply regular Danckwerts BC (Neumann BC with no flux) - * @param {Array} state - Current reactor state without enforced BCs. */ - _applyBoundaryConditions(state) { + _applyBoundaryConditions() { if (this.upstreamReactor) { - state[0] = this.upstreamReactor.state.at(-1); + for (let i = 0; i < BC_PADDING; i++) { + this.extendedState[i] = this.upstreamReactor.state.at(i-BC_PADDING); + } } else { 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_dispersion_term = (1-this.alpha)*this.D*this.A/(math.sum(this.Fs)*this.d_x); - state[0] = math.multiply(1/(1+BC_dispersion_term), math.add(BC_C_in, math.multiply(BC_dispersion_term, state[1]))); + this.extendedState[BC_PADDING] = math.multiply(1/(1+BC_dispersion_term), math.add(BC_C_in, math.multiply(BC_dispersion_term, this.extendedState[BC_PADDING+1]))); } else { - state[0] = state[1]; + this.extendedState[BC_PADDING] = this.extendedState[BC_PADDING+1]; } } if (this.downstreamReactor) { - state[this.n_x-1] = this.downstreamReactor.state[0]; + for (let i = 0; i < BC_PADDING; i++) { + this.extendedState[this.n_x+BC_PADDING+i] = this.downstreamReactor.state[i]; + } } else { // Neumann BC (no flux) - state[this.n_x-1] = state.at(-2); + for (let i = 0; i < BC_PADDING; i++) { + this.extendedState[BC_PADDING+this.n_x+i] = this.extendedState.at(-1-BC_PADDING); + } } + + this.state.forEach((row, i) => this.extendedState[i+BC_PADDING] = row); } /** * Create finite difference first derivative operator. - * @param {boolean} central - Use central difference scheme if true, otherwise use upwind scheme. - * @param {boolean} higher_order - Use higher order scheme if true, otherwise use first order scheme. * @returns {Array} - First derivative operator matrix. */ - _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; - } + _makeDoperator() { // create gradient operator + const D_size = this.n_x+2*BC_PADDING; + const I = math.resize(math.diag(Array(D_size).fill(1/12), -2), [D_size, D_size]); + const A = math.resize(math.diag(Array(D_size).fill(-2/3), -1), [D_size, D_size]); + const B = math.resize(math.diag(Array(D_size).fill(2/3), 1), [D_size, D_size]); + const C = math.resize(math.diag(Array(D_size).fill(-1/12), 2), [D_size, D_size]); + const D = math.add(I, A, B, C); + // set by BCs elsewhere + D.forEach((row, i) => i < BC_PADDING || i >= this.n_x+BC_PADDING ? row.fill(0) : row); + return D; } /** @@ -406,12 +396,13 @@ class Reactor_PFR extends Reactor { * @returns {Array} - Second derivative operator matrix. */ _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 D_size = this.n_x+2*BC_PADDING; + const I = math.diag(Array(D_size).fill(-2), 0); + const A = math.resize(math.diag(Array(D_size).fill(1), 1), [D_size, D_size]); + const B = math.resize(math.diag(Array(D_size).fill(1), -1), [D_size, D_size]); 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); + // set by BCs elsewhere + D2.forEach((row, i) => i < BC_PADDING || i >= this.n_x+BC_PADDING ? row.fill(0) : row); return D2; } }