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

Resolve merge conflicts from migration

Browse Source
This commit is contained in:
p.vanderwilt
2025-09-22 16:40:22 +02:00
parent ea35038aa1 9147a3f7d0
commit da1cff55ba
16 changed files with 3569 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

136
.gitignore vendored Normal file
View File

@@ -0,0 +1,136 @@
# Logs
logs
*.log
npm-debug.log*
yarn-debug.log*
yarn-error.log*
lerna-debug.log*
.pnpm-debug.log*
# Diagnostic reports (https://nodejs.org/api/report.html)
report.[0-9]*.[0-9]*.[0-9]*.[0-9]*.json
# Runtime data
pids
*.pid
*.seed
*.pid.lock
# Directory for instrumented libs generated by jscoverage/JSCover
lib-cov
# Coverage directory used by tools like istanbul
coverage
*.lcov
# nyc test coverage
.nyc_output
# Grunt intermediate storage (https://gruntjs.com/creating-plugins#storing-task-files)
.grunt
# Bower dependency directory (https://bower.io/)
bower_components
# node-waf configuration
.lock-wscript
# Compiled binary addons (https://nodejs.org/api/addons.html)
build/Release
# Dependency directories
node_modules/
jspm_packages/
# Snowpack dependency directory (https://snowpack.dev/)
web_modules/
# TypeScript cache
*.tsbuildinfo
# Optional npm cache directory
.npm
# Optional eslint cache
.eslintcache
# Optional stylelint cache
.stylelintcache
# Microbundle cache
.rpt2_cache/
.rts2_cache_cjs/
.rts2_cache_es/
.rts2_cache_umd/
# Optional REPL history
.node_repl_history
# Output of 'npm pack'
*.tgz
# Yarn Integrity file
.yarn-integrity
# dotenv environment variable files
.env
.env.development.local
.env.test.local
.env.production.local
.env.local
# parcel-bundler cache (https://parceljs.org/)
.cache
.parcel-cache
# Next.js build output
.next
out
# Nuxt.js build / generate output
.nuxt
dist
# Gatsby files
.cache/
# Comment in the public line in if your project uses Gatsby and not Next.js
# https://nextjs.org/blog/next-9-1#public-directory-support
# public
# vuepress build output
.vuepress/dist
# vuepress v2.x temp and cache directory
.temp
.cache
# vitepress build output
**/.vitepress/dist
# vitepress cache directory
**/.vitepress/cache
# Docusaurus cache and generated files
.docusaurus
# Serverless directories
.serverless/
# FuseBox cache
.fusebox/
# DynamoDB Local files
.dynamodb/
# TernJS port file
.tern-port
# Stores VSCode versions used for testing VSCode extensions
.vscode-test
# yarn v2
.yarn/cache
.yarn/unplugged
.yarn/build-state.yml
.yarn/install-state.gz
.pnp.*

29
README.md
View File

@@ -1,16 +1,17 @@
# reactor
Reactor: Advanced Hydraulic Tank & Biological Process Simulator
A comprehensive reactor class for wastewater treatment simulation featuring plug flow hydraulics, ASM1-ASM3 biological modeling, and multi-sectional concentration tracking. Implements hydraulic retention time calculations, dispersion modeling, and real-time biological reaction kinetics for accurate process simulation.
Key Features:
Plug Flow Hydraulics: Multi-section reactor with configurable sectioning factor and dispersion modeling
ASM1 Integration: Complete biological nutrient removal modeling with 13 state variables (COD, nitrogen, phosphorus)
Dynamic Volume Control: Automatic section management with overflow handling and retention time calculations
Oxygen Transfer: Saturation-limited O2 transfer with Fick's law slowdown effects and solubility curves
Real-time Kinetics: Continuous biological reaction rate calculations with configurable time acceleration
Weighted Averaging: Volume-based concentration mixing for accurate mass balance calculations
Child Registration: Integration with diffuser systems and upstream/downstream reactor networks
Supports complex biological treatment train modeling with temperature compensation, sludge calculations, and comprehensive process monitoring for wastewater treatment plant optimization and regulatory compliance.
Reactor: Advanced Hydraulic Tank & Biological Process Simulator
A comprehensive reactor class for wastewater treatment simulation featuring plug flow hydraulics, ASM1-ASM3 biological modeling, and multi-sectional concentration tracking. Implements hydraulic retention time calculations, dispersion modeling, and real-time biological reaction kinetics for accurate process simulation.
Key Features:
Plug Flow Hydraulics: Multi-section reactor with configurable sectioning factor and dispersion modeling
ASM1 Integration: Complete biological nutrient removal modeling with 13 state variables (COD, nitrogen, phosphorus)
Dynamic Volume Control: Automatic section management with overflow handling and retention time calculations
Oxygen Transfer: Saturation-limited O2 transfer with Fick's law slowdown effects and solubility curves
Real-time Kinetics: Continuous biological reaction rate calculations with configurable time acceleration
Weighted Averaging: Volume-based concentration mixing for accurate mass balance calculations
Child Registration: Integration with diffuser systems and upstream/downstream reactor networks
Supports complex biological treatment train modeling with temperature compensation, sludge calculations, and comprehensive process monitoring for wastewater treatment plant optimization and regulatory compliance.

57
additional_nodes/recirculation-pump.html Normal file
View File

@@ -0,0 +1,57 @@
<script type="text/javascript">
RED.nodes.registerType("recirculation-pump", {
category: "WWTP",
color: "#e4a363",
defaults: {
name: { value: "" },
F2: { value: 0, required: true },
inlet: { value: 1, required: true }
},
inputs: 1,
outputs: 2,
outputLabels: ["Main effluent", "Recirculation effluent"],
icon: "font-awesome/fa-random",
label: function() {
return this.name || "Recirculation pump";
},
oneditprepare: function() {
$("#node-input-F2").typedInput({
type:"num",
types:["num"]
});
$("#node-input-inlet").typedInput({
type:"num",
types:["num"]
});
},
oneditsave: function() {
let debit = parseFloat($("#node-input-F2").typedInput("value"));
if (isNaN(debit) || debit < 0) {
RED.notify("Debit is not set correctly", {type: "error"});
}
let inlet = parseInt($("#node-input-n_inlets").typedInput("value"));
if (inlet < 1) {
RED.notify("Number of inlets not set correctly", {type: "error"});
}
}
});
</script>
<script type="text/html" data-template-name="recirculation-pump">
<div class="form-row">
<label for="node-input-name"><i class="fa fa-tag"></i> Name</label>
<input type="text" id="node-input-name" placeholder="Name">
</div>
<div class="form-row">
<label for="node-input-F2"><i class="fa fa-tag"></i> Recirculation debit [m3 d-1]</label>
<input type="text" id="node-input-F2" placeholder="m3 s-1">
</div>
<div class="form-row">
<label for="node-input-inlet"><i class="fa fa-tag"></i> Assigned inlet recirculation</label>
<input type="text" id="node-input-inlet" placeholder="#">
</div>
</script>
<script type="text/html" data-help-name="recirculation-pump">
<p>Recirculation-pump for splitting streams</p>
</script>

