Skip to content

Commit

Permalink
add EPOXY_LINEAR_EXP model
Browse files Browse the repository at this point in the history
  • Loading branch information
@wortiz
wortiz committed Nov 13, 2024
1 parent dab0eb4 commit 76643b4
Show file tree
Hide file tree
Showing 8 changed files with 252 additions and 83 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 @@ -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
*
Expand Down
3 changes: 3 additions & 0 deletions include/mm_std_models.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
Expand Down
14 changes: 12 additions & 2 deletions src/density.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
}
}
}
Expand Down
161 changes: 86 additions & 75 deletions src/mm_fill_continuity.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) &&
Expand Down
20 changes: 20 additions & 0 deletions src/mm_fill_species.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]) {
Expand Down
29 changes: 27 additions & 2 deletions src/mm_input_mp.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
13 changes: 11 additions & 2 deletions src/mm_matrl.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
Expand All @@ -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 {
Expand Down
Loading

0 comments on commit 76643b4

Please sign in to comment.