Update with ARRHENIUS version

wortiz · Jan 9, 2025 · ae1082e · ae1082e
1 parent debaea8
commit ae1082e
Show file tree

Hide file tree

Showing 5 changed files with 145 additions and 2 deletions.
diff --git a/include/mm_mp_const.h b/include/mm_mp_const.h
@@ -339,6 +339,7 @@ extern int Num_Var_Init_Mat[MAX_NUMBER_MATLS]; /* number of variables to overwri
 #define TURBULENT_K_OMEGA    54
 
 #define EPOXY_LINEAR_EXP 55
+#define EPOXY_ARRHENIUS_EXP 56
 /*
  *  Heat source modeling
  *

diff --git a/include/mm_std_models.h b/include/mm_std_models.h
@@ -68,6 +68,9 @@ EXTERN int epoxy_species_source /* mm_std_models.c                           */
 int epoxy_linear_exp_species_source(int species_no, /* Current species number */
                                     double *param); /* pointer to user-defined parameter list */
 
+int epoxy_arrhenius_exp_species_source(int species_no, /* Current species number */
+                                    double *param); /* pointer to user-defined parameter list */
+
 EXTERN int bond_species_source /* mm_std_models.c                           */
     (int,                      /* species_no - Current species number       */
      double *);                /* param - user-defined parameter list       */