40
additional_nodes/recirculation-pump.js Normal file
View File

@@ -0,0 +1,40 @@
module.exports = function(RED) {
function recirculation(config) {
RED.nodes.createNode(this, config);
var node = this;
let name = config.name;
let F2 = parseFloat(config.F2);
const inlet_F2 = parseInt(config.inlet);
node.on('input', function(msg, send, done) {
switch (msg.topic) {
case "Fluent":
// conserve volume flow debit
let F_in = msg.payload.F;
let F1 = Math.max(F_in - F2, 0);
let F2_corr = F_in < F2 ? F_in : F2;
let msg_F1 = structuredClone(msg);
msg_F1.payload.F = F1;
let msg_F2 = {...msg};
msg_F2.payload.F = F2_corr;
msg_F2.payload.inlet = inlet_F2;
send([msg_F1, msg_F2]);
break;
case "clock":
break;
default:
console.log("Unknown topic: " + msg.topic);
}
if (done) {
done();
}
});
}
RED.nodes.registerType("recirculation-pump", recirculation);
};

57
additional_nodes/settling-basin.html Normal file
View File

@@ -0,0 +1,57 @@
<script type="text/javascript">
RED.nodes.registerType("settling-basin", {
category: "WWTP",
color: "#e4a363",
defaults: {
name: { value: "" },
TS_set: { value: 0.1, required: true },
inlet: { value: 1, required: true }
},
inputs: 1,
outputs: 2,
outputLabels: ["Main effluent", "Sludge effluent"],
icon: "font-awesome/fa-random",
label: function() {
return this.name || "Settling basin";
},
oneditprepare: function() {
$("#node-input-TS_set").typedInput({
type:"num",
types:["num"]
});
$("#node-input-inlet").typedInput({
type:"num",
types:["num"]
});
},
oneditsave: function() {
let TS_set = parseFloat($("#node-input-TS_set").typedInput("value"));
if (isNaN(TS_set) || TS_set < 0) {
RED.notify("TS is not set correctly", {type: "error"});
}
let inlet = parseInt($("#node-input-n_inlets").typedInput("value"));
if (inlet < 1) {
RED.notify("Number of inlets not set correctly", {type: "error"});
}
}
});
</script>
<script type="text/html" data-template-name="settling-basin">
<div class="form-row">
<label for="node-input-name"><i class="fa fa-tag"></i> Name</label>
<input type="text" id="node-input-name" placeholder="Name">
</div>
<div class="form-row">
<label for="node-input-TS_set"><i class="fa fa-tag"></i> Total Solids set point [g m-3]</label>
<input type="text" id="node-input-TS_set" placeholder="">
</div>
<div class="form-row">
<label for="node-input-inlet"><i class="fa fa-tag"></i> Assigned inlet return line</label>
<input type="text" id="node-input-inlet" placeholder="#">
</div>
</script>
<script type="text/html" data-help-name="settling-basin">
<p>Settling tank</p>
</script>

57
additional_nodes/settling-basin.js Normal file
View File

@@ -0,0 +1,57 @@
module.exports = function(RED) {
function settler(config) {
RED.nodes.createNode(this, config);
var node = this;
let name = config.name;
let TS_set = parseFloat(config.TS_set);
const inlet_sludge = parseInt(config.inlet);
node.on('input', function(msg, send, done) {
switch (msg.topic) {
case "Fluent":
// conserve volume flow debit
let F_in = msg.payload.F;
let C_in = msg.payload.C;
let F2 = (F_in * C_in[12]) / TS_set;
let F1 = Math.max(F_in - F2, 0);
let F2_corr = F_in < F2 ? F_in : F2;
let msg_F1 = structuredClone(msg);
msg_F1.payload.F = F1;
msg_F1.payload.C[7] = 0;
msg_F1.payload.C[8] = 0;
msg_F1.payload.C[9] = 0;
msg_F1.payload.C[10] = 0;
msg_F1.payload.C[11] = 0;
msg_F1.payload.C[12] = 0;
let msg_F2 = {...msg};
msg_F2.payload.F = F2_corr;
if (F2_corr > 0) {
msg_F2.payload.C[7] = F_in * C_in[7] / F2;
msg_F2.payload.C[8] = F_in * C_in[8] / F2;
msg_F2.payload.C[9] = F_in * C_in[9] / F2;
msg_F2.payload.C[10] = F_in * C_in[10] / F2;
msg_F2.payload.C[11] = F_in * C_in[11] / F2;
msg_F2.payload.C[12] = F_in * C_in[12] / F2;
}
msg_F2.payload.inlet = inlet_sludge;
send([msg_F1, msg_F2]);
break;
case "clock":
break;
default:
console.log("Unknown topic: " + msg.topic);
}
if (done) {
done();
}
});
}
RED.nodes.registerType("settling-basin", settler);
};

1763
flows/asm3_flows.json Normal file
View File

File diff suppressed because it is too large Load Diff

119
package-lock.json generated Normal file
View File

