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

Merge branch 'boundary-conditions'

Browse Source
This commit is contained in:
HorriblePerson555
2025-10-16 15:36:37 +02:00
parent 442ddc60ed f44bac9aab
commit 670c4deacb
1 changed files with 97 additions and 78 deletions
Download Patch File Download Diff File Expand all files Collapse all files

175
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 {
@@ -27,6 +28,10 @@ class Reactor {
this.measurements = new MeasurementContainer();
this.upstreamReactor = null;
this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
this.parent = []; // Gets assigned via child registration
this.upstreamReactor = null;
this.downstreamReactor = null;
this.asm = new ASM3();
@@ -41,7 +46,7 @@ class Reactor {
this.currentTime = Date.now(); // milliseconds since epoch [ms]
this.timeStep = 1 / (24*60*60) * this.config.timeStep; // time step in seconds, converted to days.
this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second
this.speedUpFactor = 100; // speed up factor for simulation, 60 means 1 minute per simulated second
}
/**
@@ -113,24 +118,18 @@ class Reactor {
}
}
_connectMeasurement(measurement) {
if (!measurement) {
_connectMeasurement(measurementChild) {
if (!measurementChild) {
this.logger.warn("Invalid measurement provided.");
return;
}
let position;
if (measurement.config.functionality.distance !== 'undefined') {
position = measurement.config.functionality.distance;
} else {
position = measurement.config.functionality.positionVsParent;
}
const measurementType = measurement.config.asset.type;
const key = `${measurementType}_${position}`;
const position = measurementChild.config.functionality.positionVsParent;
const measurementType = measurementChild.config.asset.type;
const eventName = `${measurementType}.measured.${position}`;
// Register event listener for measurement updates
measurement.measurements.emitter.on(eventName, (eventData) => {
measurementChild.measurements.emitter.on(eventName, (eventData) => {
this.logger.debug(`${position} ${measurementType} from ${eventData.childName}: ${eventData.value} ${eventData.unit}`);
// Store directly in parent's measurement container
@@ -145,15 +144,24 @@ class Reactor {
}
_connectReactor(reactor) {
if (!reactor) {
_connectReactor(reactorChild) {
if (!reactorChild) {
this.logger.warn("Invalid reactor provided.");
return;
}
this.upstreamReactor = reactor;
if (reactorChild.functionality.positionVsParent != "upstream") {
this.logger.warn("Reactor children of reactors should always be upstream.");
}
reactor.emitter.on("stateChange", (data) => {
if (math.abs(reactorChild.d_x - this.d_x) / this.d_x < 0.025) {
this.logger.warn("Significant grid sizing discrepancies between adjacent reactors! Change resolutions to match reactors grid step, or implement boundary value interpolation.");
}
this.upstreamReactor = reactorChild;
reactorChild.downstreamReactor = this;
reactorChild.emitter.on("stateChange", (data) => {
this.logger.debug(`State change of upstream reactor detected.`);
this.updateState(data);
});
@@ -246,11 +254,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();
@@ -262,15 +274,16 @@ class Reactor_PFR extends Reactor {
* @param {object} input - Input object (msg) containing payload with dispersion value [m2 d-1].
*/
set setDispersion(input) {
this.D = input.payload;
this.D = this._constrainDispersion(input.payload);
}
updateState(newTime) {
super.updateState(newTime);
let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
// let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
this.D = this._constrainDispersion(this.D);
let Co_D = this.D*this.timeStep/(this.d_x*this.d_x);
(Pe_local >= 2) && this.logger.warn(`Local Péclet number (${Pe_local}) is too high! Increase reactor resolution.`);
// (Pe_local >= 2) && this.logger.warn(`Local Péclet number (${Pe_local}) is too high! Increase reactor resolution.`);
(Co_D >= 0.5) && this.logger.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`);
if(DEBUG) {
@@ -288,25 +301,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.temperature));
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++) {
transfer[i][S_O_INDEX] = this._calcOTR(this.state[i][S_O_INDEX], this.temperature) * this.n_x/(this.n_x-2);
for (let i = BC_PADDING+1; i < BC_PADDING+this.n_x - 1; i++) {
transfer[i][S_O_INDEX] = this._calcOTR(this.extendedState[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");
@@ -317,13 +331,14 @@ class Reactor_PFR extends Reactor {
}
this.state = this._arrayClip2Zero(stateNew);
this.state.forEach((row, i) => this.extendedState[i+BC_PADDING] = row);
return stateNew;
}
_updateMeasurement(measurementType, value, position, context) {
switch(measurementType) {
case "quantity (oxygen)":
let grid_pos = Math.round(position / this.config.length * this.n_x);
let grid_pos = Math.round(context.distance / this.config.length * this.n_x);
this.state[grid_pos][S_O_INDEX] = value; // naive approach for reconciling measurements and simulation
break;
default:
@@ -335,57 +350,51 @@ 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) {
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])));
_applyBoundaryConditions() {
if (this.upstreamReactor) {
for (let i = 0; i < BC_PADDING; i++) {
this.extendedState[i] = this.upstreamReactor.state.at(i-BC_PADDING);
}
} else {
state[0] = state[1];
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);
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])));
this.extendedState[BC_PADDING-1] = math.add(math.multiply(2, this.extendedState[BC_PADDING]), math.multiply(-2, this.extendedState[BC_PADDING+2]), this.extendedState[BC_PADDING+3]);
} else {
for (let i = 0; i < BC_PADDING; i++) {
this.extendedState[i] = this.extendedState[BC_PADDING];
}
}
}
if (this.downstreamReactor) {
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)
for (let i = 0; i < BC_PADDING; i++) {
this.extendedState[BC_PADDING+this.n_x+i] = this.extendedState.at(-1-BC_PADDING);
}
}
// Neumann BC (no flux)
state[this.n_x-1] = state[this.n_x-2];
}
/**
* 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;
}
/**
@@ -393,14 +402,24 @@ 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;
}
_constrainDispersion(D) {
const Dmin = math.sum(this.Fs) * this.d_x / (1.999 * this.A);
if (D < Dmin) {
this.logger.warn(`Local Péclet number too high! Constraining given dispersion (${D}) to minimal value (${Dmin}).`);
return Dmin;
}
return D;
}
}
module.exports = { Reactor_CSTR, Reactor_PFR };