Skip to content

Commit

Permalink
Lubrication Logarithm Interpolation (goma#468)
Browse files Browse the repository at this point in the history 

Added Logarithmic interpolation across the level-set interface for shell lubrication equations. This seems to make the viscosity behavior for shear-thinning liquids much smoother across the LS interface.

It's invoked via
Lubrication Level-Set Interpolation = LOGARITHMIC
i.e., interpolated value pp = p0^1-H(F) * p1^H(F)

The default is the standard linear interpolation, pp = p0 * (1-H(F)) + p1 * H(F)

At present, only the lubrication flow q(gradP) can be interpolated this way.
  • Loading branch information
@rbsecor
rbsecor authored Jul 24, 2024
1 parent 5230681 commit 86061fd
Show file tree
Hide file tree
Showing 6 changed files with 163 additions and 48 deletions.
1 change: 1 addition & 0 deletions include/mm_mp_structs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -773,6 +773,7 @@ struct Material_Properties {
double Lub_Kwt_func;
int Lub_Curv_MassLump;
int Lub_Curv_Modulation;
int Lub_LS_Interpolation;

int Lub_Heat_XferModel;
double Lub_Heat_Xfer;
Expand Down
1 change: 1 addition & 0 deletions include/rf_fem_const.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,7 @@
#define LUB_VISCINT_ANALYTICAL 26
#define LUB_VISCINT_POWERLAW 27
#define MAX_LUB_NGP 5
#define LOGARITHMIC 29

/* Residence time kernel functions */
#define LINEAR_TIMETEMP 1110
Expand Down
1 change: 1 addition & 0 deletions src/dp_vif.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1581,6 +1581,7 @@ void noahs_ark(void) {
ddd_add_member(n, &mp_glob[i]->Lub_Curv_Modulation, 1, MPI_INT);
ddd_add_member(n, &mp_glob[i]->Lub_Heat_XferModel, 1, MPI_INT);
ddd_add_member(n, &mp_glob[i]->Lub_Heat_TambModel, 1, MPI_INT);
ddd_add_member(n, &mp_glob[i]->Lub_LS_Interpolation, 1, MPI_INT);
ddd_add_member(n, &mp_glob[i]->PorousShellClosedPorosityModel, 1, MPI_INT);
ddd_add_member(n, &mp_glob[i]->PorousShellClosedRadiusModel, 1, MPI_INT);
ddd_add_member(n, &mp_glob[i]->PorousShellClosedHeightModel, 1, MPI_INT);
Expand Down
1 change: 0 additions & 1 deletion src/mm_fill_species.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8852,7 +8852,6 @@ int get_convection_velocity(
volso -= fv->c[w] * mp->specific_volume[w];
}
GOMA_WH(-1, "negative nonvolatile volume fraction %g %g", volsolid, fv->c[0]);
fprintf(stderr, "density %g %g\n", mp->density, mp->specific_volume[0]);
}

for (p = 0; p < VIM; p++) {
Expand Down
20 changes: 19 additions & 1 deletion src/mm_input_mp.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10385,6 +10385,7 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file)
}
ECHO(es, echo_file);
}

strcpy(search_string, "Lubrication Curvature Normal");
model_read = look_for_mat_prop(imp, search_string, NULL, NULL, NO_USER, NULL, model_name,
SCALAR_INPUT, &NO_SPECIES, es);
Expand All @@ -10400,7 +10401,24 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file)
SPF(es, "\t(%s = %s)", search_string, "GRADH");
}
ECHO(es, echo_file);
}

strcpy(search_string, "Lubrication Level-Set Interpolation");
model_read = look_for_mat_prop(imp, search_string, NULL, NULL, NO_USER, NULL, model_name,
SCALAR_INPUT, &NO_SPECIES, es);
if (!strcasecmp(model_name, "LINEAR") || !strcasecmp(model_name, "linear")) {
mat_ptr->Lub_LS_Interpolation = LINEAR;
SPF(es, "\t(%s = %s)", search_string, "LINEAR");
} else if (!strcasecmp(model_name, "LOGARITHMIC") || !strcasecmp(model_name, "LOG")) {
mat_ptr->Lub_LS_Interpolation = LOGARITHMIC;
SPF(es, "\t(%s = %s)", search_string, "LOGARITHMIC");
} else {
mat_ptr->Lub_LS_Interpolation = LINEAR;
GOMA_WH(model_read, "Defaulting on Lubrication LS Interpolation");
SPF(es, "\t(%s = %s)", search_string, "LINEAR");
}
ECHO(es, echo_file);

} /* End of Lubrication Main Section */

