From 76643b4286f21427f831c1064f25c97cf78872a9 Mon Sep 17 00:00:00 2001 From: Weston Ortiz Date: Wed, 13 Nov 2024 15:49:30 -0700 Subject: [PATCH] add EPOXY_LINEAR_EXP model --- include/mm_mp_const.h | 1 + include/mm_std_models.h | 3 + src/density.c | 14 +++- src/mm_fill_continuity.c | 161 +++++++++++++++++++++------------------ src/mm_fill_species.c | 20 +++++ src/mm_input_mp.c | 29 ++++++- src/mm_matrl.c | 13 +++- src/mm_std_models.c | 94 ++++++++++++++++++++++- 8 files changed, 252 insertions(+), 83 deletions(-) diff --git a/include/mm_mp_const.h b/include/mm_mp_const.h index 676929a69..6e385efda 100644 --- a/include/mm_mp_const.h +++ b/include/mm_mp_const.h @@ -338,6 +338,7 @@ extern int Num_Var_Init_Mat[MAX_NUMBER_MATLS]; /* number of variables to overwri #define TURBULENT_SA_DYNAMIC 53 /* Spallart Allmaras */ #define TURBULENT_K_OMEGA 54 +#define EPOXY_LINEAR_EXP 55 /* * Heat source modeling * diff --git a/include/mm_std_models.h b/include/mm_std_models.h index 865aeed45..91ed0b111 100644 --- a/include/mm_std_models.h +++ b/include/mm_std_models.h @@ -65,6 +65,9 @@ EXTERN int epoxy_species_source /* mm_std_models.c */ (int, /* species_no - Current species number */ double *); /* param - user-defined parameter list */ +int epoxy_linear_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/density.c b/src/density.c index 9d693fa5f..2f12c853f 100644 --- a/src/density.c +++ b/src/density.c @@ -807,14 +807,24 @@ double density(DENSITY_DEPENDENCE_STRUCT *d_rho, double time) dbl alpha_g = mp->u_density[3]; dbl cure_enable = mp->u_density[4]; - if (fv->c[0] >= cure_enable) { + int ConstitutiveEquation = gn_glob[ei[pg->imtrx]->mn]->ConstitutiveEquation; + int species_no = 0; + + if (ConstitutiveEquation == CURE || ConstitutiveEquation == EPOXY || + ConstitutiveEquation == FILLED_EPOXY || ConstitutiveEquation == FOAM_PMDI_10) { + species_no = gn_glob[ei[pg->imtrx]->mn]->cure_species_no; + } else { + GOMA_EH(GOMA_ERROR, "Unknown constitutive equation for density model CURE_SHRINKAGE"); + } + + if (fv->c[species_no] >= cure_enable) { rho = rho_l + ((rho_s - rho_l) / (alpha_m - alpha_g)) * (fv->c[0] - alpha_g); if (d_rho != NULL) { var = MASS_FRACTION; if (pd->v[pg->imtrx][var]) { for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { - d_rho->C[0][j] = ((rho_s - rho_l) / (alpha_m - alpha_g)) * bf[var]->phi[j]; + d_rho->C[species_no][j] = ((rho_s - rho_l) / (alpha_m - alpha_g)) * bf[var]->phi[j]; } } } diff --git a/src/mm_fill_continuity.c b/src/mm_fill_continuity.c index 0fa0be89a..438a353ec 100644 --- a/src/mm_fill_continuity.c +++ b/src/mm_fill_continuity.c @@ -2159,111 +2159,122 @@ double FoamVolumeSource(double time, } } } else if (mp->DensityModel == DENSITY_CURE_SHRINKAGE) { - if (fv->c[0] >= mp->u_density[4]) { + int ConstitutiveEquation = gn_glob[ei[pg->imtrx]->mn]->ConstitutiveEquation; + int species_no = 0; - DENSITY_DEPENDENCE_STRUCT d_rho_struct; - DENSITY_DEPENDENCE_STRUCT *d_rho = &d_rho_struct; + if (ConstitutiveEquation == CURE || ConstitutiveEquation == EPOXY || + ConstitutiveEquation == FILLED_EPOXY || ConstitutiveEquation == FOAM_PMDI_10) { + species_no = gn_glob[ei[pg->imtrx]->mn]->cure_species_no; + } else { + GOMA_EH(GOMA_ERROR, "Unknown constitutive equation for density model CURE_SHRINKAGE"); + } - rho = density(d_rho, time); + if (fv->c[species_no] >= mp->u_density[4]) { - dbl d_grad_rho_dC[MAX_CONC][MDE] = {{0.0}}; - dbl d_grad_rho_dX[DIM][DIM][MDE] = {{{0.0}}}; + DENSITY_DEPENDENCE_STRUCT d_rho_struct; + DENSITY_DEPENDENCE_STRUCT *d_rho = &d_rho_struct; - dbl rho_l = mp->u_density[0]; - dbl rho_s = mp->u_density[1]; - dbl alpha_m = mp->u_density[2]; - dbl alpha_g = mp->u_density[3]; + rho = density(d_rho, time); - double rho_dot = ((rho_s - rho_l) / (alpha_m - alpha_g)) * (fv_dot->c[0]); + dbl d_grad_rho_dC[MAX_CONC][MDE] = {{0.0}}; + dbl d_grad_rho_dX[DIM][DIM][MDE] = {{{0.0}}}; - double grad_rho[DIM] = {0.0}; + dbl rho_l = mp->u_density[0]; + dbl rho_s = mp->u_density[1]; + dbl alpha_m = mp->u_density[2]; + dbl alpha_g = mp->u_density[3]; - for (a = 0; a < dim; a++) { - grad_rho[a] = ((rho_s - rho_l) / (alpha_m - alpha_g)) * (fv->grad_c[0][a]); - if (af->Assemble_Jacobian) { - var = MASS_FRACTION; - for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { - d_grad_rho_dC[a][j] += - ((rho_s - rho_l) / (alpha_m - alpha_g)) * (bf[var]->grad_phi[j][a]); - } - for (int b = 0; b < dim; b++) { - var = MESH_DISPLACEMENT1 + b; + double rho_dot = ((rho_s - rho_l) / (alpha_m - alpha_g)) * (fv_dot->c[0]); + + double grad_rho[DIM] = {0.0}; + + for (a = 0; a < dim; a++) { + grad_rho[a] = ((rho_s - rho_l) / (alpha_m - alpha_g)) * (fv->grad_c[0][a]); + if (af->Assemble_Jacobian) { + var = MASS_FRACTION; for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { - d_grad_rho_dX[a][b][j] += - ((rho_s - rho_l) / (alpha_m - alpha_g)) * (fv->d_grad_c_dmesh[a][0][b][j]); + d_grad_rho_dC[a][j] += + ((rho_s - rho_l) / (alpha_m - alpha_g)) * (bf[var]->grad_phi[j][a]); + } + for (int b = 0; b < dim; b++) { + var = MESH_DISPLACEMENT1 + b; + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + d_grad_rho_dX[a][b][j] += + ((rho_s - rho_l) / (alpha_m - alpha_g)) * (fv->d_grad_c_dmesh[a][0][b][j]); + } } } } - } - double inv_rho = 1 / rho; + double inv_rho = 1 / rho; - dbl include_xdot = 1; + dbl include_xdot = 1; - source = rho_dot; - for (a = 0; a < dim; a++) { - source += (fv->v[a] - include_xdot * fv_dot->x[a]) * grad_rho[a]; - } - source *= inv_rho; - - if (af->Assemble_Jacobian) { - var = TEMPERATURE; - if (pd->v[pg->imtrx][var]) { - for (j = 0; pd->v[pg->imtrx][var] && j < ei[pg->imtrx]->dof[var]; j++) { - dFVS_dT[j] = 0; - } + source = rho_dot; + for (a = 0; a < dim; a++) { + source += (fv->v[a] - include_xdot * fv_dot->x[a]) * grad_rho[a]; } - } + source *= inv_rho; - var = VELOCITY1; - if (pd->v[pg->imtrx][var]) { - for (a = 0; a < dim; a++) { - var = VELOCITY1 + a; - for (j = 0; pd->v[pg->imtrx][var] && j < ei[pg->imtrx]->dof[var]; j++) { - dFVS_dv[a][j] = inv_rho * (bf[var]->phi[j] * grad_rho[a]); + if (af->Assemble_Jacobian) { + var = TEMPERATURE; + if (pd->v[pg->imtrx][var]) { + for (j = 0; pd->v[pg->imtrx][var] && j < ei[pg->imtrx]->dof[var]; j++) { + dFVS_dT[j] = 0; + } } } - } - var = MESH_DISPLACEMENT1; - if (pd->v[pg->imtrx][var]) { - for (a = 0; a < dim; a++) { - for (int b = 0; b < dim; b++) { - var = MESH_DISPLACEMENT1 + b; + var = VELOCITY1; + if (pd->v[pg->imtrx][var]) { + for (a = 0; a < dim; a++) { + var = VELOCITY1 + a; for (j = 0; pd->v[pg->imtrx][var] && j < ei[pg->imtrx]->dof[var]; j++) { - dFVS_dx[b][j] += - inv_rho * ((-(1.0 + 2.0 * tt) / dt) * include_xdot * bf[var]->phi[j] * grad_rho[a] + - (fv->v[a] - include_xdot * fv_dot->x[a]) * d_grad_rho_dX[a][b][j]); + dFVS_dv[a][j] = inv_rho * (bf[var]->phi[j] * grad_rho[a]); } } } - } - var = MASS_FRACTION; - if (pd->v[pg->imtrx][var]) { - for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { - dFVS_dC[0][j] = inv_rho * ((rho_s - rho_l) / (alpha_m - alpha_g)) * - ((1.0 + 2.0 * tt) / dt) * bf[var]->phi[j] + - (-d_rho->C[0][j] * inv_rho * inv_rho) * rho_dot; + var = MESH_DISPLACEMENT1; + if (pd->v[pg->imtrx][var]) { for (a = 0; a < dim; a++) { - dFVS_dC[0][j] += - inv_rho * ((fv->v[a] - include_xdot * fv_dot->x[a]) * d_grad_rho_dC[a][j]) + - (-d_rho->C[0][j] * inv_rho * inv_rho) * (fv->v[a] - include_xdot * fv_dot->x[a]) * - grad_rho[a]; + for (int b = 0; b < dim; b++) { + var = MESH_DISPLACEMENT1 + b; + for (j = 0; pd->v[pg->imtrx][var] && j < ei[pg->imtrx]->dof[var]; j++) { + dFVS_dx[b][j] += + inv_rho * + ((-(1.0 + 2.0 * tt) / dt) * include_xdot * bf[var]->phi[j] * grad_rho[a] + + (fv->v[a] - include_xdot * fv_dot->x[a]) * d_grad_rho_dX[a][b][j]); + } + } } } - } - var = FILL; - if (pd->v[pg->imtrx][var]) { - for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { - source_c = inv_rho * d_rho->F[j] * source; - dFVS_dF[j] = source_c; + var = MASS_FRACTION; + if (pd->v[pg->imtrx][var]) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + dFVS_dC[species_no][j] = inv_rho * ((rho_s - rho_l) / (alpha_m - alpha_g)) * + ((1.0 + 2.0 * tt) / dt) * bf[var]->phi[j] + + (-d_rho->C[0][j] * inv_rho * inv_rho) * rho_dot; + for (a = 0; a < dim; a++) { + dFVS_dC[species_no][j] += + inv_rho * ((fv->v[a] - include_xdot * fv_dot->x[a]) * d_grad_rho_dC[a][j]) + + (-d_rho->C[0][j] * inv_rho * inv_rho) * (fv->v[a] - include_xdot * fv_dot->x[a]) * + grad_rho[a]; + } + } } + var = FILL; + if (pd->v[pg->imtrx][var]) { + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + source_c = inv_rho * d_rho->F[j] * source; + + dFVS_dF[j] = source_c; + } + } + } else { + source = 0; } - } else { - source = 0; - } } if (ls != NULL && mp->mp2nd != NULL && (mp->DensityModel != LEVEL_SET) && diff --git a/src/mm_fill_species.c b/src/mm_fill_species.c index d77acaec1..adedb21d4 100644 --- a/src/mm_fill_species.c +++ b/src/mm_fill_species.c @@ -9897,6 +9897,26 @@ int get_continuous_species_terms(struct Species_Conservation_Terms *st, err = epoxy_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_LINEAR_EXP) { + 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]) { diff --git a/src/mm_input_mp.c b/src/mm_input_mp.c index 79e40569b..69cafe7a3 100644 --- a/src/mm_input_mp.c +++ b/src/mm_input_mp.c @@ -8523,9 +8523,34 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file) mat_ptr->u_species_source[species_no][4] = a4; /* prefactor for k2, mid T */ SPF_DBL_VEC(endofstring(es), 5, mat_ptr->u_species_source[species_no]); - } + } else if (!strcmp(model_name, "EPOXY_LINEAR_EXP")) { + SpeciesSourceModel = EPOXY_LINEAR_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_LINEAR_EXP"); + GOMA_EH(GOMA_ERROR, err_msg); + } + + mat_ptr->u_species_source[species_no] = (dbl *)array_alloc(1, 8, sizeof(dbl)); - else if (!strcmp(model_name, "SSM_BOND")) { + 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 * T + B_m + mat_ptr->u_species_source[species_no][4] = a4; /* A_m */ + mat_ptr->u_species_source[species_no][5] = a5; /* B_m */ + // n = A_n * T + B_n + 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]); + } else if (!strcmp(model_name, "SSM_BOND")) { SpeciesSourceModel = SSM_BOND; model_read = 1; mat_ptr->SpeciesSourceModel[species_no] = SpeciesSourceModel; diff --git a/src/mm_matrl.c b/src/mm_matrl.c index 7a3fdef72..005d9d000 100644 --- a/src/mm_matrl.c +++ b/src/mm_matrl.c @@ -553,8 +553,17 @@ calc_density(MATRL_PROP_STRUCT *matrl, int doJac, PROPERTYJAC_STRUCT *densityJac dbl alpha_m = mp->u_density[2]; dbl alpha_g = mp->u_density[3]; dbl cure_enable = mp->u_density[4]; + int ConstitutiveEquation = gn_glob[ei[pg->imtrx]->mn]->ConstitutiveEquation; + int species_no = 0; - if (fv->c[0] >= cure_enable) { + if (ConstitutiveEquation == CURE || ConstitutiveEquation == EPOXY || + ConstitutiveEquation == FILLED_EPOXY || ConstitutiveEquation == FOAM_PMDI_10) { + species_no = gn_glob[ei[pg->imtrx]->mn]->cure_species_no; + } else { + GOMA_EH(GOMA_ERROR, "Unknown constitutive equation for density model CURE_SHRINKAGE"); + } + + if (fv->c[species_no] >= cure_enable) { rho = rho_l + ((rho_s - rho_l) / (alpha_m - alpha_g)) * (fv->c[0] - alpha_g); /* Now do sensitivies */ @@ -563,7 +572,7 @@ calc_density(MATRL_PROP_STRUCT *matrl, int doJac, PROPERTYJAC_STRUCT *densityJac if (doJac) { if (pd->v[pg->imtrx][var]) { double drhodC = ((rho_s - rho_l) / (alpha_m - alpha_g)); - propertyJac_addEnd(densityJac, MASS_FRACTION, matID, 0, drhodC, rho); + propertyJac_addEnd(densityJac, MASS_FRACTION, matID, species_no, drhodC, rho); } } } else { diff --git a/src/mm_std_models.c b/src/mm_std_models.c index a220539b3..3f6ae7a80 100644 --- a/src/mm_std_models.c +++ b/src/mm_std_models.c @@ -838,6 +838,94 @@ int epoxy_species_source(int species_no, /* Current species number */ } return 0; } + +int epoxy_linear_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 * T + B_m; + dbl mt = A_m; + dbl n = A_n * T + B_n; + dbl nt = A_n; + + 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 + 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; +} /*****************************************************************************/ /* END of routine epoxy_species_source */ /*****************************************************************************/ @@ -1541,7 +1629,8 @@ double epoxy_heat_source(HEAT_SOURCE_DEPENDENCE_STRUCT *d_h, /* find equation that has extent of reaction */ for (w = 0; w < pd->Num_Species_Eqn; w++) { - if (mp->SpeciesSourceModel[w] == EPOXY || mp->SpeciesSourceModel[w] == EPOXY_DEA) { + if (mp->SpeciesSourceModel[w] == EPOXY || mp->SpeciesSourceModel[w] == EPOXY_DEA || + mp->SpeciesSourceModel[w] == EPOXY_LINEAR_EXP) { /* extent of reaction, alpha */ species_no = w; } @@ -1563,7 +1652,8 @@ double epoxy_heat_source(HEAT_SOURCE_DEPENDENCE_STRUCT *d_h, var = MASS_FRACTION; if (d_h != NULL && pd->v[pg->imtrx][var]) { if (mp->SpeciesSourceModel[species_no] == EPOXY || - (mp->SpeciesSourceModel[species_no] == EPOXY_DEA)) { + (mp->SpeciesSourceModel[species_no] == EPOXY_DEA) || + (mp->SpeciesSourceModel[species_no] == EPOXY_LINEAR_EXP)) { for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { d_h->C[species_no][j] += delta_h * (1 + 2. * tt) / dt * bf[var]->phi[j]; }