Skip to content

Commit

Permalink
Add Density CURE_SHRINKAGE model
Browse files Browse the repository at this point in the history
  • Loading branch information
@wortiz
wortiz committed Oct 17, 2024
1 parent 011cd3d commit a401498
Show file tree
Hide file tree
Showing 5 changed files with 152 additions and 3 deletions.
1 change: 1 addition & 0 deletions include/mm_mp_const.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
/**********************************************************************************/

/*
Expand Down
17 changes: 17 additions & 0 deletions src/density.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
}
Expand Down
110 changes: 107 additions & 3 deletions src/mm_fill_continuity.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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];
}
Expand Down Expand Up @@ -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) &&
Expand Down
10 changes: 10 additions & 0 deletions src/mm_input_mp.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
17 changes: 17 additions & 0 deletions src/mm_matrl.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit a401498

Please sign in to comment.