diff --git a/src/mm_fill_species.c b/src/mm_fill_species.c
@@ -9917,6 +9917,26 @@ int get_continuous_species_terms(struct Species_Conservation_Terms *st,
         err = epoxy_linear_exp_species_source(w, mp->u_species_source[w]);
         st->MassSource[w] = mp->species_source[w];
 
+        if (af->Assemble_Jacobian) {
+          var = TEMPERATURE;
+          if (pd->v[pg->imtrx][var]) {
+            for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) {
+              st->d_MassSource_dT[w][j] = mp->d_species_source[var] * bf[var]->phi[j];
+            }
+          }
+
+          var = MASS_FRACTION;
+          if (pd->v[pg->imtrx][MASS_FRACTION]) {
+            var_offset = MAX_VARIABLE_TYPES + w;
+            for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) {
+              st->d_MassSource_dc[w][w][j] = mp->d_species_source[var_offset] * bf[var]->phi[j];
+            }
+          }
+        }
+      } else if (mp->SpeciesSourceModel[w] == EPOXY_ARRHENIUS_EXP) {
+        err = epoxy_arrhenius_exp_species_source(w, mp->u_species_source[w]);
+        st->MassSource[w] = mp->species_source[w];
+
         if (af->Assemble_Jacobian) {
           var = TEMPERATURE;
           if (pd->v[pg->imtrx][var]) {

diff --git a/src/mm_input_mp.c b/src/mm_input_mp.c
@@ -8563,7 +8563,34 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file)
       mat_ptr->u_species_source[species_no][6] = a6; /* A_n */
       mat_ptr->u_species_source[species_no][7] = a7; /* B_n */
 
-      SPF_DBL_VEC(endofstring(es), 6, mat_ptr->u_species_source[species_no]);
+      SPF_DBL_VEC(endofstring(es), 8, mat_ptr->u_species_source[species_no]);
+    } else if (!strcmp(model_name, "EPOXY_ARRHENIUS_EXP")) {
+      SpeciesSourceModel = EPOXY_ARRHENIUS_EXP;
+      model_read = 1;
+      mat_ptr->SpeciesSourceModel[species_no] = SpeciesSourceModel;
+      if (fscanf(imp, "%lf %lf %lf %lf %lf %lf %lf %lf", &a0, &a1, &a2, &a3, &a4, &a5, &a6, &a7) !=
+          8) {
+        sprintf(err_msg, "Matl %s needs 8 constants for %s %s model.\n", pd_glob[mn]->MaterialName,
+                "Species Source", "EPOXY_ARRHENIUS_EXP");
+        GOMA_EH(GOMA_ERROR, err_msg);
+      }
+
+      mat_ptr->u_species_source[species_no] = (dbl *)array_alloc(1, 8, sizeof(dbl));
+
+      mat_ptr->len_u_species_source[species_no] = 8;
+
+      mat_ptr->u_species_source[species_no][0] = a0; /* prefactor for k1 */
+      mat_ptr->u_species_source[species_no][1] = a1; /* exponent for k1 */
+      mat_ptr->u_species_source[species_no][2] = a2; /* prefactor for k2 */
+      mat_ptr->u_species_source[species_no][3] = a3; /* exponent for k2 */
+      // m = A_m * exp(B_m/T)
+      mat_ptr->u_species_source[species_no][4] = a4; /* A_m */
+      mat_ptr->u_species_source[species_no][5] = a5; /* B_m */
+      // m = A_n * exp(B_n/T)
+      mat_ptr->u_species_source[species_no][6] = a6; /* A_n */
+      mat_ptr->u_species_source[species_no][7] = a7; /* B_n */
+
+      SPF_DBL_VEC(endofstring(es), 8, mat_ptr->u_species_source[species_no]);
     } else if (!strcmp(model_name, "SSM_BOND")) {
       SpeciesSourceModel = SSM_BOND;
       model_read = 1;

diff --git a/src/mm_std_models.c b/src/mm_std_models.c
@@ -918,7 +918,99 @@ int epoxy_linear_exp_species_source(int species_no, /* Current species number */
         // k1 and k2 portion
         mp->d_species_source[var] = (k1 * E1 + k2 * E2 * alpha_m) * alpha_n / (T * T);
         // k2 * alpha_m portion
-        mp->d_species_source[var] += k2 * mt * log(alpha) * alpha_m * alpha_n;
+        if (alpha > 0) {
+          mp->d_species_source[var] += k2 * mt * log(alpha) * alpha_m * alpha_n;
+        }
+        //  alpha_n portion
+        mp->d_species_source[var] += nt * log(1.0 - alpha) * alpha_n * (k1 + k2 * alpha_m);
+      }
+    }
+  }
+  return 0;
+}
+
+int epoxy_arrhenius_exp_species_source(int species_no, /* Current species number */
+                                    double *param)  /* pointer to user-defined parameter list */
+
+{
+  /* Local Variables */
+  int eqn, var, var_offset, imtrx;
+  /*  int p, q, a, b, c;*/
+
+  /*  int v,w;*/
+
+  /*  int mdofs,vdofs;*/
+
+  /*  dbl C[MAX_CONC]; Convenient local variables */
+  dbl T; /* temperature for rate constants */
+  dbl A1, E1, A2, E2;
+  dbl A_m, B_m, A_n, B_n;
+  dbl k1, k2;
+  dbl alpha, alpha_m, alpha_m1, alpha_n, alpha_n1;
+
+  /* Begin Execution */
+
+  if (pd->gv[TEMPERATURE]) {
+    T = fv->T;
+  } else {
+    T = upd->Process_Temperature;
+  }
+
+  /* extent of reaction, alpha */
+  alpha = fv->c[species_no];
+  /*  if(alpha <= 0.) alpha = 0.0001; */
+  A1 = param[0];
+  E1 = param[1];
+  A2 = param[2];
+  E2 = param[3];
+  A_m = param[4];
+  B_m = param[5];
+  A_n = param[6];
+  B_n = param[7];
+
+  /* Arhenius type rate constants for extent of reaction model */
+  k1 = A1 * exp(-E1 / T);
+  k2 = A2 * exp(-E2 / T);
+
+  dbl m = A_m * exp(-B_m/T);
+  dbl mt = B_m * m / (T * T);
+  dbl n = A_n * exp(-B_n/T);
+  dbl nt = B_n * n / (T*T);
+
+  if (alpha > 0.0) {
+    alpha_m = pow(alpha, m);
+    alpha_m1 = pow(alpha, m - 1);
+  } else {
+    alpha_m = 0.;
+    alpha_m1 = 0.;
+  }
+
+  alpha_n = pow((1. - alpha), n);
+  alpha_n1 = pow((1. - alpha), n - 1);
+
+  /**********************************************************/
+
+  /* Species piece */
+  eqn = MASS_FRACTION;
+  for (imtrx = 0; imtrx < upd->Total_Num_Matrices; imtrx++) {
+    if (pd->e[imtrx][eqn] & T_SOURCE) {
+      mp->species_source[species_no] = (k1 + k2 * alpha_m) * alpha_n;
+
+      /* Jacobian entries for source term */
+      var = MASS_FRACTION;
+      if (pd->v[pg->imtrx][var]) {
+        var_offset = MAX_VARIABLE_TYPES + species_no;
+        mp->d_species_source[var_offset] =
+            (m * k2 * alpha_m1) * alpha_n - (k1 + k2 * alpha_m) * n * alpha_n1;
+      }
+      var = TEMPERATURE;
+      if (pd->v[pg->imtrx][var]) {
+        // k1 and k2 portion
+        mp->d_species_source[var] = (k1 * E1 + k2 * E2 * alpha_m) * alpha_n / (T * T);
+        // k2 * alpha_m portion
+        if (alpha > 0) {
+          mp->d_species_source[var] += k2 * mt * log(alpha) * alpha_m * alpha_n;
+        } 
         //  alpha_n portion
         mp->d_species_source[var] += nt * log(1.0 - alpha) * alpha_n * (k1 + k2 * alpha_m);
       }