From 2e76f733a8fed39a17f3f699a2c6b0506017a315 Mon Sep 17 00:00:00 2001 From: "p.vanderwilt" Date: Tue, 24 Jun 2025 16:23:33 +0200 Subject: [PATCH] Working on fixing the Derivative operators and BCs --- dependencies/reactor_class.js | 40 ++++++++++++++++++++--------------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js index e030637..9e02fca 100644 --- a/dependencies/reactor_class.js +++ b/dependencies/reactor_class.js @@ -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; // }