diff --git a/reactor.html b/reactor.html index e8c57c5..1b122db 100644 --- a/reactor.html +++ b/reactor.html @@ -102,6 +102,19 @@ } else { $(".PFR").show(); } + + const updateDx = () => { + const length = parseFloat($("#node-input-length").val()) || 0; + const resolution = parseFloat($("#node-input-resolution_L").val()) || 1; + const dx = resolution > 0 ? (length / resolution).toFixed(6) : "N/A"; + $("#dx-output").text(dx + " m"); + }; + + // Set up event listeners for real-time updates + $("#node-input-length, #node-input-resolution_L").on("change keyup", updateDx); + + // Initial calculation + updateDx(); }, oneditsave: function() { // save logger fields @@ -144,6 +157,10 @@ +
+ + -- +

Internal mass transfer calculation (optional)

diff --git a/src/nodeClass.js b/src/nodeClass.js index cc42911..f98b5fe 100644 --- a/src/nodeClass.js +++ b/src/nodeClass.js @@ -34,7 +34,6 @@ class nodeClass { switch (msg.topic) { case "clock": this.source.updateState(msg.timestamp); - send([msg, null, null]); break; case "Fluent": this.source.setInfluent = msg; @@ -42,9 +41,6 @@ class nodeClass { case "OTR": this.source.setOTR = msg; break; - case "Temperature": - this.source.setTemperature = msg; - break; case "Dispersion": this.source.setDispersion = msg; break; diff --git a/src/reaction_modules/asm3_class Koch.js b/src/reaction_modules/asm3_class Koch.js index 741c3e7..ff6a8d2 100644 --- a/src/reaction_modules/asm3_class Koch.js +++ b/src/reaction_modules/asm3_class Koch.js @@ -7,6 +7,62 @@ const ASM_CONSTANTS = { NUM_SPECIES: 13 }; +const KINETIC_CONSTANTS = { + // Hydrolysis + k_H: 9., // 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: 12., // storage rate constant [g S_S g-1 X_H d-1] + nu_NO: 0.5, // 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: 10., // saturation constant S_s [g COD m-3] + K_STO: 0.1, // saturation constant X_STO [g X_STO g-1 X_H] + mu_H_max: 3., // 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.3, // aerobic respiration rate [d-1] + b_H_NO: 0.15, // anoxic respiration rate [d-1] + b_STO_O: 0.3, // aerobic respitation rate X_STO [d-1] + b_STO_NO: 0.15, // anoxic respitation rate X_STO [d-1] + // Autotrophs + mu_A_max: 1.3, // maximum specific growth rate [d-1] + K_A_NH: 1.4, // 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.20, // aerobic respiration rate [d-1] + b_A_NO: 0.10 // anoxic respiration rate [d-1] +}; + +const STOICHIOMETRIC_CONSTANTS = { + // Fractions + 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.80, // aerobic yield X_STO per S_S [g X_STO g-1 S_S] + Y_STO_NO: 0.70, // anoxic yield X_STO per S_S [g X_STO g-1 S_S] + Y_H_O: 0.80, // aerobic yield X_H per X_STO [g X_H g-1 X_STO] + Y_H_NO: 0.65, // 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.04, // nitrogen content X_I [g N g-1 X_I] + i_NXS: 0.03, // 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] +}; + /** * ASM3 class for the Activated Sludge Model No. 3 (ASM3). Using Koch et al. 2000 parameters. */ @@ -17,65 +73,13 @@ class ASM3 { * Kinetic parameters for ASM3 at 20 C. Using Koch et al. 2000 parameters. * @property {Object} kin_params - Kinetic parameters */ - this.kin_params = { - // Hydrolysis - k_H: 9., // 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: 12., // storage rate constant [g S_S g-1 X_H d-1] - nu_NO: 0.5, // 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: 10., // saturation constant S_s [g COD m-3] - K_STO: 0.1, // saturation constant X_STO [g X_STO g-1 X_H] - mu_H_max: 3., // 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.3, // aerobic respiration rate [d-1] - b_H_NO: 0.15, // anoxic respiration rate [d-1] - b_STO_O: 0.3, // aerobic respitation rate X_STO [d-1] - b_STO_NO: 0.15, // anoxic respitation rate X_STO [d-1] - // Autotrophs - mu_A_max: 1.3, // maximum specific growth rate [d-1] - K_A_NH: 1.4, // 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.20, // aerobic respiration rate [d-1] - b_A_NO: 0.10 // anoxic respiration rate [d-1] - }; + this.kin_params = KINETIC_CONSTANTS; /** * Stoichiometric and composition parameters for ASM3. Using Koch et al. 2000 parameters. * @property {Object} stoi_params - Stoichiometric parameters */ - this.stoi_params = { - // Fractions - 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.80, // aerobic yield X_STO per S_S [g X_STO g-1 S_S] - Y_STO_NO: 0.70, // anoxic yield X_STO per S_S [g X_STO g-1 S_S] - Y_H_O: 0.80, // aerobic yield X_H per X_STO [g X_H g-1 X_STO] - Y_H_NO: 0.65, // 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.04, // nitrogen content X_I [g N g-1 X_I] - i_NXS: 0.03, // 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_params = STOICHIOMETRIC_CONSTANTS; /** * Temperature theta parameters for ASM3. Using Koch et al. 2000 parameters. @@ -215,4 +219,4 @@ class ASM3 { } } -module.exports = { ASM3, ASM_CONSTANTS }; \ No newline at end of file +module.exports = { ASM3, ASM_CONSTANTS, KINETIC_CONSTANTS, STOICHIOMETRIC_CONSTANTS }; \ No newline at end of file diff --git a/src/reaction_modules/asm3_class.js b/src/reaction_modules/asm3_class.js index db39890..53756bf 100644 --- a/src/reaction_modules/asm3_class.js +++ b/src/reaction_modules/asm3_class.js @@ -7,6 +7,62 @@ const ASM_CONSTANTS = { NUM_SPECIES: 13 }; +const KINETIC_CONSTANTS = { + // 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] +}; + +const STOICHIOMETRIC_CONSTANTS = { + // Fractions + 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] +}; + /** * ASM3 class for the Activated Sludge Model No. 3 (ASM3). */ @@ -17,65 +73,13 @@ class ASM3 { * Kinetic parameters for ASM3 at 20 C. * @property {Object} kin_params - Kinetic parameters */ - this.kin_params = { - // 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.kin_params = KINETIC_CONSTANTS; /** * Stoichiometric and composition parameters for ASM3. * @property {Object} stoi_params - Stoichiometric parameters */ - this.stoi_params = { - // Fractions - 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_params = STOICHIOMETRIC_CONSTANTS; /** * Temperature theta parameters for ASM3. @@ -215,4 +219,4 @@ class ASM3 { } } -module.exports = { ASM3, ASM_CONSTANTS }; \ No newline at end of file +module.exports = { ASM3, ASM_CONSTANTS, KINETIC_CONSTANTS, STOICHIOMETRIC_CONSTANTS }; \ No newline at end of file diff --git a/src/specificClass.js b/src/specificClass.js index 23c1289..5f6e2c7 100644 --- a/src/specificClass.js +++ b/src/specificClass.js @@ -10,9 +10,9 @@ const mathConfig = { const math = create(all, mathConfig); -const BC_PADDING = 2; +const BC_PADDING = 2; // Boundary Condition padding for open boundaries in extendedState variable const DEBUG = false; -const DAY2MS = 1000 * 60 * 60 * 24; +const DAY2MS = 1000 * 60 * 60 * 24; // convert between days and milliseconds class Reactor { /** @@ -25,13 +25,14 @@ class Reactor { this.logger = new logger(this.config.general.logging.enabled, this.config.general.logging.logLevel, config.general.name); this.emitter = new EventEmitter(); this.measurements = new MeasurementContainer(); - this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility + this.childRegistrationUtils = new childRegistrationUtils(this); // child registration utility + // placeholder variables for children and parents this.upstreamReactor = null; this.downstreamReactor = null; this.returnPump = null; - this.asm = new ASM3(); + this.asm = new ASM3(); // Reaction model this.volume = config.volume; // fluid volume reactor [m3] @@ -42,9 +43,9 @@ class Reactor { this.kla = config.kla; // if NaN, use externaly provided OTR [d-1] - this.currentTime = Date.now(); // milliseconds since epoch [ms] + this.currentTime = null; // milliseconds since epoch [ms] this.timeStep = 1 / (24*60*60) * this.config.timeStep; // time step in seconds, converted to days. - this.speedUpFactor = 100; // speed up factor for simulation, 60 means 1 minute per simulated second + this.speedUpFactor = 1; // speed up factor for simulation, 60 means 1 minute per simulated second } /** @@ -113,6 +114,11 @@ class Reactor { } } + /** + * Register child function required for child registration. + * @param {object} child + * @param {string} softwareType + */ registerChild(child, softwareType) { if(!child) { this.logger.error(`Invalid ${softwareType} child provided.`); @@ -161,18 +167,14 @@ class Reactor { _connectReactor(reactorChild) { if (reactorChild.config.functionality.positionVsParent != "upstream") { - this.logger.warn("Reactor children of reactors should always be upstream."); - } - - 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.logger.warn("Reactor children of other reactors should always be upstream!"); } // set upstream and downstream reactor variable in current and child nodes respectively for easy access this.upstreamReactor = reactorChild; reactorChild.downstreamReactor = this; - reactorChild.emitter.on("stateChange", (eventData) => { + reactorChild.emitter.on("stateChange", (eventData) => { // Triggers state update in downstream reactor. this.logger.debug(`State change of upstream reactor detected.`); this.updateState(eventData); }); @@ -203,20 +205,32 @@ class Reactor { * Update the reactor state based on the new time. * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch. */ - updateState(newTime) { // expect update with timestamp - if (this.upstreamReactor) { - this.setInfluent = this.upstreamReactor.getEffluent[0]; // grab main effluent upstream reactor + updateState(newTime) { + if (!this.currentTime) { // initialise currentTime variable + this.currentTime = newTime; + return; + } + + if (this.upstreamReactor) { // grab main effluent upstream reactor + this.setInfluent = this.upstreamReactor.getEffluent[0]; } const n_iter = Math.floor(this.speedUpFactor * (newTime-this.currentTime) / (this.timeStep*DAY2MS)); - if (n_iter) { - let n = 0; - while (n < n_iter) { - this.tick(this.timeStep); - n += 1; - } - this.currentTime += n_iter * this.timeStep * DAY2MS / this.speedUpFactor; - this.emitter.emit("stateChange", this.currentTime); + + if (n_iter == 0) { // no update required, change in currentTime smaller than time step + return; + } + + let n = 0; + while (n < n_iter) { + this.tick(this.timeStep); + n += 1; + } + this.currentTime += n_iter * this.timeStep * DAY2MS / this.speedUpFactor; + this.emitter.emit("stateChange", this.currentTime); // update downstream reactors + + if (this.returnPump) { // update recirculation pump state + this.returnPump.updateSourceSink(); } } } @@ -231,6 +245,23 @@ class Reactor_CSTR extends Reactor { this.state = config.initialState; } + _updateMeasurement(measurementType, value, position, context) { + + switch(measurementType) { + case "quantity (oxygen)": + this.state[ASM_CONSTANTS.S_O_INDEX] = value; + break; + case "quantity (ammonium)": + this.state[ASM_CONSTANTS.S_NH_INDEX] = value; + break; + case "quantity (nox)": + this.state[ASM_CONSTANTS.S_NO_INDEX] = value; + break; + default: + super._updateMeasurement(measurementType, value, position, context); + } + } + /** * Tick the reactor state using the forward Euler method. * @param {number} time_step - Time step for the simulation [d]. @@ -241,7 +272,7 @@ class Reactor_CSTR extends Reactor { const outflow = math.multiply(-1 * math.sum(this.Fs) / this.volume, this.state); const reaction = this.asm.compute_dC(this.state, this.temperature); const transfer = Array(ASM_CONSTANTS.NUM_SPECIES).fill(0.0); - transfer[ASM_CONSTANTS.S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[S_O_INDEX], this.temperature); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR + transfer[ASM_CONSTANTS.S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[ASM_CONSTANTS.S_O_INDEX], this.temperature); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step) this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations @@ -290,13 +321,23 @@ class Reactor_PFR extends Reactor { this.D = this._constrainDispersion(input.payload); } + _connectReactor(reactorChild) { + 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."); + } + super._connectReactor(reactorChild); + } + + /** + * Update the reactor state based on the new time. Performs checks specific to PFR. + * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch. + */ updateState(newTime) { super.updateState(newTime); - // let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A) - this.D = this._constrainDispersion(this.D); + + this.D = this._constrainDispersion(this.D); // constrains D to minimum dispersion, so that local Péclet number is always above 2 const 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.`); (Co_D >= 0.5) && this.logger.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`); if(DEBUG) { @@ -375,8 +416,8 @@ class Reactor_PFR extends Reactor { */ _applyBoundaryConditions() { // Upstream BC - if (this.upstreamReactor) { - // Open boundary + if (this.upstreamReactor && this.upstreamReactor.config.reactor_type == "PFR") { + // Open boundary, if upstream reactor is PFR this.extendedState.splice(0, BC_PADDING, ...this.upstreamReactor.state.slice(-BC_PADDING)); } else { if (math.sum(this.Fs) > 0) { @@ -384,7 +425,7 @@ class Reactor_PFR extends Reactor { const BC_C_in = math.multiply(1 / math.sum(this.Fs), [this.Fs], this.Cs_in)[0]; const BC_dispersion_term = 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]))); - // Numerical boundary condition + // Numerical boundary condition (first-order accurate) 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 { // Neumann BC (no flux) @@ -393,8 +434,8 @@ class Reactor_PFR extends Reactor { } // Downstream BC - if (this.downstreamReactor) { - // Open boundary + if (this.downstreamReactor && this.downstreamReactor.config.reactor_type == "PFR") { + // Open boundary, if downstream reactor is PFR this.extendedState.splice(this.n_x+BC_PADDING, BC_PADDING, ...this.downstreamReactor.state.slice(0, BC_PADDING)); } else { // Neumann BC (no flux) @@ -404,7 +445,6 @@ class Reactor_PFR extends Reactor { /** * Create finite difference first derivative operator. - * @returns {Array} - First derivative operator matrix. */ _makeDoperator() { // create gradient operator const D_size = this.n_x+2*BC_PADDING; @@ -420,7 +460,6 @@ class Reactor_PFR extends Reactor { /** * Create central finite difference second derivative operator. - * @returns {Array} - Second derivative operator matrix. */ _makeD2operator() { // create the central second derivative operator const D_size = this.n_x+2*BC_PADDING; @@ -433,6 +472,9 @@ class Reactor_PFR extends Reactor { return D2; } + /** + * Constrains dispersion so that local Péclet number stays below 2. Needed for stable central differencing method. + */ _constrainDispersion(D) { const Dmin = math.sum(this.Fs) * this.d_x / (1.999 * this.A); if (D < Dmin) {