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

Working on fixing the Derivative operators and BCs

Browse Source
This commit is contained in:
p.vanderwilt
2025-06-24 16:23:33 +02:00
parent f2d94b26c5
commit 2e76f733a8
1 changed files with 23 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

40
dependencies/reactor_class.js vendored
View File

@@ -69,7 +69,8 @@ class Reactor_CSTR {
const dC_total = math.multiply(math.add(dC_in, dC_out, r, t_O), time_step);
this.state = math.abs(math.add(this.state, dC_total)); // make sure that concentrations do not go negative
// clip value element-wise to each subarray to avoid negative concentrations
this.state = math.add(this.state, dC_total).map(val => val < 0 ? 0 : val);
return this.state;
}
}
@@ -99,8 +100,8 @@ class Reactor_PFR {
this.kla = kla; // if NaN, use external OTR [d-1]
this.currentTime = Date.now(); // milliseconds since epoch [ms]
this.timeStep = 1/(24*60*60); // time step [d]
this.speedUpFactor = 1;
this.timeStep = 1/(24*60*15); // time step [d]
this.speedUpFactor = 60;
this.D_op = this.makeDoperator();
this.D2_op = this.makeD2operator();
@@ -147,7 +148,7 @@ class Reactor_PFR {
tick_fe(time_step) { // tick reactor state using forward Euler method
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(-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));
reaction[0] = Array(13).fill(0.0);
const transfer = Array.from(Array(this.n_x), () => new Array(13).fill(0.0));
@@ -158,36 +159,40 @@ class Reactor_PFR {
transfer.forEach((x, i) => { x[0] = this.calcOTR(this.state[i][0]); });
}
transfer[0][0] = 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_C_in = math.multiply(1/math.sum(this.Fs),[this.Fs], this.Cs_in)[0];
const BC_gradient = Array(this.n_x).fill(0.0);
BC_gradient[0] = 1;
BC_gradient[1] = -1;
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)[0];
this.state[0] = math.add(BC_influx, BC_dispersion);
this.state[0] = math.add(BC_C_in, BC_dispersion);
console.log(BC_dispersion);
}
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
// clip value element-wise to each subarray to avoid negative concentrations
this.state = math.add(this.state, dC_total).map(row => row.map(val => val < 0 ? 0 : val));
return this.state;
}
makeDoperator() { // create the upwind scheme gradient operator
const I = math.identity(this.n_x);
const A = math.resize(math.diag(Array(this.n_x).fill(-1), 1), [this.n_x, this.n_x]);
const I = math.diag(Array(this.n_x).fill(-1), 0);
const A = math.resize(math.diag(Array(this.n_x).fill(1), 1), [this.n_x, this.n_x]);
I[0][0] = 0;
I[0][1] = 1;
I[0][1] = -1;
I[this.n_x-1][this.n_x-1] = 0; // Neumann boundary condition at x=L
return math.add(I, A);
}
makeD2operator() { // create the upwind scheme second derivative operator
const I = math.identity(this.n_x);
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 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]);
I[0][0] = 0;
I[0][1] = 1;
I[0][1] = -1;
return math.add(I, A, B);
}
}
@@ -199,8 +204,9 @@ class Reactor_PFR {
// 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 < 500) {
// while (N < 5000) {
// console.log(Reactor.tick_fe(0.001));
// N += 1;
// }