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

Implement Danckwerts boundary condition in tick_fe method for Reactor_PFR

Browse Source
This commit is contained in:
p.vanderwilt
2025-06-23 17:46:55 +02:00
parent 70531a3a59
commit e6c1e21c16
1 changed files with 15 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

21
dependencies/reactor_class.js vendored
View File

@@ -137,14 +137,11 @@ class Reactor_PFR {
} }
tick_fe(time_step) { // tick reactor state using forward Euler method tick_fe(time_step) { // tick reactor state using forward Euler method
if (math.sum(this.Fs) > 0) {
this.state[0] = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0] // Dichelet boundary condition
}
const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state); const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.state);
const advection = math.multiply(math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state); const advection = math.multiply(math.sum(this.Fs)/(this.A*this.d_x), this.D_op, this.state);
const reaction = this.state.map(this.asm.compute_dC); const reaction = this.state.map(this.asm.compute_dC);
const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0)) reaction[0] = Array(13).fill(0.0);
const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
transfer.forEach((x) => { x[0] = this.OTR; }); transfer.forEach((x) => { x[0] = this.OTR; });
@@ -152,6 +149,15 @@ class Reactor_PFR {
transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); }); transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); });
} }
if (math.sum(this.Fs) > 0) { // Danckwerts BC
const BC_influx = math.multiply(math.divide([this.Fs], this.A), this.Cs_in)[0];
const BC_gradient = Array(this.n_x).fill(0.0);
BC_gradient[0] = 1;
BC_gradient[1] = -1;
const BC_dispersion = math.multiply(this.D * this.A / (math.sum(this.Fs)*this.d_x), [BC_gradient], this.state);
this.state[0] = math.add(BC_influx, BC_dispersion);
}
const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step); const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer), time_step);
this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
@@ -161,6 +167,8 @@ class Reactor_PFR {
makeDoperator() { // create the upwind scheme gradient operator makeDoperator() { // create the upwind scheme gradient operator
const I = math.identity(this.n_x); const I = math.identity(this.n_x);
const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]); const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
I[0, 0] = 0;
I[0, 1] = 0;
I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L
return math.add(I, A); return math.add(I, A);
} }
@@ -169,7 +177,8 @@ class Reactor_PFR {
const I = math.diag(Array(this.n_x).fill(2), 0); const I = math.diag(Array(this.n_x).fill(2), 0);
const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]); const A = math.diag(Array(this.n_x).fill(-1), 1).resize([this.n_x, this.n_x]);
const B = math.diag(Array(this.n_x).fill(-1), -1).resize([this.n_x, this.n_x]); const B = math.diag(Array(this.n_x).fill(-1), -1).resize([this.n_x, this.n_x]);
I[0, 0] = 1; I[0, 0] = 0;
I[0, 1] = 0;
return math.add(I, A, B); return math.add(I, A, B);
} }
} }