@@ -0,0 +1,119 @@
{
"name": "asm3",
"version": "0.0.1",
"lockfileVersion": 3,
"requires": true,
"packages": {
"": {
"name": "asm3",
"version": "0.0.1",
"license": "SEE LICENSE",
"dependencies": {
"generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git#fix-missing-references",
"mathjs": "^14.5.2"
}
},
"node_modules/@babel/runtime": {
"version": "7.28.3",
"resolved": "https://registry.npmjs.org/@babel/runtime/-/runtime-7.28.3.tgz",
"integrity": "sha512-9uIQ10o0WGdpP6GDhXcdOJPJuDgFtIDtN/9+ArJQ2NAfAmiuhTQdzkaTGR33v43GYS2UrSA0eX2pPPHoFVvpxA==",
"license": "MIT",
"engines": {
"node": ">=6.9.0"
}
},
"node_modules/complex.js": {
"version": "2.4.2",
"resolved": "https://registry.npmjs.org/complex.js/-/complex.js-2.4.2.tgz",
"integrity": "sha512-qtx7HRhPGSCBtGiST4/WGHuW+zeaND/6Ld+db6PbrulIB1i2Ev/2UPiqcmpQNPSyfBKraC0EOvOKCB5dGZKt3g==",
"license": "MIT",
"engines": {
"node": "*"
},
"funding": {
"type": "github",
"url": "https://github.com/sponsors/rawify"
}
},
"node_modules/decimal.js": {
"version": "10.6.0",
"resolved": "https://registry.npmjs.org/decimal.js/-/decimal.js-10.6.0.tgz",
"integrity": "sha512-YpgQiITW3JXGntzdUmyUR1V812Hn8T1YVXhCu+wO3OpS4eU9l4YdD3qjyiKdV6mvV29zapkMeD390UVEf2lkUg==",
"license": "MIT"
},
"node_modules/escape-latex": {
"version": "1.2.0",
"resolved": "https://registry.npmjs.org/escape-latex/-/escape-latex-1.2.0.tgz",
"integrity": "sha512-nV5aVWW1K0wEiUIEdZ4erkGGH8mDxGyxSeqPzRNtWP7ataw+/olFObw7hujFWlVjNsaDFw5VZ5NzVSIqRgfTiw==",
"license": "MIT"
},
"node_modules/fraction.js": {
"version": "5.3.4",
"resolved": "https://registry.npmjs.org/fraction.js/-/fraction.js-5.3.4.tgz",
"integrity": "sha512-1X1NTtiJphryn/uLQz3whtY6jK3fTqoE3ohKs0tT+Ujr1W59oopxmoEh7Lu5p6vBaPbgoM0bzveAW4Qi5RyWDQ==",
"license": "MIT",
"engines": {
"node": "*"
},
"funding": {
"type": "github",
"url": "https://github.com/sponsors/rawify"
}
},
"node_modules/generalFunctions": {
"version": "1.0.0",
"resolved": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git#302e12238745766a679ef11ca6ed5f4ea1548f87",
"license": "SEE LICENSE"
},
"node_modules/javascript-natural-sort": {
"version": "0.7.1",
"resolved": "https://registry.npmjs.org/javascript-natural-sort/-/javascript-natural-sort-0.7.1.tgz",
"integrity": "sha512-nO6jcEfZWQXDhOiBtG2KvKyEptz7RVbpGP4vTD2hLBdmNQSsCiicO2Ioinv6UI4y9ukqnBpy+XZ9H6uLNgJTlw==",
"license": "MIT"
},
"node_modules/mathjs": {
"version": "14.7.0",
"resolved": "https://registry.npmjs.org/mathjs/-/mathjs-14.7.0.tgz",
"integrity": "sha512-RaMhb+9MSESjDZNox/FzzuFpIUI+oxGLyOy1t3BMoW53pGWnTzZtlucJ5cvbit0dIMYlCq00gNbW1giZX4/1Rg==",
"license": "Apache-2.0",
"dependencies": {
"@babel/runtime": "^7.26.10",
"complex.js": "^2.2.5",
"decimal.js": "^10.4.3",
"escape-latex": "^1.2.0",
"fraction.js": "^5.2.1",
"javascript-natural-sort": "^0.7.1",
"seedrandom": "^3.0.5",
"tiny-emitter": "^2.1.0",
"typed-function": "^4.2.1"
},
"bin": {
"mathjs": "bin/cli.js"
},
"engines": {
"node": ">= 18"
}
},
"node_modules/seedrandom": {
"version": "3.0.5",
"resolved": "https://registry.npmjs.org/seedrandom/-/seedrandom-3.0.5.tgz",
"integrity": "sha512-8OwmbklUNzwezjGInmZ+2clQmExQPvomqjL7LFqOYqtmuxRgQYqOD3mHaU+MvZn5FLUeVxVfQjwLZW/n/JFuqg==",
"license": "MIT"
},
"node_modules/tiny-emitter": {
"version": "2.1.0",
"resolved": "https://registry.npmjs.org/tiny-emitter/-/tiny-emitter-2.1.0.tgz",
"integrity": "sha512-NB6Dk1A9xgQPMoGqC5CVXn123gWyte215ONT5Pp5a0yt4nlEoO1ZWeCwpncaekPHXO60i47ihFnZPiRPjRMq4Q==",
"license": "MIT"
},
"node_modules/typed-function": {
"version": "4.2.1",
"resolved": "https://registry.npmjs.org/typed-function/-/typed-function-4.2.1.tgz",
"integrity": "sha512-EGjWssW7Tsk4DGfE+5yluuljS1OGYWiI1J6e8puZz9nTMM51Oug8CD5Zo4gWMsOhq5BI+1bF+rWTm4Vbj3ivRA==",
"license": "MIT",
"engines": {
"node": ">= 18"
}
}
}
}

33
package.json Normal file
View File

@@ -0,0 +1,33 @@
{
"name": "asm3",
"version": "0.0.1",
"description": "Implementation of the asm3 model for Node-Red",
"repository": {
"type": "git",
"url": "https://gitea.centraal.wbd-rd.nl/RnD/reactor.git"
},
"keywords": [
"asm3",
"activated sludge",
"wastewater",
"biological model",
"node-red"
],
"license": "SEE LICENSE",
"author": "P.R. van der Wilt",
"main": "reactor.js",
"scripts": {
"test": "node reactor.js"
},
"node-red": {
"nodes": {
"reactor": "reactor.js",
"recirculation-pump": "additional_nodes/recirculation-pump.js",
"settling-basin": "additional_nodes/settling-basin.js"
}
},
"dependencies": {
"generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/p.vanderwilt/generalFunctions.git#fix-missing-references",
"mathjs": "^14.5.2"
}
}

248
reactor.html Normal file
View File