/* Shell Lubrication Curvature Diffusion Term */
if (pd_glob[mn]->gv[R_SHELL_LUB_CURV] || pd_glob[mn]->gv[R_SHELL_LUB_CURV_2]) {
Expand Down
187 changes: 141 additions & 46 deletions src/mm_shell_util.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4168,6 +4168,7 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co
so need to exclude those */
if (pd->v[pg->imtrx][VAR] && (nonmoving_model || gn->ConstitutiveEquation == POWER_LAW ||
gn->ConstitutiveEquation == HERSCHEL_BULKLEY)) {
load_lsi(ls->Length_Scale);
if (mp->mp2nd->ViscosityModel == RATIO) {
ratio = 1. / mp->mp2nd->viscosity; /* Assuming model = RATIO for now */
q_mag2 = q_mag * ratio;
Expand All @@ -4183,20 +4184,55 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co
vis_w /= factor;
} else if (nonmoving_model && (mp->mp2nd->ViscosityModel == CONSTANT ||
mp->mp2nd->ViscosityModel == TIME_RAMP)) {
double dq_gradp2, pre_delP2;
k_turb = 12.;
dq_gradp2 = pre_delP2 = -CUBE(H) / (k_turb * mp->mp2nd->viscosity);
q_mag2 = pre_delP2 * pgrad;
q_mag = ls_modulate_property(q_mag, q_mag2, ls->Length_Scale,
(double)mp->mp2nd->viscositymask[0],
(double)mp->mp2nd->viscositymask[1], dqmag_dF, &factor);
dq_gradp = dq_gradp * factor + dq_gradp2 * (1. - factor);
pre_delP = pre_delP * factor + pre_delP2 * (1. - factor);
dq_dH = dq_dH * factor +
(1. - factor) * (-3. * SQUARE(H) / (k_turb * mp->mp2nd->viscosity) * pgrad);
dq_dT *= factor; // mp2nd->viscosity is independent of Temperature
srate = srate * factor + (1. - factor) * fabs(tau_w / mp->mp2nd->viscosity);
vis_w = vis_w * factor + (1. - factor) * mp->mp2nd->viscosity;
if (mp->Lub_LS_Interpolation == LOGARITHMIC) {
if (lsi->near || (fv->F > 0 && mp->mp2nd->viscositymask[1]) ||
(fv->F < 0 && mp->mp2nd->viscositymask[0])) {
double dq_gradp2, pre_delP2, dq_dH2, srate2, qmag_log;
k_turb = 12.;
dq_gradp2 = pre_delP2 = -CUBE(H) / (k_turb * mp->mp2nd->viscosity);
q_mag2 = pre_delP2 * pgrad;
dq_dH2 = -3. * SQUARE(H) / (k_turb * mp->mp2nd->viscosity) * pgrad;
srate2 = tau_w / mp->mp2nd->viscosity;
qmag_log = (DOUBLE_NONZERO(q_mag) ? log(q_mag2 / q_mag) : 0.0);
if (fabs(fv->F) > ls->Length_Scale) {
q_mag = q_mag2;
dq_gradp = dq_gradp2;
pre_delP = pre_delP2;
dq_dH = dq_dH2;
srate = srate2;
vis_w = mp->mp2nd->viscosity;
} else {
factor = (mp->mp2nd->viscositymask[1] ? (1.0 - lsi->H) : lsi->H);
q_mag = -pow(-q_mag, factor) * pow(-q_mag2, 1.0 - factor);
dq_gradp = -pow(-dq_gradp, factor) * pow(-dq_gradp2, 1.0 - factor);
pre_delP = -pow(-pre_delP, factor) * pow(-pre_delP2, 1.0 - factor);
dq_dH = -pow(-dq_dH, factor) * pow(-dq_dH2, 1.0 - factor);
srate = pow(srate, factor) * pow(srate2, 1.0 - factor);
vis_w = pow(vis_w, factor) * pow(mp->mp2nd->viscosity, 1.0 - factor);
for (j = 0; j < ei[pg->imtrx]->dof[VAR]; j++) {
dqmag_dF[j] += q_mag * qmag_log * lsi->d_H_dF[j];
}
}
}
} else if (mp->Lub_LS_Interpolation == LINEAR) {
double dq_gradp2, pre_delP2, dq_dH2, srate2;
k_turb = 12.;
dq_gradp2 = pre_delP2 = -CUBE(H) / (k_turb * mp->mp2nd->viscosity);
q_mag2 = pre_delP2 * pgrad;
dq_dH2 = -3. * SQUARE(H) / (k_turb * mp->mp2nd->viscosity) * pgrad;
srate2 = tau_w / mp->mp2nd->viscosity;
q_mag = ls_modulate_property(q_mag, q_mag2, ls->Length_Scale,
(double)mp->mp2nd->viscositymask[0],
(double)mp->mp2nd->viscositymask[1], dqmag_dF, &factor);
dq_gradp = dq_gradp * factor + dq_gradp2 * (1. - factor);
pre_delP = pre_delP * factor + pre_delP2 * (1. - factor);
dq_dH = dq_dH * factor + (1. - factor) * dq_dH2;
dq_dT *= factor; // mp2nd->viscosity is independent of Temperature
srate = srate * factor + (1. - factor) * srate2;
vis_w = vis_w * factor + (1. - factor) * mp->mp2nd->viscosity;
} else {
GOMA_WH(GOMA_ERROR, "mp->Lub_LS_Interpolation needs to be LOG or LINEAR...\n");
}
} else {
GOMA_WH(GOMA_ERROR, "mp2nd->ViscosityModel needs to be RATIO or CONSTANT...\n");
}
Expand Down Expand Up @@ -4262,20 +4298,55 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co
vis_w /= factor;
} else if (mp->mp2nd->ViscosityModel == CONSTANT ||
mp->mp2nd->ViscosityModel == TIME_RAMP) {
double dq_gradp2, pre_delP2;
k_turb = 12.;
dq_gradp2 = pre_delP2 = -CUBE(H) / (k_turb * mp->mp2nd->viscosity);
q_mag2 = pre_delP2 * pgrad;
q_mag = ls_modulate_property(q_mag, q_mag2, ls->Length_Scale,
(double)mp->mp2nd->viscositymask[0],
(double)mp->mp2nd->viscositymask[1], dqmag_dF, &factor);
dq_gradp = dq_gradp * factor + dq_gradp2 * (1. - factor);
pre_delP = pre_delP * factor + pre_delP2 * (1. - factor);
dq_dH = dq_dH * factor +
(1. - factor) * (-3. * SQUARE(H) / (k_turb * mp->mp2nd->viscosity) * pgrad);
dq_dT *= factor; // mp2nd->viscosity is independent of Temperature
srate = srate * factor + (1. - factor) * fabs(tau_w / mp->mp2nd->viscosity);
vis_w = vis_w * factor + (1. - factor) * mp->mp2nd->viscosity;
if (mp->Lub_LS_Interpolation == LOGARITHMIC) {
if (lsi->near || (fv->F > 0 && mp->mp2nd->viscositymask[1]) ||
(fv->F < 0 && mp->mp2nd->viscositymask[0])) {
double dq_gradp2, pre_delP2, dq_dH2, srate2, qmag_log;
k_turb = 12.;
dq_gradp2 = pre_delP2 = -CUBE(H) / (k_turb * mp->mp2nd->viscosity);
q_mag2 = pre_delP2 * pgrad;
dq_dH2 = -3. * SQUARE(H) / (k_turb * mp->mp2nd->viscosity) * pgrad;
srate2 = tau_w / mp->mp2nd->viscosity;
qmag_log = (DOUBLE_NONZERO(q_mag) ? log(q_mag2 / q_mag) : 0.0);
if (fabs(fv->F) > ls->Length_Scale) {
q_mag = q_mag2;
dq_gradp = dq_gradp2;
pre_delP = pre_delP2;
dq_dH = dq_dH2;
srate = srate2;
vis_w = mp->mp2nd->viscosity;
} else {
factor = (mp->mp2nd->viscositymask[1] ? (1.0 - lsi->H) : lsi->H);
q_mag = -pow(-q_mag, factor) * pow(-q_mag2, 1.0 - factor);
dq_gradp = -pow(-dq_gradp, factor) * pow(-dq_gradp2, 1.0 - factor);
pre_delP = -pow(-pre_delP, factor) * pow(-pre_delP2, 1.0 - factor);
dq_dH = -pow(-dq_dH, factor) * pow(-dq_dH2, 1.0 - factor);
srate = pow(srate, factor) * pow(srate2, 1.0 - factor);
vis_w = pow(vis_w, factor) * pow(mp->mp2nd->viscosity, 1.0 - factor);
for (j = 0; j < ei[pg->imtrx]->dof[VAR]; j++) {
dqmag_dF[j] += q_mag * qmag_log * lsi->d_H_dF[j];
}
}
}
} else if (mp->Lub_LS_Interpolation == LINEAR) {
double dq_gradp2, pre_delP2, dq_dH2, srate2;
k_turb = 12.;
dq_gradp2 = pre_delP2 = -CUBE(H) / (k_turb * mp->mp2nd->viscosity);
q_mag2 = pre_delP2 * pgrad;
dq_dH2 = -3. * SQUARE(H) / (k_turb * mp->mp2nd->viscosity) * pgrad;
srate2 = tau_w / mp->mp2nd->viscosity;
q_mag = ls_modulate_property(q_mag, q_mag2, ls->Length_Scale,
(double)mp->mp2nd->viscositymask[0],
(double)mp->mp2nd->viscositymask[1], dqmag_dF, &factor);
dq_gradp = dq_gradp * factor + dq_gradp2 * (1. - factor);
pre_delP = pre_delP * factor + pre_delP2 * (1. - factor);
dq_dH = dq_dH * factor + (1. - factor) * dq_dH2;
dq_dT *= factor; // mp2nd->viscosity is independent of Temperature
srate = srate * factor + (1. - factor) * srate2;
vis_w = vis_w * factor + (1. - factor) * mp->mp2nd->viscosity;
} else {
GOMA_WH(GOMA_ERROR, "mp->Lub_LS_Interpolation needs to be LOG or LINEAR...\n");
}
} else {
GOMA_WH(GOMA_ERROR, "mp2nd->ViscosityModel needs to be RATIO or CONSTANT...\n");
}
Expand Down Expand Up @@ -5555,11 +5626,35 @@ void calculate_lub_q_v_old(
(double)mp->mp2nd->viscositymask[0],
(double)mp->mp2nd->viscositymask[1], dqmag_dF, &factor);
} else if (mp->mp2nd->ViscosityModel == CONSTANT || mp->mp2nd->ViscosityModel == TIME_RAMP) {
k_turb = 12.;
q_mag2 = -CUBE(H_old) / (k_turb * mp->mp2nd->viscosity) * pgrad;
q_mag = ls_modulate_property(q_mag, q_mag2, ls->Length_Scale,
(double)mp->mp2nd->viscositymask[0],
(double)mp->mp2nd->viscositymask[1], dqmag_dF, &factor);
if (mp->Lub_LS_Interpolation == LOGARITHMIC) {
if (lsi->near || (fv->F > 0 && mp->mp2nd->viscositymask[1]) ||
(fv->F < 0 && mp->mp2nd->viscositymask[0])) {
double pre_delP2, qmag_log;
k_turb = 12.;
pre_delP2 = -CUBE(H_old) / (k_turb * mp->mp2nd->viscosity);
q_mag2 = pre_delP2 * pgrad;
qmag_log = (DOUBLE_NONZERO(q_mag) ? log(q_mag2 / q_mag) : 0.0);
if (fabs(fv->F) > ls->Length_Scale) {
q_mag = q_mag2;
} else {
factor = (mp->mp2nd->viscositymask[1] ? (1.0 - lsi->H) : lsi->H);
q_mag = -pow(-q_mag, factor) * pow(-q_mag2, 1.0 - factor);
for (j = 0; j < ei[pg->imtrx]->dof[VAR]; j++) {
dqmag_dF[j] += q_mag * qmag_log * lsi->d_H_dF[j];
}
}
}
} else if (mp->Lub_LS_Interpolation == LINEAR) {
double pre_delP2;
k_turb = 12.;
pre_delP2 = -CUBE(H_old) / (k_turb * mp->mp2nd->viscosity);
q_mag2 = pre_delP2 * pgrad;
q_mag = ls_modulate_property(q_mag, q_mag2, ls->Length_Scale,
(double)mp->mp2nd->viscositymask[0],
(double)mp->mp2nd->viscositymask[1], dqmag_dF, &factor);
} else {
GOMA_WH(GOMA_ERROR, "mp->Lub_LS_Interpolation needs to be LOG or LINEAR...\n");
}
} else {
GOMA_WH(GOMA_ERROR, "mp2nd->ViscosityModel needs to be RATIO or CONSTANT...\n");
}
Expand Down Expand Up @@ -6601,24 +6696,24 @@ int lub_viscosity_integrate(const double strs,
switch (gn->ConstitutiveEquation) {
case BINGHAM:
case BINGHAM_WLF: {
double shrF = 1. / F, shrY = pow(yield * pow(lam, 1. - nexp) / (mu0 - muinf), 1. / nexp);
double shrF = 1. / F; /* the below assumes muinf = 0 */
double shrY = pow((yield + mu0 * shrF) * pow(lam, 1. - nexp) / mu0, 1. / nexp);
double tmp_pl3 = 1. / CUBE(lam * shrw);
double tmp_cy1 = pow(1. + CUBE(lam * shrw), (2. + nexp) / 3.);
double tmp_cyY = pow(1. + CUBE(lam * shrY), (2. * nexp) / 3.);
double tmp_cy2 = pow(1. + CUBE(lam * shrw), (1. + 2. * nexp) / 3.);
double tmp_cyY2 = pow(1. + CUBE(lam * shrY), (1. + 2. * nexp) / 3.);
double shr1 = MIN(shrw, shrF), shr2 = MIN(shrw, shrY);

xint_a = SQUARE(F * yield + mu0) * CUBE(shr1 / shrw) / 3.;
xint_a = SQUARE(F * yield) * CUBE(shr1 / shrw) / 3.;
if (shrw > shrF) {
xint_a += (SQUARE(yield) * (shr2 - shrF) + mu0 * yield * (SQUARE(shr2) - SQUARE(shrF)) +
SQUARE(mu0) / 3. * (CUBE(shr2) - CUBE(shrF))) /
CUBE(shrw);
xint_a += SQUARE(yield + mu0 * shrF) * (shr2 - shrF);
}
xint_a /= CUBE(shrw);
if (shrw > shrY) {
xint_a += SQUARE(muinf) / 3. * (1. - CUBE(shrY / shrw)) +
2. * muinf * tmp_pl3 * (mu0 - muinf) / (2. + nexp) * (tmp_cy1 - tmp_cyY) +
tmp_pl3 * SQUARE(mu0 - muinf) / (1. + 2. * nexp) * (tmp_cy2 - tmp_cyY2);
xint_a += SQUARE(muinf) / 3. * tmp_pl3 * (CUBE(shrw) - CUBE(shrY)) +
2. * muinf * (mu0 - muinf) / (2. + nexp) * tmp_pl3 * (tmp_cy1 - tmp_cyY) +
SQUARE(mu0 - muinf) / (1. + 2. * nexp) * tmp_pl3 * (tmp_cy2 - tmp_cyY2);
}
} break;
case CARREAU:
Expand All @@ -6634,8 +6729,6 @@ int lub_viscosity_integrate(const double strs,
for (jdi = 0; jdi < JDI_MAX; jdi++) {
double cee, x0, delx, vis = 1., jdiv, xfact, tmp, tpe, tp2, P_sig;
int idiv, l;
double shrF = 1. / F;
double shrY = pow(yield * pow(lam, 1. - nexp) / (mu0 - muinf), 1. / nexp);
jdiv = pow(2., jdi);
delx = 1. / jdiv;
x0 = 0.0;
Expand All @@ -6652,19 +6745,21 @@ int lub_viscosity_integrate(const double strs,
visc_a = muinf + (mu0 - muinf) / pow(1. + CUBE(cee * lam * shrw), (1. - nexp) / 3.);
break;
case BINGHAM:
case BINGHAM_WLF:
case BINGHAM_WLF: {
double shrF = 1. / F;
double shrY = pow((yield + mu0 * shrF) * pow(lam, 1. - nexp) / mu0, 1. / nexp);
tp2 = F * shrw;
P_sig = pow(1. + tp2, P_eps);
tpe = (1. - exp(-tp2)) / shrw * P_sig;
vis = muinf + (mu0 - muinf + yield * tpe) * tmp;
if (cee * shrw < shrF) {
visc_a = F * yield + mu0;
} else if (cee * shrw < shrY) {
visc_a = mu0 + yield / (cee * shrw);
visc_a = (mu0 * shrF + yield) / (cee * shrw);
} else {
visc_a = muinf + (mu0 - muinf) / pow(1. + CUBE(cee * lam * shrw), (1. - nexp) / 3.);
}
break;
} break;
default:
GOMA_EH(GOMA_ERROR, "Missing Lub Viscosity model!");
}
Expand Down

0 comments on commit 86061fd

Please sign in to comment.