From a40149841fe2971eb8d5bd546c832df4c1a22fbd Mon Sep 17 00:00:00 2001 From: Weston Ortiz Date: Thu, 17 Oct 2024 15:56:18 -0600 Subject: [PATCH] Add Density CURE_SHRINKAGE model --- include/mm_mp_const.h | 1 + src/density.c | 17 ++++++ src/mm_fill_continuity.c | 110 +++++++++++++++++++++++++++++++++++++-- src/mm_input_mp.c | 10 ++++ src/mm_matrl.c | 17 ++++++ 5 files changed, 152 insertions(+), 3 deletions(-) diff --git a/include/mm_mp_const.h b/include/mm_mp_const.h index 473fca262..676929a69 100644 --- a/include/mm_mp_const.h +++ b/include/mm_mp_const.h @@ -229,6 +229,7 @@ extern int Num_Var_Init_Mat[MAX_NUMBER_MATLS]; /* number of variables to overwri #define DENSITY_FOAM_PBE_EQN 35 #define DENSITY_FOAM_PMDI_10 20 #define DENSITY_MOMENT_BASED 21 +#define DENSITY_CURE_SHRINKAGE 906 /**********************************************************************************/ /* diff --git a/src/density.c b/src/density.c index d73fc5ae5..5c8545ea5 100644 --- a/src/density.c +++ b/src/density.c @@ -799,6 +799,23 @@ double density(DENSITY_DEPENDENCE_STRUCT *d_rho, double time) } } } + } else if (mp->DensityModel == DENSITY_CURE_SHRINKAGE) { + // parameters + 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 = 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]; + } + } + } } else { GOMA_EH(GOMA_ERROR, "Unrecognized density model"); } diff --git a/src/mm_fill_continuity.c b/src/mm_fill_continuity.c index 89280eb8f..ace6e064a 100644 --- a/src/mm_fill_continuity.c +++ b/src/mm_fill_continuity.c @@ -490,7 +490,8 @@ int assemble_continuity(dbl time_value, /* current time */ if (mp->DensityModel == DENSITY_FOAM || mp->DensityModel == DENSITY_FOAM_CONC || mp->DensityModel == DENSITY_FOAM_TIME || mp->DensityModel == DENSITY_FOAM_TIME_TEMP || mp->DensityModel == DENSITY_MOMENT_BASED || - mp->DensityModel == DENSITY_FOAM_PMDI_10) { + mp->DensityModel == DENSITY_FOAM_PMDI_10 || + mp->DensityModel == DENSITY_CURE_SHRINKAGE) { /* These density models locally permit a time and spatially varying density. Consequently, the Lagrangian derivative of the density terms in the continuity equation are not zero and are @@ -1085,7 +1086,8 @@ int assemble_continuity(dbl time_value, /* current time */ if (mp->DensityModel == DENSITY_FOAM || mp->DensityModel == DENSITY_FOAM_CONC || mp->DensityModel == DENSITY_FOAM_TIME || mp->DensityModel == DENSITY_FOAM_PBE || mp->DensityModel == DENSITY_FOAM_TIME_TEMP || - mp->DensityModel == DENSITY_FOAM_PMDI_10) { + mp->DensityModel == DENSITY_FOAM_PMDI_10 || + mp->DensityModel == DENSITY_CURE_SHRINKAGE) { source = sourceBase * d_h3detJ_dmesh_bj * wt + dFVS_dx[b][j] * d_area; source *= phi_i * source_etm; } else if (mp->DensityModel == REACTIVE_FOAM) { @@ -1193,7 +1195,8 @@ int assemble_continuity(dbl time_value, /* current time */ /* Foaming volume source term */ if (mp->DensityModel == REACTIVE_FOAM || mp->DensityModel == DENSITY_FOAM || - mp->DensityModel == DENSITY_FOAM_PBE) { + mp->DensityModel == DENSITY_FOAM_PBE || + mp->DensityModel == DENSITY_CURE_SHRINKAGE) { source += dFVS_dC[w][j]; source *= phi_i * h3 * det_J * wt * pd->etm[pg->imtrx][eqn][LOG2_SOURCE]; } @@ -2155,6 +2158,107 @@ double FoamVolumeSource(double time, } } } + } else if (mp->DensityModel == DENSITY_CURE_SHRINKAGE) { + DENSITY_DEPENDENCE_STRUCT d_rho_struct; + DENSITY_DEPENDENCE_STRUCT *d_rho = &d_rho_struct; + + rho = density(d_rho, time); + + dbl d_grad_rho_dC[MAX_CONC][MDE] = {{0.0}}; + dbl d_grad_rho_dX[DIM][DIM][MDE] = {{{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]; + + 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_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; + + 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; + } + } + } + + 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]); + } + } + } + + 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; + 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 = 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; + 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]; + } + } + } + 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; + } + } } if (ls != NULL && mp->mp2nd != NULL && (mp->DensityModel != LEVEL_SET) && diff --git a/src/mm_input_mp.c b/src/mm_input_mp.c index e2d6e4e6b..f841924f0 100644 --- a/src/mm_input_mp.c +++ b/src/mm_input_mp.c @@ -632,6 +632,16 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file) } mat_ptr->len_u_density = num_const; SPF_DBL_VEC(endofstring(es), num_const, mat_ptr->u_density); + } else if (model_read == -1 && !strcmp(model_name, "CURE_SHRINKAGE")) { + mat_ptr->DensityModel = DENSITY_CURE_SHRINKAGE; + num_const = read_constants(imp, &(mat_ptr->u_density), 0); + if (num_const != 4) { + sprintf(err_msg, "Material %s - expected 4 constants for %s %s model.\n", + pd_glob[mn]->MaterialName, "Density", "CURE_SHRINKAGE"); + GOMA_EH(GOMA_ERROR, err_msg); + } + mat_ptr->len_u_density = num_const; + SPF_DBL_VEC(endofstring(es), num_const, mat_ptr->u_density); } else { sprintf(err_msg, "Material %s - unrecognized model for %s \"%s\" ???\n", pd_glob[mn]->MaterialName, "Density", model_name); diff --git a/src/mm_matrl.c b/src/mm_matrl.c index 26ca9ce88..cb5a331df 100644 --- a/src/mm_matrl.c +++ b/src/mm_matrl.c @@ -547,6 +547,23 @@ calc_density(MATRL_PROP_STRUCT *matrl, int doJac, PROPERTYJAC_STRUCT *densityJac propertyJac_addEnd(densityJac, TEMPERATURE, matID, 0, drhoDT, rho); } } + } else if (matrl->DensityModel == DENSITY_CURE_SHRINKAGE) { + 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 = rho_l + ((rho_s - rho_l) / (alpha_m - alpha_g)) * (fv->c[0] - alpha_g); + + /* Now do sensitivies */ + + int var = MASS_FRACTION; + 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); + } + } } else if (matrl->DensityModel == DENSITY_MOMENT_BASED) { int var;