@@ -0,0 +1,248 @@
<script src="/advancedReactor/menu.js"></script>
<script type="text/javascript">
RED.nodes.registerType("advancedReactor", {
category: "WWTP",
color: "#c4cce0",
defaults: {
name: { value: "" },
reactor_type: { value: "CSTR", required: true },
volume: { value: 0., required: true },
length: { value: 0.},
resolution_L: { value: 0.},
alpha: {value: 0},
n_inlets: { value: 1, required: true},
kla: { value: null },
S_O_init: { value: 0., required: true },
S_I_init: { value: 30., required: true },
S_S_init: { value: 100., required: true },
S_NH_init: { value: 16., required: true },
S_N2_init: { value: 0., required: true },
S_NO_init: { value: 0., required: true },
S_HCO_init: { value: 5., required: true },
X_I_init: { value: 25., required: true },
X_S_init: { value: 75., required: true },
X_H_init: { value: 30., required: true },
X_STO_init: { value: 0., required: true },
X_A_init: { value: 0.001, required: true },
X_TS_init: { value: 125.0009, required: true },
timeStep: { value: 1, required: true },
enableLog: { value: false },
logLevel: { value: "error" },
positionVsParent: { value: "" },
},
inputs: 1,
outputs: 3,
inputLabels: ["input"],
outputLabels: ["process", "dbase", "parent"],
icon: "font-awesome/fa-recycle",
label: function() {
return this.name || "advancedReactor";
},
oneditprepare: function() {
// wait for the menu scripts to load
const waitForMenuData = () => {
if (window.EVOLV?.nodes?.advancedReactor?.initEditor) {
window.EVOLV.nodes.advancedReactor.initEditor(this);
} else {
setTimeout(waitForMenuData, 50);
}
};
waitForMenuData();
$("#node-input-volume").typedInput({
type:"num",
types:["num"]
});
$("#node-input-n_inlets").typedInput({
type:"num",
types:["num"]
});
$("#node-input-length").typedInput({
type:"num",
types:["num"]
});
$("#node-input-resolution_L").typedInput({
type:"num",
types:["num"]
});
$("#node-input-kla").typedInput({
type:"num",
types:["num"]
});
$(".concentrations").typedInput({
type:"num",
types:["num"]
});
$("#node-input-reactor_type").typedInput({
types: [
{
value: "CSTR",
options: [
{ value: "CSTR", label: "CSTR"},
{ value: "PFR", label: "PFR"}
]
}
]
})
$("#node-input-reactor_type").on("change", function() {
const type = $("#node-input-reactor_type").typedInput("value");
if (type === "CSTR") {
$(".PFR").hide();
} else {
$(".PFR").show();
}
});
$("#node-input-alpha").typedInput({
type:"num",
types:["num"]
})
$("#node-input-timeStep").typedInput({
type:"num",
types:["num"]
})
// Set initial visibility on dialog open
const initialType = $("#node-input-reactor_type").typedInput("value");
if (initialType === "CSTR") {
$(".PFR").hide();
} else {
$(".PFR").show();
}
},
oneditsave: function() {
// save logger fields
if (window.EVOLV?.nodes?.['advancedReactor']?.loggerMenu?.saveEditor) {
window.EVOLV.nodes['advancedReactor'].loggerMenu.saveEditor(this);
}
// save position field
if (window.EVOLV?.nodes?.measurement?.positionMenu?.saveEditor) {
window.EVOLV.nodes.rotatingMachine.positionMenu.saveEditor(this);
}
let volume = parseFloat($("#node-input-volume").typedInput("value"));
if (isNaN(volume) || volume <= 0) {
RED.notify("Fluid volume not set correctly", {type: "error"});
}
let n_inlets = parseInt($("#node-input-n_inlets").typedInput("value"));
if (isNaN(n_inlets) || n_inlets < 1) {
RED.notify("Number of inlets not set correctly", {type: "error"});
}
}
});
</script>
<script type="text/html" data-template-name="advancedReactor">
<div class="form-row">
<label for="node-input-name"><i class="fa fa-tag"></i> Name</label>
<input type="text" id="node-input-name" placeholder="Name">
</div>
<h2> Reactor properties </h2>
<div class="form-row">
<label for="node-input-reactor_type"><i class="fa fa-tag"></i> Reactor type</label>
<input type="text" id="node-input-reactor_type">
</div>
<div class="form-row">
<label for="node-input-volume"><i class="fa fa-tag"></i> Fluid volume [m3]</label>
<input type="text" id="node-input-volume" placeholder="m3">
</div>
<div class="form-row PFR">
<label for="node-input-length"><i class="fa fa-tag"></i> Reactor length [m]</label>
<input type="text" id="node-input-length" placeholder="m">
</div>
<div class="form-row PFR">
<label for="node-input-resolution_L"><i class="fa fa-tag"></i> Resolution</label>
<input type="text" id="node-input-resolution_L" placeholder="#">
</div>
<div class="PFR">
<p> Inlet boundary condition parameter &alpha; (&alpha; = 0: Danckwerts BC / &alpha; = 1: Dirichlet BC) </p>
<div class="form-row">
<label for="node-input-alpha"><i class="fa fa-tag"></i>Adjustable parameter BC</label>
<input type="text" id="node-input-alpha">
</div>
</div>
<div class="form-row">
<label for="node-input-n_inlets"><i class="fa fa-tag"></i> Number of inlets</label>
<input type="text" id="node-input-n_inlets" placeholder="#">
</div>
<h3> Internal mass transfer calculation (optional) </h3>
<div class="form-row">
<label for="node-input-kla"><i class="fa fa-tag"></i> kLa [d-1]</label>
<input type="text" id="node-input-kla" placeholder="d-1">
</div>
<h2> Dissolved components </h2>
<div class="form-row">
<label for="node-input-S_O_init"><i class="fa fa-tag"></i> Initial dissolved oxygen [g O2 m-3]</label>
<input type="text" id="node-input-S_O_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-S_I_init"><i class="fa fa-tag"></i> Initial soluble inert organics [g COD m-3]</label>
<input type="text" id="node-input-S_I_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-S_S_init"><i class="fa fa-tag"></i> Initial readily biodegrable substrates [g COD m-3]</label>
<input type="text" id="node-input-S_S_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-S_NH_init"><i class="fa fa-tag"></i> Initial ammonium / ammonia [g N m-3]</label>
<input type="text" id="node-input-S_NH_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-S_N2_init"><i class="fa fa-tag"></i> Initial dinitrogen, released by denitrification [g N m-3]</label>
<input type="text" id="node-input-S_N2_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-S_NO_init"><i class="fa fa-tag"></i> Initial nitrite + nitrate [g N m-3]</label>
<input type="text" id="node-input-S_NO_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-S_HCO_init"><i class="fa fa-tag"></i> Initial alkalinity, bicarbonate [mole HCO3- m-3]</label>
<input type="text" id="node-input-S_HCO_init" class="concentrations">
</div>
<h2> Particulate components </h2>
<div class="form-row">
<label for="node-input-X_I_init"><i class="fa fa-tag"></i> Initial inert particulate organics [g COD m-3]</label>
<input type="text" id="node-input-X_I_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-X_S_init"><i class="fa fa-tag"></i> Initial slowly biodegrable substrates [g COD m-3]</label>
<input type="text" id="node-input-X_S_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-X_H_init"><i class="fa fa-tag"></i> Initial heterotrophic biomass [g COD m-3]</label>
<input type="text" id="node-input-X_H_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-X_STO_init"><i class="fa fa-tag"></i> Initial Organics stored by heterotrophs [g COD m-3]</label>
<input type="text" id="node-input-X_STO_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-X_A_init"><i class="fa fa-tag"></i> Initial autotrophic, nitrifying biomass [g COD m-3]</label>
<input type="text" id="node-input-X_A_init" class="concentrations">
</div>
<div class="form-row">
<label for="node-input-X_TS_init"><i class="fa fa-tag"></i> Initial total suspended solids [g TSS m-3]</label>
<input type="text" id="node-input-X_TS_init" class="concentrations">
</div>
<h2> Simulation parameters </h2>
<div class="form-row">
<label for="node-input-timeStep"><i class="fa fa-tag"></i> Time step [s]</label>
<input type="text" id="node-input-timeStep" placeholder="s">
</div>
<!-- Logger fields injected here -->
<div id="logger-fields-placeholder"></div>
<!-- Position fields will be injected here -->
<div id="position-fields-placeholder"></div>
</script>
<script type="text/html" data-help-name="advancedReactor">
<p>New reactor node</p>
</script>

26
reactor.js Normal file
View File

