Merge pull request 'Added temperature-dependent rate calculation to ASM3 model' (#12) from experimental into main

Reviewed-on: p.vanderwilt/asm3#12

advanced-reactor.html

@@ -8,6 +8,7 @@
             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 },
@@ -75,6 +76,10 @@
                     $(".PFR").show();
                 }
             });
+            $("#node-input-alpha").typedInput({
+                type:"num",
+                types:["num"]
+            })
             // Set initial visibility on dialog open
             const initialType = $("#node-input-reactor_type").typedInput("value");
             if (initialType === "CSTR") {
@@ -118,6 +123,13 @@
         <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="#">

src/nodeClass.js

@@ -70,6 +70,7 @@ class nodeClass {
             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: [

src/reaction_modules/asm3_class.js

@@ -69,6 +69,28 @@ class ASM3 {
       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();
   }
 
@@ -118,45 +140,71 @@ class ASM3 {
   }
 
   /**
-   * Computes the reaction rates for each process reaction based on the current state.
+   * 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) {
+  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 : k_H * this._monod(X_S / X_H, K_X) * X_H;
+    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] = k_STO * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
+    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] = k_STO * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * 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 : mu_H_max * 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[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 : mu_H_max * 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[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] = b_H_O * this._monod(S_O, K_O) * X_H;
+    rates[5] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
-    rates[6] = b_H_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * 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] = b_STO_O * this._monod(S_O, K_O) * X_H;
+    rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
-    rates[8] = b_STO_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
+    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] = mu_A_max * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
+    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] = b_A_O * this._monod(S_O, K_O) * X_A;
+    rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
-    rates[11] = b_A_NO * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * 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.
+   * 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) { // compute changes in 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));
+    return math.multiply(this.stoi_matrix, this.compute_rates(state, T));
   }
 }
 

src/specificClass.js

@@ -10,7 +10,7 @@ const math = create(all, config);
 
 const S_O_INDEX = 0;
 const NUM_SPECIES = 13;
-const DEBUG = false;
+const DEBUG = true;
 
 class Reactor {
     /**
@@ -30,7 +30,7 @@ class Reactor {
 
         this.currentTime = Date.now(); // milliseconds since epoch [ms]
         this.timeStep = 1 / (24*60*15); // time step [d]
-        this.speedUpFactor = 60; // speed up factor for simulation, 60 means 1 minute per simulated second
+        this.speedUpFactor = 1; // speed up factor for simulation, 60 means 1 minute per simulated second
     }
 
     /**
@@ -149,6 +149,8 @@ class Reactor_PFR extends Reactor {
         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())
 
         // console.log("Initial State: ")
@@ -178,6 +180,7 @@ class Reactor_PFR extends Reactor {
     set setInfluent(input) {
         super.setInfluent = input;
         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 " + this.d_x*math.sum(this.Fs)/(this.D*this.A));
             console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
@@ -205,9 +208,7 @@ class Reactor_PFR extends Reactor {
             const BC_gradient = Array(this.n_x).fill(0);
             BC_gradient[0] = -1;
             BC_gradient[1] = 1;
-            let Pe = this.length * math.sum(this.Fs) / (this.D * this.A);
+            const BC_dispersion = math.multiply((1 - this.alpha) * this.D*this.A / (math.sum(this.Fs) * this.d_x), [BC_gradient], state)[0];
-            let residence_time = this.volume/math.sum(this.Fs);
-            const BC_dispersion = math.multiply((1 - (1 + 4*residence_time/Pe)^0.5) / (Pe*this.d_x), [BC_gradient], state)[0];
             state[0] = math.add(BC_C_in, BC_dispersion).map(val => val < 0 ? 0 : val);
         } else { // Neumann BC (no flux)
             state[0] = state[1];