From 86061fde8d806def890e181fdd640e326c02e291 Mon Sep 17 00:00:00 2001 From: Robert Secor Date: Wed, 24 Jul 2024 14:28:16 -0500 Subject: [PATCH] Lubrication Logarithm Interpolation (#468) 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. --- include/mm_mp_structs.h | 1 + include/rf_fem_const.h | 1 + src/dp_vif.c | 1 + src/mm_fill_species.c | 1 - src/mm_input_mp.c | 20 ++++- src/mm_shell_util.c | 187 ++++++++++++++++++++++++++++++---------- 6 files changed, 163 insertions(+), 48 deletions(-) diff --git a/include/mm_mp_structs.h b/include/mm_mp_structs.h index ba3d0995..4aa23a2b 100755 --- a/include/mm_mp_structs.h +++ b/include/mm_mp_structs.h @@ -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; diff --git a/include/rf_fem_const.h b/include/rf_fem_const.h index 46cdcd5f..7464cfd5 100644 --- a/include/rf_fem_const.h +++ b/include/rf_fem_const.h @@ -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 diff --git a/src/dp_vif.c b/src/dp_vif.c index da433544..e65d6fd8 100644 --- a/src/dp_vif.c +++ b/src/dp_vif.c @@ -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); diff --git a/src/mm_fill_species.c b/src/mm_fill_species.c index 88f42e25..d77acaec 100644 --- a/src/mm_fill_species.c +++ b/src/mm_fill_species.c @@ -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++) { diff --git a/src/mm_input_mp.c b/src/mm_input_mp.c index 789958ea..866d4e7f 100644 --- a/src/mm_input_mp.c +++ b/src/mm_input_mp.c @@ -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); @@ -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]) { diff --git a/src/mm_shell_util.c b/src/mm_shell_util.c index 12f0c0f8..422a7e23 100644 --- a/src/mm_shell_util.c +++ b/src/mm_shell_util.c @@ -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; @@ -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"); } @@ -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"); } @@ -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"); } @@ -6601,7 +6696,8 @@ 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.); @@ -6609,16 +6705,15 @@ int lub_viscosity_integrate(const double strs, 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: @@ -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; @@ -6652,7 +6745,9 @@ 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; @@ -6660,11 +6755,11 @@ int lub_viscosity_integrate(const double strs, 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!"); }