@@ -0,0 +1,26 @@
const nameOfNode = "advancedReactor"; // name of the node, should match file name and node type in Node-RED
const nodeClass = require('./src/nodeClass.js'); // node class
const { MenuManager } = require('generalFunctions');
module.exports = function (RED) {
// Register the node type
RED.nodes.registerType(nameOfNode, function (config) {
// Initialize the Node-RED node first
RED.nodes.createNode(this, config);
// Then create your custom class and attach it
this.nodeClass = new nodeClass(config, RED, this, nameOfNode);
});
const menuMgr = new MenuManager();
// Serve /advancedReactor/menu.js
RED.httpAdmin.get(`/${nameOfNode}/menu.js`, (req, res) => {
try {
const script = menuMgr.createEndpoint(nameOfNode, ['logger', 'position']);
res.type('application/javascript').send(script);
} catch (err) {
res.status(500).send(`// Error generating menu: ${err.message}`);
}
});
};

165
src/nodeClass.js Normal file
View File

@@ -0,0 +1,165 @@
const { Reactor_CSTR, Reactor_PFR } = require('./specificClass.js');
class nodeClass {
/**
* Node-RED node class for advanced-reactor.
* @param {object} uiConfig - Node-RED node configuration
* @param {object} RED - Node-RED runtime API
* @param {object} nodeInstance - Node-RED node instance
* @param {string} nameOfNode - Name of the node
*/
constructor(uiConfig, RED, nodeInstance, nameOfNode) {
// Preserve RED reference for HTTP endpoints if needed
this.node = nodeInstance;
this.RED = RED;
this.name = nameOfNode;
this.source = null;
this._loadConfig(uiConfig)
this._setupClass();
this._attachInputHandler();
this._registerChild();
this._startTickLoop();
this._attachCloseHandler();
}
/**
* Handle node-red input messages
*/
_attachInputHandler() {
this.node.on('input', (msg, send, done) => {
switch (msg.topic) {
case "clock":
this.source.updateState(msg.timestamp);
send([msg, null, null]);
break;
case "Fluent":
this.source.setInfluent = msg;
break;
case "OTR":
this.source.setOTR = msg;
break;
case "Temperature":
this.source.setTemperature = msg;
break;
case "Dispersion":
this.source.setDispersion = msg;
break;
case 'registerChild':
// Register this node as a parent of the child node
const childId = msg.payload;
const childObj = this.RED.nodes.getNode(childId);
this.source.childRegistrationUtils.registerChild(childObj.source, msg.positionVsParent);
break;
default:
console.log("Unknown topic: " + msg.topic);
}
if (done) {
done();
}
});
}
/**
* Parse node configuration
* @param {object} uiConfig Config set in UI in node-red
*/
_loadConfig(uiConfig) {
this.config = {
general: {
name: uiConfig.name || this.name,
id: this.node.id,
unit: null,
logging: {
enabled: uiConfig.enableLog,
logLevel: uiConfig.logLevel
}
},
functionality: {
positionVsParent: uiConfig.positionVsParent || 'atEquipment', // Default to 'atEquipment' if not specified
softwareType: "reactor" // should be set in config manager
},
reactor_type: uiConfig.reactor_type,
volume: parseFloat(uiConfig.volume),
length: parseFloat(uiConfig.length),
resolution_L: parseInt(uiConfig.resolution_L),
alpha: parseFloat(uiConfig.alpha),
n_inlets: parseInt(uiConfig.n_inlets),
kla: parseFloat(uiConfig.kla),
initialState: [
parseFloat(uiConfig.S_O_init),
parseFloat(uiConfig.S_I_init),
parseFloat(uiConfig.S_S_init),
parseFloat(uiConfig.S_NH_init),
parseFloat(uiConfig.S_N2_init),
parseFloat(uiConfig.S_NO_init),
parseFloat(uiConfig.S_HCO_init),
parseFloat(uiConfig.X_I_init),
parseFloat(uiConfig.X_S_init),
parseFloat(uiConfig.X_H_init),
parseFloat(uiConfig.X_STO_init),
parseFloat(uiConfig.X_A_init),
parseFloat(uiConfig.X_TS_init)
],
timeStep: parseFloat(uiConfig.timeStep)
}
}
/**
* Register this node as a child upstream and downstream.
* Delayed to avoid Node-RED startup race conditions.
*/
_registerChild() {
setTimeout(() => {
this.node.send([
null,
null,
{ topic: 'registerChild', payload: this.node.id, positionVsParent: this.config?.functionality?.positionVsParent || 'atEquipment' }
]);
}, 100);
}
/**
* Setup reactor class based on config
*/
_setupClass() {
let new_reactor;
switch (this.config.reactor_type) {
case "CSTR":
new_reactor = new Reactor_CSTR(this.config);
break;
case "PFR":
new_reactor = new Reactor_PFR(this.config);
break;
default:
console.warn("Unknown reactor type: " + uiConfig.reactor_type);
}
this.source = new_reactor; // protect from reassignment
this.node.source = this.source;
}
_startTickLoop() {
setTimeout(() => {
this._tickInterval = setInterval(() => this._tick(), 1000);
}, 1000);
}
_tick(){
this.node.send([this.source.getEffluent, null, null]);
}
_attachCloseHandler() {
this.node.on('close', (done) => {
clearInterval(this._tickInterval);
done();
});
}
}
module.exports = nodeClass;

211
src/reaction_modules/asm3_class Koch.js Normal file
View File

