RnD/reactor
3
0
Fork 0
Code Issues 1 Pull Requests Actions Packages Projects Releases Wiki Activity

Rewrite for improved boundary condition

Browse Source
This commit is contained in:
p.vanderwilt
2025-10-01 16:50:48 +02:00
parent d9511dc3c7
commit 3aea0e55c4
1 changed files with 46 additions and 55 deletions
Download Patch File Download Diff File Expand all files Collapse all files

101
src/specificClass.js
View File

@@ -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;
}
}