From e6c1e21c166ac90751e8e320344bf741b9b4bfc8 Mon Sep 17 00:00:00 2001
From: "p.vanderwilt" <p.van.der.wilt@brabantsedelta.nl>
Date: Mon, 23 Jun 2025 17:46:55 +0200
Subject: [PATCH] Implement Danckwerts boundary condition in tick_fe method for
 Reactor_PFR

---
 dependencies/reactor_class.js | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/dependencies/reactor_class.js b/dependencies/reactor_class.js
index aa7f9dd..bb5a96c 100644
--- a/dependencies/reactor_class.js
+++ b/dependencies/reactor_class.js
@@ -137,14 +137,11 @@ class Reactor_PFR {
     }
 
     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 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 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
             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]); });
         }
 
+        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);
 
         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
         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]);
+        I[0, 0] = 0;
+        I[0, 1] = 0;
         I[this.n_x-1, this.n_x-1] = 0; // Neumann boundary condition at x=L
         return math.add(I, A);
     }
@@ -169,7 +177,8 @@ class Reactor_PFR {
         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 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);
     }
 }