@@ -0,0 +1,211 @@
const math = require('mathjs')
/**
* ASM3 class for the Activated Sludge Model No. 3 (ASM3). Using Koch et al. 2000 parameters.
*/
class ASM3 {
constructor() {
/**
* 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]
};
/**
* 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]
};
/**
* Temperature theta parameters for ASM3. Using Koch et al. 2000 parameters.
* These parameters are used to adjust reaction rates based on temperature.
* @property {Object} temp_params - Temperature theta parameters
*/
this.temp_params = {
// Hydrolysis
theta_H: 0.04,
// Heterotrophs
theta_STO: 0.07,
theta_mu_H: 0.07,
theta_b_H_O: 0.07,
theta_b_H_NO: 0.07,
theta_b_STO_O: this._compute_theta(0.1, 0.3, 10, 20),
theta_b_STO_NO: this._compute_theta(0.05, 0.15, 10, 20),
// Autotrophs
theta_mu_A: 0.105,
theta_b_A_O: 0.105,
theta_b_A_NO: 0.105
};
this.stoi_matrix = this._initialise_stoi_matrix();
}
/**
* Initialises the stoichiometric matrix for ASM3.
* @returns {Array} - The stoichiometric matrix for ASM3. (2D array)
*/
_initialise_stoi_matrix() { // initialise stoichiometric matrix
const { f_SI, f_XI, Y_STO_O, Y_STO_NO, Y_H_O, Y_H_NO, Y_A, i_CODN, i_CODNO, i_NSI, i_NSS, i_NXI, i_NXS, i_NBM, i_TSXI, i_TSXS, i_TSBM, i_TSSTO, i_cNH, i_cNO } = this.stoi_params;
const stoi_matrix = Array(12);
// S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
stoi_matrix[0] = [0., f_SI, 1.-f_SI, i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI, 0., 0., (i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI)*i_cNH, 0., -1., 0., 0., 0., -i_TSXS];
stoi_matrix[1] = [-(1.-Y_STO_O), 0, -1., i_NSS, 0., 0., i_NSS*i_cNH, 0., 0., 0., Y_STO_O, 0., Y_STO_O*i_TSSTO];
stoi_matrix[2] = [0., 0., -1., i_NSS, -(1.-Y_STO_NO)/(i_CODNO-i_CODN), (1.-Y_STO_NO)/(i_CODNO-i_CODN), i_NSS*i_cNH + (1.-Y_STO_NO)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., Y_STO_NO, 0., Y_STO_NO*i_TSSTO];
stoi_matrix[3] = [-(1.-Y_H_O)/Y_H_O, 0., 0., -i_NBM, 0., 0., -i_NBM*i_cNH, 0., 0., 1., -1./Y_H_O, 0., i_TSBM-i_TSSTO/Y_H_O];
stoi_matrix[4] = [0., 0., 0., -i_NBM, -(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), (1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), -i_NBM*i_cNH+(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN))*i_cNO, 0., 0., 1., -1./Y_H_NO, 0., i_TSBM-i_TSSTO/Y_H_NO];
stoi_matrix[5] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
stoi_matrix[6] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
stoi_matrix[7] = [-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., -i_TSSTO];
stoi_matrix[8] = [0., 0., 0., 0., -1./(i_CODNO-i_CODN), 1./(i_CODNO-i_CODN), i_cNO/(i_CODNO-i_CODN), 0., 0., 0., -1., 0., -i_TSSTO];
stoi_matrix[9] = [1.+i_CODNO/Y_A, 0., 0., -1./Y_A-i_NBM, 0., 1./Y_A, (-1./Y_A-i_NBM)*i_cNH+i_cNO/Y_A, 0., 0., 0., 0., 1., i_TSBM];
stoi_matrix[10] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
stoi_matrix[11] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
return stoi_matrix[0].map((col, i) => stoi_matrix.map(row => row[i])); // transpose matrix
}
/**
* Computes the Monod equation rate value for a given concentration and half-saturation constant.
* @param {number} c - Concentration of reaction species.
* @param {number} K - Half-saturation constant for the reaction species.
* @returns {number} - Monod equation rate value for the given concentration and half-saturation constant.
*/
_monod(c, K) {
return c / (K + c);
}
/**
* Computes the inverse Monod equation rate value for a given concentration and half-saturation constant. Used for inhibition.
* @param {number} c - Concentration of reaction species.
* @param {number} K - Half-saturation constant for the reaction species.
* @returns {number} - Inverse Monod equation rate value for the given concentration and half-saturation constant.
*/
_inv_monod(c, K) {
return K / (K + c);
}
/**
* Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter.
* @param {number} k - Rate constant at 20 degrees Celcius.
* @param {number} theta - Theta parameter.
* @param {number} T - Temperature in Celcius.
* @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation.
*/
_arrhenius(k, theta, T) {
return k * Math.exp(theta*(T-20));
}
/**
* Computes the temperature theta parameter based on two rate constants and their corresponding temperatures.
* @param {number} k1 - Rate constant at temperature T1.
* @param {number} k2 - Rate constant at temperature T2.
* @param {number} T1 - Temperature T1 in Celcius.
* @param {number} T2 - Temperature T2 in Celcius.
* @returns {number} - Theta parameter.
*/
_compute_theta(k1, k2, T1, T2) {
return Math.log(k1/k2)/(T1-T2);
}
/**
* Computes the reaction rates for each process reaction based on the current state and temperature.
* @param {Array} state - State vector containing concentrations of reaction species.
* @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
* @returns {Array} - Reaction rates for each process reaction.
*/
compute_rates(state, T = 20) {
// state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
const rates = Array(12);
const [S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS] = state;
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;
const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params;
// Hydrolysis
rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * this._monod(X_S / X_H, K_X) * X_H;
// Heterotrophs
rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
rates[2] = this._arrhenius(k_STO, theta_STO, T) * 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] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
// Autotrophs
rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
return rates;
}
/**
* Computes the change in concentrations of reaction species based on the current state and temperature.
* @param {Array} state - State vector containing concentrations of reaction species.
* @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
* @returns {Array} - Change in reaction species concentrations.
*/
compute_dC(state, T = 20) { // compute changes in concentrations
// state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
return math.multiply(this.stoi_matrix, this.compute_rates(state, T));
}
}
module.exports = ASM3;

211
src/reaction_modules/asm3_class.js Normal file
View File

@@ -0,0 +1,211 @@
const math = require('mathjs')
/**
* ASM3 class for the Activated Sludge Model No. 3 (ASM3).
*/
class ASM3 {
constructor() {
/**
* 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]
};
/**
* 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]
};
/**
* Temperature theta parameters for ASM3.
* These parameters are used to adjust reaction rates based on temperature.
* @property {Object} temp_params - Temperature theta parameters
*/
this.temp_params = {
// Hydrolysis
theta_H: this._compute_theta(2, 3, 10, 20),
// Heterotrophs
theta_STO: this._compute_theta(2.5, 5, 10, 20),
theta_mu_H: this._compute_theta(1, 2, 10, 20),
theta_b_H_O: this._compute_theta(0.1, 0.2, 10, 20),
theta_b_H_NO: this._compute_theta(0.05, 0.1, 10, 20),
theta_b_STO_O: this._compute_theta(0.1, 0.2, 10, 20),
theta_b_STO_NO: this._compute_theta(0.05, 0.1, 10, 20),
// Autotrophs
theta_mu_A: this._compute_theta(0.35, 1, 10, 20),
theta_b_A_O: this._compute_theta(0.05, 0.15, 10, 20),
theta_b_A_NO: this._compute_theta(0.02, 0.05, 10, 20)
};
this.stoi_matrix = this._initialise_stoi_matrix();
}
/**
* Initialises the stoichiometric matrix for ASM3.
* @returns {Array} - The stoichiometric matrix for ASM3. (2D array)
*/
_initialise_stoi_matrix() { // initialise stoichiometric matrix
const { f_SI, f_XI, Y_STO_O, Y_STO_NO, Y_H_O, Y_H_NO, Y_A, i_CODN, i_CODNO, i_NSI, i_NSS, i_NXI, i_NXS, i_NBM, i_TSXI, i_TSXS, i_TSBM, i_TSSTO, i_cNH, i_cNO } = this.stoi_params;
const stoi_matrix = Array(12);
// S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
stoi_matrix[0] = [0., f_SI, 1.-f_SI, i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI, 0., 0., (i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI)*i_cNH, 0., -1., 0., 0., 0., -i_TSXS];
stoi_matrix[1] = [-(1.-Y_STO_O), 0, -1., i_NSS, 0., 0., i_NSS*i_cNH, 0., 0., 0., Y_STO_O, 0., Y_STO_O*i_TSSTO];
stoi_matrix[2] = [0., 0., -1., i_NSS, -(1.-Y_STO_NO)/(i_CODNO-i_CODN), (1.-Y_STO_NO)/(i_CODNO-i_CODN), i_NSS*i_cNH + (1.-Y_STO_NO)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., Y_STO_NO, 0., Y_STO_NO*i_TSSTO];
stoi_matrix[3] = [-(1.-Y_H_O)/Y_H_O, 0., 0., -i_NBM, 0., 0., -i_NBM*i_cNH, 0., 0., 1., -1./Y_H_O, 0., i_TSBM-i_TSSTO/Y_H_O];
stoi_matrix[4] = [0., 0., 0., -i_NBM, -(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), (1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), -i_NBM*i_cNH+(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN))*i_cNO, 0., 0., 1., -1./Y_H_NO, 0., i_TSBM-i_TSSTO/Y_H_NO];
stoi_matrix[5] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
stoi_matrix[6] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
stoi_matrix[7] = [-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., -i_TSSTO];
stoi_matrix[8] = [0., 0., 0., 0., -1./(i_CODNO-i_CODN), 1./(i_CODNO-i_CODN), i_cNO/(i_CODNO-i_CODN), 0., 0., 0., -1., 0., -i_TSSTO];
stoi_matrix[9] = [1.+i_CODNO/Y_A, 0., 0., -1./Y_A-i_NBM, 0., 1./Y_A, (-1./Y_A-i_NBM)*i_cNH+i_cNO/Y_A, 0., 0., 0., 0., 1., i_TSBM];
stoi_matrix[10] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
stoi_matrix[11] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
return stoi_matrix[0].map((col, i) => stoi_matrix.map(row => row[i])); // transpose matrix
}
/**
* Computes the Monod equation rate value for a given concentration and half-saturation constant.
* @param {number} c - Concentration of reaction species.
* @param {number} K - Half-saturation constant for the reaction species.
* @returns {number} - Monod equation rate value for the given concentration and half-saturation constant.
*/
_monod(c, K) {
return c / (K + c);
}
/**
* Computes the inverse Monod equation rate value for a given concentration and half-saturation constant. Used for inhibition.
* @param {number} c - Concentration of reaction species.
* @param {number} K - Half-saturation constant for the reaction species.
* @returns {number} - Inverse Monod equation rate value for the given concentration and half-saturation constant.
*/
_inv_monod(c, K) {
return K / (K + c);
}
/**
* Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter.
* @param {number} k - Rate constant at 20 degrees Celcius.
* @param {number} theta - Theta parameter.
* @param {number} T - Temperature in Celcius.
* @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation.
*/
_arrhenius(k, theta, T) {
return k * Math.exp(theta*(T-20));
}
/**
* Computes the temperature theta parameter based on two rate constants and their corresponding temperatures.
* @param {number} k1 - Rate constant at temperature T1.
* @param {number} k2 - Rate constant at temperature T2.
* @param {number} T1 - Temperature T1 in Celcius.
* @param {number} T2 - Temperature T2 in Celcius.
* @returns {number} - Theta parameter.
*/
_compute_theta(k1, k2, T1, T2) {
return Math.log(k1/k2)/(T1-T2);
}
/**
* Computes the reaction rates for each process reaction based on the current state and temperature.
* @param {Array} state - State vector containing concentrations of reaction species.
* @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
* @returns {Array} - Reaction rates for each process reaction.
*/
compute_rates(state, T = 20) {
// state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
const rates = Array(12);
const [S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS] = state;
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;
const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params;
// Hydrolysis
rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * this._monod(X_S / X_H, K_X) * X_H;
// Heterotrophs
rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
rates[2] = this._arrhenius(k_STO, theta_STO, T) * 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] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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 : this._arrhenius(mu_H_max, theta_mu_H, T) * 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] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
// Autotrophs
rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
return rates;
}
/**
* Computes the change in concentrations of reaction species based on the current state and temperature.
* @param {Array} state - State vector containing concentrations of reaction species.
* @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
* @returns {Array} - Change in reaction species concentrations.
*/
compute_dC(state, T = 20) { // compute changes in concentrations
// state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
return math.multiply(this.stoi_matrix, this.compute_rates(state, T));
}
}
module.exports = ASM3;

413
src/specificClass.js Normal file
View File

@@ -0,0 +1,413 @@
const ASM3 = require('./reaction_modules/asm3_class.js');
const { create, all, isArray } = require('mathjs');
const { assertNoNaN } = require('./utils.js');
const { childRegistrationUtils, logger, MeasurementContainer } = require('generalFunctions');
const EventEmitter = require('events');
const mathConfig = {
matrix: 'Array' // use Array as the matrix type
};
const math = create(all, mathConfig);
const S_O_INDEX = 0;
const NUM_SPECIES = 13;
const DEBUG = false;
class Reactor {
/**
* Reactor base class.
* @param {object} config - Configuration object containing reactor parameters.
*/
constructor(config) {
this.config = config;
// EVOLV stuff
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.upstreamReactor = null;
this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
this.asm = new ASM3();
this.volume = config.volume; // fluid volume reactor [m3]
this.Fs = Array(config.n_inlets).fill(0); // fluid debits per inlet [m3 d-1]
this.Cs_in = Array.from(Array(config.n_inlets), () => new Array(NUM_SPECIES).fill(0)); // composition influents
this.OTR = 0.0; // oxygen transfer rate [g O2 d-1 m-3]
this.temperature = 20; // temperature [C]
this.kla = config.kla; // if NaN, use externaly provided OTR [d-1]
this.currentTime = Date.now(); // milliseconds since epoch [ms]
this.timeStep = 1 / (24*60*60) * this.config.timeStep; // time step [d]
this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second
}
/**
* Setter for influent data.
* @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations.
*/
set setInfluent(input) {
let index_in = input.payload.inlet;
this.Fs[index_in] = input.payload.F;
this.Cs_in[index_in] = input.payload.C;
}
/**
* Setter for OTR (Oxygen Transfer Rate).
* @param {object} input - Input object (msg) containing payload with OTR value [g O2 d-1 m-3].
*/
set setOTR(input) {
this.OTR = input.payload;
}
/**
* Getter for effluent data.
* @returns {object} Effluent data object (msg), defaults to inlet 0.
*/
get getEffluent() { // getter for Effluent, defaults to inlet 0
if (isArray(this.state.at(-1))) {
return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state.at(-1) }, timestamp: this.currentTime };
}
return { topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: this.state }, timestamp: this.currentTime };
}
/**
* Calculate the oxygen transfer rate (OTR) based on the dissolved oxygen concentration and temperature.
* @param {number} S_O - Dissolved oxygen concentration [g O2 m-3].
* @param {number} T - Temperature in Celsius, default to 20 C.
* @returns {number} - Calculated OTR [g O2 d-1 m-3].
*/
_calcOTR(S_O, T = 20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
let S_O_sat = 14.652 - 4.1022e-1 * T + 7.9910e-3 * T*T + 7.7774e-5 * T*T*T;
return this.kla * (S_O_sat - S_O);
}
/**
* Clip values in an array to zero.
* @param {Array} arr - Array of values to clip.
* @returns {Array} - New array with values clipped to zero.
*/
_arrayClip2Zero(arr) {
if (Array.isArray(arr)) {
return arr.map(x => this._arrayClip2Zero(x));
} else {
return arr < 0 ? 0 : arr;
}
}
registerChild(child, softwareType) {
switch (softwareType) {
case "measurement":
this.logger.debug(`Registering measurement child.`);
this._connectMeasurement(child);
break;
case "reactor":
this.logger.debug(`Registering reactor child.`);
this._connectReactor(child);
break;
default:
this.logger.error(`Unrecognized softwareType: ${softwareType}`);
}
}
_connectMeasurement(measurement) {
if (!measurement) {
this.logger.warn("Invalid measurement provided.");
return;
}
const position = measurement.config.functionality.positionVsParent;
const measurementType = measurement.config.asset.type;
const key = `${measurementType}_${position}`;
const eventName = `${measurementType}.measured.${position}`;
// Register event listener for measurement updates
this.measurements.emitter.on(eventName, (eventData) => {
this.logger.debug(`${position} ${measurementType} from ${eventData.childName}: ${eventData.value} ${eventData.unit}`);
// Store directly in parent's measurement container
this.measurements
.type(measurementType)
.variant("measured")
.position(position)
.value(eventData.value, eventData.timestamp, eventData.unit);
this._updateMeasurement(measurementType, eventData.value, position, eventData);
});
}
_connectReactor(reactor) {
if (!reactor) {
this.logger.warn("Invalid reactor provided.");
return;
}
this.upstreamReactor = reactor;
reactor.emitter.on("stateChange", (data) => {
this.logger.debug(`State change of upstream reactor detected.`);
this.updateState(data);
});
}
_updateMeasurement(measurementType, value, position, context) {
this.logger.debug(`---------------------- updating ${measurementType} ------------------ `);
switch (measurementType) {
case "temperature":
if (position == "atEquipment") {
this.temperature = value;
}
break;
default:
this.logger.error(`Type '${measurementType}' not recognized for measured update.`);
return;
}
}
/**
* 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 = Date.now()) { // expect update with timestamp
const day2ms = 1000 * 60 * 60 * 24;
if (this.upstreamReactor) {
this.setInfluent = this.upstreamReactor.getEffluent;
}
let 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);
}
}
}
class Reactor_CSTR extends Reactor {
/**
* Reactor_CSTR class for Continuous Stirred Tank Reactor.
* @param {object} config - Configuration object containing reactor parameters.
*/
constructor(config) {
super(config);
this.state = config.initialState;
}
/**
* Tick the reactor state using the forward Euler method.
* @param {number} time_step - Time step for the simulation [d].
* @returns {Array} - New reactor state.
*/
tick(time_step) { // tick reactor state using forward Euler method
const inflow = math.multiply(math.divide([this.Fs], this.volume), this.Cs_in)[0];
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(NUM_SPECIES).fill(0.0);
transfer[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
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
if(DEBUG){
assertNoNaN(dC_total, "change in state");
assertNoNaN(this.state, "new state");
}
return this.state;
}
}
class Reactor_PFR extends Reactor {
/**
* Reactor_PFR class for Plug Flow Reactor.
* @param {object} config - Configuration object containing reactor parameters.
*/
constructor(config) {
super(config);
this.length = config.length; // reactor length [m]
this.n_x = config.resolution_L; // number of slices
this.d_x = this.length / this.n_x;
this.A = this.volume / this.length; // crosssectional area [m2]
this.alpha = config.alpha;
this.state = Array.from(Array(this.n_x), () => config.initialState.slice())
this.D = 0.0; // axial dispersion [m2 d-1]
this.D_op = this._makeDoperator(true, true);
assertNoNaN(this.D_op, "Derivative operator");
this.D2_op = this._makeD2operator();
assertNoNaN(this.D2_op, "Second derivative operator");
}
/**
* Setter for axial dispersion.
* @param {object} input - Input object (msg) containing payload with dispersion value [m2 d-1].
*/
set setDispersion(input) {
this.D = input.payload;
}
updateState(newTime) {
super.updateState(newTime);
let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
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.`);
(Co_D >= 0.5) && this.logger.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`);
if(DEBUG) {
console.log("Inlet state max " + math.max(this.state[0]))
console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A));
console.log("Pe local " + Pe_local);
console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
console.log("Co D " + Co_D);
}
}
/**
* Tick the reactor state using explicit finite difference method.
* @param {number} time_step - Time step for the simulation [d].
* @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));
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++) {
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);
}
}
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
const stateNew = math.add(this.state, dC_total);
this._applyBoundaryConditions(stateNew);
if (DEBUG) {
assertNoNaN(dispersion, "dispersion");
assertNoNaN(advection, "advection");
assertNoNaN(reaction, "reaction");
assertNoNaN(dC_total, "change in state");
assertNoNaN(stateNew, "new state post BC");
}
this.state = this._arrayClip2Zero(stateNew);
return stateNew;
}
_updateMeasurement(measurementType, value, position, context) {
switch(measurementType) {
case "oxygen":
grid_pos = Math.round(position * this.n_x);
this.state[grid_pos][S_O_INDEX] = value; // naive approach for reconciling measurements and simulation
break;
}
super._updateMeasurement(measurementType, value, position, context);
}
/**
* 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])));
} else {
state[0] = state[1];
}
// 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;
}
}
/**
* Create central finite difference second derivative operator.
* @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 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;
}
}
module.exports = { Reactor_CSTR, Reactor_PFR };
// DEBUG
// state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
// let initial_state = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1];
// 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 < 5000) {
// console.log(Reactor.tick(0.001));
// N += 1;
// }

18
src/utils.js Normal file
View File

@@ -0,0 +1,18 @@
/**
* Assert that no NaN values are present in an array.
* @param {Array} arr
* @param {string} label
*/
function assertNoNaN(arr, label = "array") {
if (Array.isArray(arr)) {
for (const el of arr) {
assertNoNaN(el, label);
}
} else {
if (Number.isNaN(arr)) {
throw new Error(`NaN detected in ${label}!`);
}
}
}
module.exports = { assertNoNaN };