From 8d859361ac7cfea257d53e5872a8ab6977a0f054 Mon Sep 17 00:00:00 2001 From: Robert Secor Date: Wed, 28 Aug 2024 08:42:57 -0500 Subject: [PATCH] Lub log ls2 (#470) * git fetch update * Curvature updates * Curvature updates * Curvature updates * lsi_near * lsi_near * Curvature updates * Curvature updates * Curvature updates * Curvature updates * Curvature updates * Curvature updates * Curvature updates * Curvature updates * Curvature updates * Curvature updates * Curvature updates * Shell Energy * Shell Energy * Lub Temp Conv * Lub Temp Conv * Lub Temp Conv * Lub Temp Conv * Lub Temp Conv * Lub Temp Conv * Lub Temp Conv * shell species update * shell species update * shell species update * shell species update * shell species update * shell species update * shell species update * shell species update * shell species update * shell species update * shell species update * shell species update * shell species update * Lub LogLS * Lub LogLS * Lub LogLS * Lub LogLS * more LS_Log * more LS_Log * more LS_Log * more LS_Log * more LS_Log * more LS_Log * more LS_Log * some shell temp updates * some shell temp updates * some shell temp updates * some shell temp updates * some shell temp updates * some shell temp updates * some shell temp updates * some shell temp updates * some shell temp updates * some shell temp updates * Fix some address sanitizer issues, maybe d_LSnormal_dF will change convergence --------- Co-authored-by: Robert Secor Co-authored-by: Robert Secor Co-authored-by: Weston Ortiz --- include/mm_as_structs.h | 7 +++ include/mm_mp_structs.h | 1 + include/mm_post_def.h | 6 +- include/rf_fem_const.h | 2 +- src/dp_vif.c | 1 + src/mm_fill_shell.c | 85 +++++++++++++++++++--------- src/mm_flux.c | 116 +++++++++++++++++++++++--------------- src/mm_input_mp.c | 50 +++++++++++++--- src/mm_shell_bc.c | 23 ++++++-- src/mm_shell_util.c | 76 ++++++++++++++++++------- src/mm_std_models_shell.c | 34 +++++++++-- 11 files changed, 285 insertions(+), 116 deletions(-) diff --git a/include/mm_as_structs.h b/include/mm_as_structs.h index 9239b992..e2b201e7 100644 --- a/include/mm_as_structs.h +++ b/include/mm_as_structs.h @@ -3222,8 +3222,11 @@ struct Lubrication_Auxiliaries { double gradP_tangent[DIM]; /* Tangent vector of the pressure gradient */ double gradP_normal[DIM]; /* Unit vector perpendicular to the pressure */ double H; /* Lubrication Gap Height */ + double H_cap; /* Lubrication Gap Height unaffected by wall effects */ double srate; /* Lubrication Characteristic Shear Rate */ double mu_star; /* Lubrication Characteristic Viscosity */ + double op_curv; /* Lubrication Out-of-plane Curvature */ + double visc_diss; /* Lubrication Integrated Viscous Dissipation */ double dgradP_mag_dP; /* Pressure gradient magnitude sensitivities w.r.t. pressure */ @@ -3289,6 +3292,10 @@ struct Lubrication_Auxiliaries { solid */ double dH_dp; /* lubrication gap sensitivities w.r.t. pressure */ double dH_ddh; /* lubrication gap sensitivities w.r.t. added height */ + double dop_curv_dx[DIM][MDE]; /* Out-of-plane Curvature sensitivities w.r.t. mesh deformation */ + double dop_curv_df[MDE]; /* Out-of-plane Curvature sensitivities w.r.t. level set */ + double dvisc_diss_dT; /* Lubrication Integrated Viscous Dissipation Sensitivity to Temperature */ + double dvisc_diss_dpgrad; /* Lubrication Integrated Viscous Dissipation Sensitivity to Pgrad */ }; typedef struct Lubrication_Auxiliaries LUBRICATION_AUXILIARIES_STRUCT; diff --git a/include/mm_mp_structs.h b/include/mm_mp_structs.h index 4aa23a2b..7d20aaa0 100755 --- a/include/mm_mp_structs.h +++ b/include/mm_mp_structs.h @@ -774,6 +774,7 @@ struct Material_Properties { int Lub_Curv_MassLump; int Lub_Curv_Modulation; int Lub_LS_Interpolation; + int Lub_Curv_Combine; int Lub_Heat_XferModel; double Lub_Heat_Xfer; diff --git a/include/mm_post_def.h b/include/mm_post_def.h index bef7c8ba..5eca1ef6 100644 --- a/include/mm_post_def.h +++ b/include/mm_post_def.h @@ -136,6 +136,8 @@ #define I_POROUS_LIQUID_INV_2 49 #define I_POROUS_LIQUID_INV_3 50 #define I_EM_ABSORB_CROSS_SECTION 51 +#define I_HEAT_ENERGY_NEG 52 +#define I_HEAT_ENERGY_POS 53 #ifdef GOMA_MM_POST_PROC_C struct Post_Processing_Flux_Names { char *name; /* flux string */ @@ -263,7 +265,9 @@ VOL_NAME_STRUCT pp_vol_names[] = {{"VOLUME", I_VOLUME}, {"LAMB_MAG", I_LAMB_MAG}, {"HELICITY", I_HELICITY}, {"Q_FCN", I_Q_FCN}, - {"ABSORPTION_CROSS_SECTION", I_EM_ABSORB_CROSS_SECTION}}; + {"ABSORPTION_CROSS_SECTION", I_EM_ABSORB_CROSS_SECTION}, + {"HEAT_ENERGY_NEG", I_HEAT_ENERGY_NEG}, + {"HEAT_ENERGY_POS", I_HEAT_ENERGY_POS}}; int Num_Vol_Names = sizeof(pp_vol_names) / sizeof(VOL_NAME_STRUCT); diff --git a/include/rf_fem_const.h b/include/rf_fem_const.h index fa4b9f02..f1ead3bc 100644 --- a/include/rf_fem_const.h +++ b/include/rf_fem_const.h @@ -149,7 +149,7 @@ #define LUB_VISCINT_GAUSSIAN 25 #define LUB_VISCINT_ANALYTICAL 26 #define LUB_VISCINT_POWERLAW 27 -#define MAX_LUB_NGP 5 +#define MAX_LUB_NGP 6 #define LOGARITHMIC 29 /* Residence time kernel functions */ diff --git a/src/dp_vif.c b/src/dp_vif.c index f99ddb67..60cb99bf 100644 --- a/src/dp_vif.c +++ b/src/dp_vif.c @@ -1590,6 +1590,7 @@ void noahs_ark(void) { 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]->Lub_Curv_Combine, 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_shell.c b/src/mm_fill_shell.c index 867d2001..64dd647c 100644 --- a/src/mm_fill_shell.c +++ b/src/mm_fill_shell.c @@ -7247,9 +7247,6 @@ int assemble_shell_energy(double time, /* present time value */ /* * Load material property constants */ -#if 0 - mu = viscosity(gn, NULL, d_mu); -#endif rho = density(d_rho, time); @@ -7271,7 +7268,8 @@ int assemble_shell_energy(double time, /* present time value */ if (pd->gv[R_LUBP]) { calculate_lub_q_v(R_LUBP, time, dt, xi, exo); } - H = LubAux->H; + /* Let's use the lubrication gap unaffected by wall effects... */ + H = LubAux->H_cap; k_eff = t_cond; /* default laminar case */ d_k_eff_dh = 0.; @@ -7328,7 +7326,11 @@ int assemble_shell_energy(double time, /* present time value */ if (pd->v[pg->imtrx][FILL]) { load_lsi(ls->Length_Scale); load_lsi_derivs(); - if_liquid = (1. - lsi->H); + if (mp->mp2nd->viscositymask[1]) { + if_liquid = (1. - lsi->H); + } else { + if_liquid = lsi->H; + } } /* Temperature gradient */ @@ -7375,10 +7377,14 @@ int assemble_shell_energy(double time, /* present time value */ /* Finall, Load up all heat source models */ /* Note all of this depend on H_lub. You need to add those sensitivities */ - /* No source terms available right now. Feel free to add your own. */ + /* No source terms available right now. Feel free to add your own. + Added heat energy is negative, so I guess this is heat loss from the channel */ const double ht_coeff = mp->Lub_Heat_Xfer; const double ht_tamb = mp->Lub_Heat_Tamb; q_tot = ht_coeff * (fv->sh_t - ht_tamb); + if (mp->HeatSourceModel == VISC_DISS) { + q_tot -= mp->u_heat_source[0] * LubAux->visc_diss; + } /*** RESIDUAL ASSEMBLY ******************************************************/ if (af->Assemble_Residual) { @@ -7403,7 +7409,7 @@ int assemble_shell_energy(double time, /* present time value */ if (pd->TimeIntegration != STEADY) { if (pd->e[pg->imtrx][eqn] & T_MASS) { mass = fv_dot->sh_t; - mass *= -H * phi_i * rho * Cp * det_J * wt; + mass *= H * phi_i * rho * Cp * det_J * wt; mass *= h3; mass *= pd->etm[pg->imtrx][eqn][(LOG2_MASS)]; } @@ -7424,7 +7430,7 @@ int assemble_shell_energy(double time, /* present time value */ advection += LubAux->q[p] * grad_T[p]; } - advection *= -rho * Cp * det_J * wt_func * wt; + advection *= rho * Cp * det_J * wt_func * wt; advection *= h3; advection *= pd->etm[pg->imtrx][eqn][(LOG2_ADVECTION)]; } @@ -7525,7 +7531,7 @@ int assemble_shell_energy(double time, /* present time value */ if (pd->TimeIntegration != STEADY) { if (pd->e[pg->imtrx][eqn] & T_MASS) { mass = (1 + 2. * tt) * phi_j / dt; - mass *= -H * phi_i * rho * Cp * det_J * wt; + mass *= H * phi_i * rho * Cp * det_J * wt; mass *= h3 * pd->etm[pg->imtrx][eqn][(LOG2_MASS)]; } } @@ -7539,7 +7545,7 @@ int assemble_shell_energy(double time, /* present time value */ advection += LubAux->dq_dT[a] * grad_T[a] * phi_j; } - advection *= -rho * Cp * det_J * wt * wt_func; + advection *= rho * Cp * det_J * wt * wt_func; advection *= h3; advection *= pd->etm[pg->imtrx][eqn][(LOG2_ADVECTION)]; } @@ -7558,7 +7564,9 @@ int assemble_shell_energy(double time, /* present time value */ if (pd->e[pg->imtrx][eqn] & T_SOURCE) { /* PRS Note these are initialized to zero in mm_input_mp.c */ source = 0.; /* add your dq_dt sensitivities here */ - source = ht_coeff * phi_j; + if (mp->HeatSourceModel == VISC_DISS) { + source = mp->u_heat_source[0] * (-LubAux->dvisc_diss_dT + ht_coeff) * phi_j; + } } source *= phi_i * if_liquid * det_J * wt * h3 * pd->etm[pg->imtrx][eqn][(LOG2_SOURCE)]; @@ -7592,7 +7600,7 @@ int assemble_shell_energy(double time, /* present time value */ if (pd->TimeIntegration != STEADY) { if (pd->e[pg->imtrx][eqn] & T_MASS) { mass = fv_dot->sh_t; - mass *= -H * phi_i * det_J * wt * (rho * d_Cp->F[j] + d_rho->F[j] * Cp); + mass *= H * phi_i * det_J * wt * (rho * d_Cp->F[j] + d_rho->F[j] * Cp); mass *= h3 * pd->etm[pg->imtrx][eqn][(LOG2_MASS)]; } } @@ -7605,7 +7613,7 @@ int assemble_shell_energy(double time, /* present time value */ advection += LubAux->dq_df[a][j] * grad_T[a]; } - advection *= -(rho * d_Cp->F[j] + d_rho->F[j] * Cp) * det_J * wt * wt_func; + advection *= rho * d_Cp->F[j] + d_rho->F[j] * Cp * det_J * wt * wt_func; advection *= h3; advection *= pd->etm[pg->imtrx][eqn][(LOG2_ADVECTION)]; } @@ -7681,7 +7689,7 @@ int assemble_shell_energy(double time, /* present time value */ if (pd->TimeIntegration != STEADY) { if (pd->e[pg->imtrx][eqn] & T_MASS) { mass = fv_dot->sh_t; - mass *= (-LubAux->dH_dmesh[b][jk] * phi_i * rho * Cp * det_J * wt - + mass *= (LubAux->dH_dmesh[b][jk] * phi_i * rho * Cp * det_J * wt + H * phi_i * rho * Cp * fv->dsurfdet_dx[b][jk] * wt); mass *= h3; mass *= pd->etm[pg->imtrx][eqn][(LOG2_MASS)]; @@ -7698,7 +7706,7 @@ int assemble_shell_energy(double time, /* present time value */ advection += fv->dsurfdet_dx[b][jk] * LubAux->q[p] * grad_T[p]; } - advection *= -rho * Cp * wt * wt_func; + advection *= rho * Cp * wt * wt_func; advection *= h3; if (supg != 0.) { @@ -7711,7 +7719,7 @@ int assemble_shell_energy(double time, /* present time value */ LubAux->dH_dmesh[b][jk]); } for (p = 0; p < dim; p++) { - advection += (-rho * Cp * d_wt_func * h3 * det_J * wt) * LubAux->q[p] * grad_T[p]; + advection += rho * Cp * d_wt_func * h3 * det_J * wt * LubAux->q[p] * grad_T[p]; } } advection *= pd->etm[pg->imtrx][eqn][(LOG2_ADVECTION)]; @@ -7789,7 +7797,7 @@ int assemble_shell_energy(double time, /* present time value */ if (pd->TimeIntegration != STEADY) { if (pd->e[pg->imtrx][eqn] & T_MASS) { mass = fv_dot->sh_t; - mass *= (-LubAux->dH_drealsolid[b][jk] * phi_i * rho * Cp * det_J * wt); + mass *= LubAux->dH_drealsolid[b][jk] * phi_i * rho * Cp * det_J * wt; mass *= h3; mass *= pd->etm[pg->imtrx][eqn][(LOG2_MASS)]; } @@ -7803,7 +7811,7 @@ int assemble_shell_energy(double time, /* present time value */ advection += det_J * LubAux->dq_drs[p][b][jk] * grad_T[p]; } - advection *= -rho * Cp * wt * wt_func; + advection *= rho * Cp * wt * wt_func; advection *= h3; if (supg != 0.) { @@ -7813,7 +7821,7 @@ int assemble_shell_energy(double time, /* present time value */ LubAux->dH_drealsolid[b][jk]); } for (p = 0; p < dim; p++) { - advection += (-rho * Cp * d_wt_func * h3 * det_J * wt) * LubAux->q[p] * grad_T[p]; + advection += rho * Cp * d_wt_func * h3 * det_J * wt * LubAux->q[p] * grad_T[p]; } } advection *= pd->etm[pg->imtrx][eqn][(LOG2_ADVECTION)]; @@ -7887,7 +7895,7 @@ int assemble_shell_energy(double time, /* present time value */ } for (p = 0; p < dim; p++) { - advection += (-rho * Cp * d_wt_func * h3 * det_J * wt) * LubAux->q[p] * grad_T[p]; + advection += rho * Cp * d_wt_func * h3 * det_J * wt * LubAux->q[p] * grad_T[p]; } } advection *= pd->etm[pg->imtrx][eqn][(LOG2_ADVECTION)]; @@ -7927,7 +7935,7 @@ int assemble_shell_energy(double time, /* present time value */ if (pd->TimeIntegration != STEADY) { if (pd->e[pg->imtrx][eqn] & T_MASS) { mass = fv_dot->sh_t; - mass *= (-phi_j * phi_i * rho * Cp * det_J * wt); + mass *= phi_j * phi_i * rho * Cp * det_J * wt; mass *= h3; mass *= pd->etm[pg->imtrx][eqn][(LOG2_MASS)]; } @@ -7949,7 +7957,7 @@ int assemble_shell_energy(double time, /* present time value */ advection += LubAux->q[p] * grad_T[p] * d_wt_func; } - advection *= -rho * Cp * det_J * wt; + advection *= rho * Cp * det_J * wt; advection *= h3; advection *= pd->etm[pg->imtrx][eqn][(LOG2_ADVECTION)]; } @@ -7989,13 +7997,21 @@ int assemble_shell_energy(double time, /* present time value */ /* Load basis functions (j) */ phi_j = bf[eqn]->phi[j]; + for (ii = 0; ii < VIM; ii++) { + grad_phi_j[ii] = bf[eqn]->grad_phi[j][ii]; + grad_II_phi_j[ii] = 0.0; + for (jj = 0; jj < VIM; jj++) { + grad_II_phi_j[ii] += (grad_phi_j[jj] * delta(ii, jj) - + grad_phi_j[jj] * (fv->snormal[ii] * fv->snormal[jj])); + } + } /* Mass term */ mass = 0.; if (pd->TimeIntegration != STEADY) { if (pd->e[pg->imtrx][eqn] & T_MASS) { mass = fv_dot->sh_t; - mass *= -LubAux->dH_dp * phi_j * phi_i * rho * Cp * det_J * wt; + mass *= LubAux->dH_dp * phi_j * phi_i * rho * Cp * det_J * wt; mass *= h3; mass *= pd->etm[pg->imtrx][eqn][(LOG2_MASS)]; } @@ -8012,7 +8028,7 @@ int assemble_shell_energy(double time, /* present time value */ } } - advection *= -rho * Cp * det_J * wt * wt_func; + advection *= rho * Cp * det_J * wt * wt_func; advection *= h3; advection *= pd->etm[pg->imtrx][eqn][(LOG2_ADVECTION)]; @@ -8032,7 +8048,7 @@ int assemble_shell_energy(double time, /* present time value */ grad_II_phi_i[p]); } } - advection_b *= -rho * Cp * det_J * wt; + advection_b *= rho * Cp * det_J * wt; } advection += advection_b; } @@ -8049,7 +8065,13 @@ int assemble_shell_energy(double time, /* present time value */ } source = 0.0; if (pd->e[pg->imtrx][eqn] & T_SOURCE) { - source += 0.; /* Add on your dq_tot_d_lubp terms here */ + /* Viscous dissipation term depends on pgrad and q_mag */ + if (mp->HeatSourceModel == VISC_DISS) { + for (p = 0; p < dim; p++) { + source += fv->grad_lubp[p] * grad_II_phi_j[p]; + } + source *= -mp->u_heat_source[0] * LubAux->dvisc_diss_dpgrad / LubAux->gradP_mag; + } } source *= phi_i * if_liquid * wt * det_J * h3 * pd->etm[pg->imtrx][eqn][(LOG2_SOURCE)]; @@ -13984,6 +14006,9 @@ int assemble_lubrication_curvature(double time, /* present time value for (a = 0; a < VIM; a++) { div += LSnormal[a] * grad_II_phi_i[a]; } + if (mp->Lub_Curv_Combine) { + div += LubAux->op_curv * phi_i; + } } div *= curvX * det_J * wt * h3 * pd->etm[pg->imtrx][eqn][(LOG2_DIVERGENCE)]; } @@ -14002,7 +14027,7 @@ int assemble_lubrication_curvature(double time, /* present time value /* Loop over DOFs (i) */ for (i = 0; i < ei[pg->imtrx]->dof[eqn]; i++) { - /* Prepare basis funcitons (i) */ + /* Prepare basis functions (i) */ ShellBF(eqn, i, &phi_i, grad_phi_i, grad_II_phi_i, d_grad_II_phi_i_dmesh, n_dof[MESH_DISPLACEMENT1], dof_map); @@ -14122,6 +14147,9 @@ int assemble_lubrication_curvature(double time, /* present time value div += LSnormal[a] * d_grad_II_phi_i_dmesh[a][b][jj] * det_J; div += LSnormal[a] * grad_II_phi_i[a] * fv->dsurfdet_dx[b][jj]; } + if (mp->Lub_Curv_Combine) { + div += LubAux->dop_curv_dx[b][jj] * phi_i; + } div *= curvX * wt * h3 * pd->etm[pg->imtrx][eqn][(LOG2_DIVERGENCE)]; } } @@ -14154,6 +14182,9 @@ int assemble_lubrication_curvature(double time, /* present time value for (a = 0; a < VIM; a++) { div += d_LSnormal_dF[a][j] * grad_II_phi_i[a]; } + if (mp->Lub_Curv_Combine) { + div += LubAux->dop_curv_df[j] * phi_i; + } div *= curvX * det_J * wt * h3 * pd->etm[pg->imtrx][eqn][(LOG2_DIVERGENCE)]; } } diff --git a/src/mm_flux.c b/src/mm_flux.c index 07e46926..de2a97c1 100644 --- a/src/mm_flux.c +++ b/src/mm_flux.c @@ -4156,36 +4156,24 @@ double evaluate_volume_integral(const Exo_DB *exo, /* ptr to basic exodus ii me int side_id[12]; double xf2D[6][2]; int nint2D[6]; - - adaptive_integration_active = - (ls != NULL && ls->AdaptIntegration && - (quantity == I_POS_FILL || quantity == I_NEG_FILL || quantity == I_MASS_NEGATIVE_FILL || - quantity == I_MASS_POSITIVE_FILL || quantity == I_POS_VOLPLANE || - quantity == I_NEG_VOLPLANE || quantity == I_POS_CENTER_X || quantity == I_POS_CENTER_Y || - quantity == I_POS_CENTER_Z || quantity == I_NEG_CENTER_X || quantity == I_NEG_CENTER_Y || - quantity == I_NEG_CENTER_Z || quantity == I_POS_VX || quantity == I_POS_VY || - quantity == I_POS_VZ || quantity == I_NEG_VX || quantity == I_NEG_VY || - quantity == I_NEG_VZ)); - + int neg_LS_integ = + (quantity == I_NEG_FILL || quantity == I_MASS_NEGATIVE_FILL || quantity == I_NEG_VOLPLANE || + quantity == I_NEG_VX || quantity == I_NEG_VY || quantity == I_NEG_VZ || + quantity == I_NEG_CENTER_X || quantity == I_NEG_CENTER_Y || quantity == I_NEG_CENTER_Z || + I_HEAT_ENERGY_NEG); + int pos_LS_integ = + (quantity == I_POS_FILL || quantity == I_MASS_POSITIVE_FILL || quantity == I_POS_VOLPLANE || + quantity == I_POS_VX || quantity == I_POS_VY || quantity == I_POS_VZ || + quantity == I_POS_CENTER_X || quantity == I_POS_CENTER_Y || quantity == I_POS_CENTER_Z || + I_HEAT_ENERGY_POS); + int do_LS_integ = (neg_LS_integ || pos_LS_integ); + + adaptive_integration_active = (ls != NULL && ls->AdaptIntegration && do_LS_integ); subgrid_integration_active = - (ls != NULL && ls->Integration_Depth > 0 && !adaptive_integration_active && - (quantity == I_POS_FILL || quantity == I_NEG_FILL || quantity == I_MASS_NEGATIVE_FILL || - quantity == I_MASS_POSITIVE_FILL || quantity == I_POS_VOLPLANE || - quantity == I_NEG_VOLPLANE || quantity == I_POS_CENTER_X || quantity == I_POS_CENTER_Y || - quantity == I_POS_CENTER_Z || quantity == I_NEG_CENTER_X || quantity == I_NEG_CENTER_Y || - quantity == I_NEG_CENTER_Z || quantity == I_LS_ARC_LENGTH || quantity == I_POS_VX || - quantity == I_POS_VY || quantity == I_POS_VZ || quantity == I_NEG_VX || - quantity == I_NEG_VY || quantity == I_NEG_VZ)); + (ls != NULL && ls->Integration_Depth > 0 && !adaptive_integration_active && do_LS_integ); subelement_vol_integration_active = (ls != NULL && ls->SubElemIntegration && !adaptive_integration_active && - (params == NULL || params[0] == 0.) && - (quantity == I_POS_FILL || quantity == I_NEG_FILL || quantity == I_MASS_NEGATIVE_FILL || - quantity == I_MASS_POSITIVE_FILL || quantity == I_POS_VOLPLANE || - quantity == I_NEG_VOLPLANE || quantity == I_POS_CENTER_X || quantity == I_POS_CENTER_Y || - quantity == I_POS_CENTER_Z || quantity == I_NEG_CENTER_X || quantity == I_NEG_CENTER_Y || - quantity == I_NEG_CENTER_Z || quantity == I_POS_VX || quantity == I_POS_VY || - quantity == I_POS_VZ || quantity == I_NEG_VX || quantity == I_NEG_VY || - quantity == I_NEG_VZ)); + (params == NULL || params[0] == 0.) && do_LS_integ); subelement_surf_integration_active = (ls != NULL && ls->SubElemIntegration && !adaptive_integration_active && (params == NULL || params[0] == 0.) && @@ -4246,7 +4234,7 @@ double evaluate_volume_integral(const Exo_DB *exo, /* ptr to basic exodus ii me ip_total = elem_info(NQUAD, ei[pg->imtrx]->ielem_type); if (subgrid_integration_active) { - double width = params == NULL ? ls->Length_Scale : 2.0 * params[0]; + double width = num_params == 0 ? ls->Length_Scale : 2.0 * params[0]; if ((Use_Subgrid_Integration = current_elem_overlaps_interface(width))) { ip_total = @@ -4254,15 +4242,9 @@ double evaluate_volume_integral(const Exo_DB *exo, /* ptr to basic exodus ii me } } else if (subelement_vol_integration_active) { if ((Use_Subelement_Integration = current_elem_on_isosurface(FILL, 0.))) { - if (quantity == I_NEG_FILL || quantity == I_MASS_NEGATIVE_FILL || - quantity == I_NEG_VOLPLANE || quantity == I_NEG_VX || quantity == I_NEG_VY || - quantity == I_NEG_VZ || quantity == I_NEG_CENTER_X || quantity == I_NEG_CENTER_Y || - quantity == I_NEG_CENTER_Z) { + if (neg_LS_integ) { ip_total = get_subelement_integration_pts(&s, &weight, NULL, 0., -2, -1); - } else if (quantity == I_POS_FILL || quantity == I_MASS_POSITIVE_FILL || - quantity == I_POS_VOLPLANE || quantity == I_POS_VX || quantity == I_POS_VY || - quantity == I_POS_VZ || quantity == I_POS_CENTER_X || - quantity == I_POS_CENTER_Y || quantity == I_POS_CENTER_Z) { + } else if (pos_LS_integ) { ip_total = get_subelement_integration_pts(&s, &weight, NULL, 0., -2, 1); } else { ip_total = get_subelement_integration_pts(&s, &weight, NULL, 0., -2, 0); @@ -4281,9 +4263,7 @@ double evaluate_volume_integral(const Exo_DB *exo, /* ptr to basic exodus ii me ls_F[i] = *esp_old->F[i]; } wt_type = 2; - if (quantity == I_NEG_FILL || quantity == I_NEG_VOLPLANE || quantity == I_NEG_CENTER_X || - quantity == I_NEG_CENTER_Y || quantity == I_NEG_CENTER_Z || quantity == I_NEG_VX || - quantity == I_NEG_VY || quantity == I_NEG_VZ) { + if (neg_LS_integ) { for (i = 0; i < ei[pg->imtrx]->num_local_nodes; i++) { ls_F[i] = -ls_F[i]; } @@ -5060,8 +5040,22 @@ int compute_volume_integrand(const int quantity, } } break; - case I_HEAT_ENERGY: { - double rho, Cp; + case I_HEAT_ENERGY: + case I_HEAT_ENERGY_POS: + case I_HEAT_ENERGY_NEG: { + double alpha, rho, Cp; + double height = 1.0, temp = fv->T; + double H_U, dH_U_dtime, H_L, dH_L_dtime; + double dH_U_dX[DIM], dH_L_dX[DIM], dH_U_dp, dH_U_ddh, dH_dF[MDE]; + memset(dH_U_dX, 0, sizeof(dbl) * DIM); + memset(dH_L_dX, 0, sizeof(dbl) * DIM); + memset(dH_dF, 0, sizeof(dbl) * MDE); + + if (ls != NULL && (quantity == I_HEAT_ENERGY_NEG || quantity == I_HEAT_ENERGY_POS)) { + alpha = params == NULL ? ls->Length_Scale : 2.0 * params[0]; + load_lsi(alpha); + } + rho = density(NULL, time); Cp = heat_capacity(NULL, time); @@ -5069,7 +5063,21 @@ int compute_volume_integrand(const int quantity, rho = 1.0; Cp = 1.0; } - *sum += rho * Cp * fv->T * weight * det; + if (is_shell) { + height = height_function_model(&H_U, &dH_U_dtime, &H_L, &dH_L_dtime, dH_U_dX, dH_L_dX, + &dH_U_dp, &dH_U_ddh, dH_dF, time, delta_t); + if (mp->HeightUFunctionModel == WALL_DISTMOD || mp->HeightUFunctionModel == WALL_DISTURB) { + height = MAX(H_U - H_L, DBL_SEMI_SMALL); + } + temp = fv->sh_t; + } + + if (quantity == I_HEAT_ENERGY_NEG) { + height *= lsi->H; + } else if (quantity == I_HEAT_ENERGY_POS) { + height *= (1.0 - lsi->H); + } + *sum += height * rho * Cp * temp * weight * det; if (J_AC != NULL) { for (p = 0; p < dim; p++) { @@ -5079,17 +5087,33 @@ int compute_volume_integrand(const int quantity, for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { J_AC[ei[pg->imtrx]->gun_list[var][j]] += - rho * Cp * fv->T * weight * + height * rho * Cp * temp * weight * (h3 * bf[pd->ShapeVar]->d_det_J_dm[p][j] + fv->dh3dq[p] * bf[var]->phi[j] * det_J); + J_AC[ei[pg->imtrx]->gun_list[var][j]] += + (dH_U_dX[p] - dH_L_dX[p]) * bf[var]->phi[j] * rho * Cp * temp * weight * det; } } } - var = TEMPERATURE; + + if (is_shell) { + var = SHELL_TEMPERATURE; + } else { + var = TEMPERATURE; + } + if (pd->v[pg->imtrx][var]) { + + for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { + J_AC[ei[pg->imtrx]->gun_list[var][j]] += + height * rho * Cp * bf[var]->phi[j] * weight * det; + } + } + + var = FILL; if (pd->v[pg->imtrx][var]) { for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { - J_AC[ei[pg->imtrx]->gun_list[var][j]] += rho * Cp * bf[var]->phi[j] * weight * det; + J_AC[ei[pg->imtrx]->gun_list[var][j]] += dH_dF[j] * rho * Cp * temp * weight * det; } } } @@ -5205,7 +5229,7 @@ int compute_volume_integrand(const int quantity, double alpha, height = 1.0; double H_U, dH_U_dtime, H_L, dH_L_dtime; double dH_U_dX[DIM], dH_L_dX[DIM], dH_U_dp, dH_U_ddh, dH_dF[MDE]; - if (pd->e[pg->imtrx][R_LUBP]) + if (is_shell) height = height_function_model(&H_U, &dH_U_dtime, &H_L, &dH_L_dtime, dH_U_dX, dH_L_dX, &dH_U_dp, &dH_U_ddh, dH_dF, time, delta_t); if (adapt_int_flag) { diff --git a/src/mm_input_mp.c b/src/mm_input_mp.c index 575d5825..e2d6e4e6 100644 --- a/src/mm_input_mp.c +++ b/src/mm_input_mp.c @@ -10319,7 +10319,22 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file) mat_ptr->Lub_gpts[1] = 0.5 * (1. - 0.538469310105683); mat_ptr->Lub_gpts[2] = 0.5; mat_ptr->Lub_gpts[3] = 0.5 * (1. + 0.538469310105683); - mat_ptr->Lub_gpts[0] = 0.5 * (1. + 0.906179845938664); + mat_ptr->Lub_gpts[4] = 0.5 * (1. + 0.906179845938664); + } else if (mat_ptr->LubInt_NGP == 6) { + mat_ptr->Lub_wts[0] = 0.5 * 0.171324492379170; + mat_ptr->Lub_wts[1] = 0.5 * 0.360761573048139; + mat_ptr->Lub_wts[2] = 0.5 * 0.467913934572691; + mat_ptr->Lub_wts[3] = mat_ptr->Lub_wts[2]; + mat_ptr->Lub_wts[4] = mat_ptr->Lub_wts[1]; + mat_ptr->Lub_wts[5] = mat_ptr->Lub_wts[0]; + mat_ptr->Lub_gpts[0] = 0.5 * (1. - 0.932469514203152); + mat_ptr->Lub_gpts[1] = 0.5 * (1. - 0.661209386466265); + mat_ptr->Lub_gpts[2] = 0.5 * (1. - 0.238619186083197); + mat_ptr->Lub_gpts[3] = 0.5 * (1. + 0.238619186083197); + mat_ptr->Lub_gpts[4] = 0.5 * (1. + 0.661209386466265); + mat_ptr->Lub_gpts[5] = 0.5 * (1. + 0.932469514203152); + } else if (mat_ptr->LubInt_NGP > MAX_LUB_NGP) { + GOMA_EH(GOMA_ERROR, "Too many Gauss points -- increase MAX_LUB_NGP!"); } else { GOMA_EH(GOMA_ERROR, "Those integration points not defined yet!"); } @@ -10395,10 +10410,10 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file) SCALAR_INPUT, &NO_SPECIES, es); if (!strcasecmp(model_name, "GRADH") || !strcasecmp(model_name, "gradient")) { mat_ptr->Lub_Curv_NormalModel = TRUE; - SPF(es, "\t(%s = %s)", search_string, "GRADH"); + SPF(es, "%s = %s", search_string, "GRADH"); } else if (!strcasecmp(model_name, "DELTA") || !strcasecmp(model_name, "delta")) { mat_ptr->Lub_Curv_NormalModel = FALSE; - SPF(es, "\t(%s = %s)", search_string, "DELTA"); + SPF(es, "%s = %s", search_string, "DELTA"); } else { mat_ptr->Lub_Curv_NormalModel = TRUE; GOMA_WH(model_read, "Defaulting on Lubrication Curvature Normal"); @@ -10411,10 +10426,10 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file) 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"); + SPF(es, "%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"); + SPF(es, "%s = %s", search_string, "LOGARITHMIC"); } else { mat_ptr->Lub_LS_Interpolation = LINEAR; GOMA_WH(model_read, "Defaulting on Lubrication LS Interpolation"); @@ -10485,10 +10500,10 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file) if (!strcasecmp(model_name, "yes") || !strcasecmp(model_name, "true")) { mat_ptr->Lub_Curv_MassLump = TRUE; - SPF(es, "\t(%s = %s)", search_string, "on"); + SPF(es, "%s = %s", search_string, "on"); } else if (!strcasecmp(model_name, "no") || !strcasecmp(model_name, "false")) { mat_ptr->Lub_Curv_MassLump = FALSE; - SPF(es, "\t(%s = %s)", search_string, "off"); + SPF(es, "%s = %s", search_string, "off"); } else { mat_ptr->Lub_Curv_MassLump = TRUE; SPF(es, "\t(%s = %s)", search_string, "on"); @@ -10502,16 +10517,33 @@ void rd_mp_specs(FILE *imp, char input[], int mn, char *echo_file) if (!strcasecmp(model_name, "yes") || !strcasecmp(model_name, "true")) { mat_ptr->Lub_Curv_Modulation = TRUE; - SPF(es, "\t(%s = %s)", search_string, "on"); + SPF(es, "%s = %s", search_string, "on"); } else if (!strcasecmp(model_name, "no") || !strcasecmp(model_name, "false")) { mat_ptr->Lub_Curv_Modulation = FALSE; - SPF(es, "\t(%s = %s)", search_string, "off"); + SPF(es, "%s = %s", search_string, "off"); } else { mat_ptr->Lub_Curv_Modulation = FALSE; SPF(es, "\t(%s = %s)", search_string, "off"); } ECHO(es, echo_file); + /* Shell Lubrication Curvature Combination */ + strcpy(search_string, "Lubrication Curvature Combination"); + model_read = look_for_mat_prop(imp, "Lubrication Curvature Combination", NULL, NULL, NO_USER, + NULL, model_name, SCALAR_INPUT, &NO_SPECIES, es); + + if (!strcasecmp(model_name, "yes") || !strcasecmp(model_name, "true")) { + mat_ptr->Lub_Curv_Combine = TRUE; + SPF(es, "%s = %s", search_string, "on"); + } else if (!strcasecmp(model_name, "no") || !strcasecmp(model_name, "false")) { + mat_ptr->Lub_Curv_Combine = FALSE; + SPF(es, "%s = %s", search_string, "off"); + } else { + mat_ptr->Lub_Curv_Combine = FALSE; + SPF(es, "\t(%s = %s)", search_string, "off"); + } + ECHO(es, echo_file); + } /* End of Lubrication Curvature Section */ /* diff --git a/src/mm_shell_bc.c b/src/mm_shell_bc.c index 8195506e..97f134b7 100644 --- a/src/mm_shell_bc.c +++ b/src/mm_shell_bc.c @@ -277,6 +277,7 @@ void shell_n_dot_curv_bc(double func[DIM], double bound_normal[DIM], bound_dnormal_dx[DIM][DIM][MDE]; int curv_near; double curvX; + const double penalty = upd->strong_penalty; int eqn = R_SHELL_LUB_CURV; if (ei[pg->imtrx]->ielem_dim == 3) @@ -376,8 +377,8 @@ void shell_n_dot_curv_bc(double func[DIM], /* Derivative of normal vector */ for (jj = 0; jj < DIM; jj++) { - d_LSnormal_dF[ii][jj] = grad_II_phi_i[jj] * LSnormal_maginv; - d_LSnormal_dF[ii][jj] -= gradII_F[jj] * d_LSnormal_mag_dF[ii] * pow(LSnormal_maginv, 2); + d_LSnormal_dF[jj][ii] = grad_II_phi_i[jj] * LSnormal_maginv; + d_LSnormal_dF[jj][ii] -= gradII_F[jj] * d_LSnormal_mag_dF[ii] * pow(LSnormal_maginv, 2); } } @@ -447,10 +448,15 @@ void shell_n_dot_curv_bc(double func[DIM], /* Loop over DOFs (j) */ for (j = 0; j < ei[pg->imtrx]->dof[var]; j++) { - - if (pd->e[pg->imtrx][var] && (ibc_flag != -1)) { - for (ii = 0; ii < pd->Num_Dim; ii++) { - d_func[0][var][j] -= curvX * d_LSnormal_dF[ii][j] * bound_normal[ii]; + if (pd->e[pg->imtrx][var]) { + if (ibc_flag > -1) { + for (ii = 0; ii < pd->Num_Dim; ii++) { + d_func[0][var][j] -= curvX * d_LSnormal_dF[ii][j] * bound_normal[ii]; + } + } else if (ibc_flag == -2) { + for (ii = 0; ii < pd->Num_Dim; ii++) { + d_func[0][var][j] -= curvX * penalty * d_LSnormal_dF[ii][j] * bound_normal[ii]; + } } } } // End of loop over DOFs (j) @@ -479,6 +485,11 @@ void shell_n_dot_curv_bc(double func[DIM], if (curv_near) { if (ibc_flag == -1) { func[0] -= curvX * cos(M_PIE * theta_deg / 180.); + } else if (ibc_flag == -2) { + for (ii = 0; ii < pd->Num_Dim; ii++) { + func[0] -= curvX * penalty * LSnormal[ii] * bound_normal[ii]; + } + func[0] += curvX * penalty * cos(M_PIE * theta_deg / 180.); } else { for (ii = 0; ii < pd->Num_Dim; ii++) { func[0] -= curvX * LSnormal[ii] * bound_normal[ii]; diff --git a/src/mm_shell_util.c b/src/mm_shell_util.c index 422a7e23..c74b23dc 100644 --- a/src/mm_shell_util.c +++ b/src/mm_shell_util.c @@ -3553,7 +3553,8 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co if (EQN == R_LUBP_2) { dmu_df = d_mu->pf[0]; } - if (gn->ConstitutiveEquation == THERMAL || gn->ConstitutiveEquation == TABLE) { + if (pd->v[pg->imtrx][SHELL_TEMPERATURE] && + (gn->ConstitutiveEquation == THERMAL || gn->ConstitutiveEquation == TABLE)) { dmu_dT = mp->d_viscosity[SHELL_TEMPERATURE]; } } @@ -3823,6 +3824,7 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co memset(dHc_U_dX, 0.0, sizeof(double) * DIM); memset(dHc_L_dX, 0.0, sizeof(double) * DIM); memset(D_Hc_DX, 0.0, sizeof(double) * DIM * MDE); + GOMA_WH(GOMA_ERROR, "Lubrication Wall Effect assumes constant capillary height for now..."); } else { H_cap = H; for (i = 0; i < dim; i++) { @@ -3849,9 +3851,12 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co this sign convention is opposite of generally accepted one for curvature i.e., 2H = grad-dot-normal_vector vs. 2H = -grad-dot-normal_vector */ CURV = -(cos(dcaU + atan(slopeU)) + cos(dcaL + atan(-slopeL))) / H_cap; + LubAux->op_curv = CURV; /* Curvature - numerical in planview direction */ - if (pd->e[pg->imtrx][SHELL_LUB_CURV]) { + if (mp->Lub_Curv_Combine && pd->e[pg->imtrx][SHELL_LUB_CURV]) { + CURV = fv->sh_l_curv; + } else if (pd->e[pg->imtrx][SHELL_LUB_CURV]) { CURV += fv->sh_l_curv; } if (pd->e[pg->imtrx][SHELL_LUB_CURV_2]) { @@ -3862,28 +3867,39 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co } /* Sensitivity to height */ - if (!pd->e[pg->imtrx][SHELL_LUB_CURV]) { - D_CURV_DH = (cos(dcaU + atan(slopeU)) + cos(dcaL + atan(-slopeL))) / (H_cap * H_cap); + D_CURV_DH = (cos(dcaU + atan(slopeU)) + cos(dcaL + atan(-slopeL))) / (H_cap * H_cap); - /* Sensitivity to level set F */ - for (i = 0; i < ei[pg->imtrx]->dof[VAR]; i++) { - for (j = 0; j < DIM; j++) { - D_CURV_DF[i] += sin(dcaU + atan(slopeU)) / (H_cap * (1 + slopeU * slopeU)) * - dHc_U_dX[j] * lsi->d_normal_dF[j][i]; - D_CURV_DF[i] += sin(dcaL + atan(-slopeL)) / (H_cap * (1 + slopeL * slopeL)) * - dHc_L_dX[j] * lsi->d_normal_dF[j][i]; - } - D_CURV_DF[i] += D_CURV_DH * dH_dF[i]; + /* Sensitivity to level set F */ + for (i = 0; i < ei[pg->imtrx]->dof[VAR]; i++) { + for (j = 0; j < DIM; j++) { + D_CURV_DF[i] += sin(dcaU + atan(slopeU)) / (H_cap * (1 + slopeU * slopeU)) * dHc_U_dX[j] * + lsi->d_normal_dF[j][i]; + D_CURV_DF[i] += sin(dcaL + atan(-slopeL)) / (H_cap * (1 + slopeL * slopeL)) * + dHc_L_dX[j] * lsi->d_normal_dF[j][i]; } + D_CURV_DF[i] += D_CURV_DH * dH_dF[i]; + } - /* Sensitivity to mesh */ + /* Sensitivity to mesh */ + for (i = 0; i < dim; i++) { + for (j = 0; j < n_dof[MESH_DISPLACEMENT1]; j++) { + D_CURV_DX[i][j] += D_CURV_DH * D_Hc_DX[i][j]; + } + } + if (mp->Lub_Curv_Combine) { + for (i = 0; i < ei[pg->imtrx]->dof[VAR]; i++) { + LubAux->dop_curv_df[i] = D_CURV_DF[i]; + } for (i = 0; i < dim; i++) { for (j = 0; j < n_dof[MESH_DISPLACEMENT1]; j++) { - D_CURV_DX[i][j] += D_CURV_DH * D_Hc_DX[i][j]; + LubAux->dop_curv_dx[i][j] = D_CURV_DX[i][j]; } } + } - /* Sensitivity to shell normal */ + /* Sensitivity to shell normal */ + if ((pd->e[pg->imtrx][R_SHELL_NORMAL1]) && (pd->e[pg->imtrx][R_SHELL_NORMAL2]) && + (pd->e[pg->imtrx][R_SHELL_NORMAL3])) { for (i = 0; i < dim; i++) { for (j = 0; j < ei[pg->imtrx]->dof[SHELL_NORMAL1]; j++) { D_CURV_DNORMAL[i][j] += D_CURV_DH * D_H_DNORMAL[i][j]; @@ -4008,7 +4024,7 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co dbl dq_dH = 0., dv_dH = 0., H_inv = 1. / H; dbl dqmag_dF[MDE], factor, ratio = 0., q_mag2; dbl q[DIM], ev[DIM], pgrad, pg_cmp[DIM], dev_dpg[DIM][DIM]; - dbl v_avg[DIM], dq_dT = 0.; + dbl v_avg[DIM], dq_dT = 0., mu_diss = 0., dmu_diss_dT = 0., dmu_diss_dpgrad = 0.; double DQ_DH[DIM]; double D_Q_DF[DIM][MDE], D_V_DF[DIM][MDE], DGRADP_DF[DIM][MDE], DGRADP_DK = 0.; double DGRADP_DX[DIM][DIM][MDE], DGRADP_DNORMAL[DIM][DIM][MDE]; @@ -4182,8 +4198,8 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co dq_dT *= factor; pre_delP *= factor; vis_w /= factor; - } else if (nonmoving_model && (mp->mp2nd->ViscosityModel == CONSTANT || - mp->mp2nd->ViscosityModel == TIME_RAMP)) { + } else if (mp->mp2nd->ViscosityModel == CONSTANT || + mp->mp2nd->ViscosityModel == TIME_RAMP) { if (mp->Lub_LS_Interpolation == LOGARITHMIC) { if (lsi->near || (fv->F > 0 && mp->mp2nd->viscositymask[1]) || (fv->F < 0 && mp->mp2nd->viscositymask[0])) { @@ -4237,6 +4253,11 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co GOMA_WH(GOMA_ERROR, "mp2nd->ViscosityModel needs to be RATIO or CONSTANT...\n"); } } + if (pd->v[pg->imtrx][SHELL_TEMPERATURE]) { + mu_diss = -q_mag * pgrad; + dmu_diss_dT = -dq_dT * pgrad; + dmu_diss_dpgrad = -q_mag - pgrad * dq_gradp; + } memset(q, 0.0, sizeof(double) * DIM); for (i = 0; i < dim; i++) { q[i] += q_mag * ev[i]; @@ -4282,7 +4303,11 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co } /* End of Viscosity Models **/ /* modulate q (moving wall part) if level-set interface present */ + /* This part needs to be redone -- no q_mag for moving models */ if (pd->v[pg->imtrx][VAR] && movwall_model) { + for (i = 0; i < dim; i++) { + q_mag += q[i] * ev[i]; + } /** This isn't quite right */ if (mp->mp2nd->ViscosityModel == RATIO) { ratio = 1. / mp->mp2nd->viscosity; /* Assuming model = RATIO for now */ q_mag2 = q_mag * ratio; @@ -4351,6 +4376,11 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co GOMA_WH(GOMA_ERROR, "mp2nd->ViscosityModel needs to be RATIO or CONSTANT...\n"); } } + if (pd->v[pg->imtrx][SHELL_TEMPERATURE]) { + mu_diss = -q_mag * pgrad; /* Need to add the drag flow part yet */ + dmu_diss_dT = -dq_dT * pgrad; + dmu_diss_dpgrad = -q_mag - pgrad * dq_gradp; + } memset(q, 0.0, sizeof(double) * DIM); for (i = 0; i < dim; i++) { q[i] += q_mag * ev[i]; @@ -4647,11 +4677,15 @@ void calculate_lub_q_v(const int EQN, double time, double dt, double xi[DIM], co memset(LubAux->dq_dnormal, 0.0, sizeof(double) * DIM * DIM * MDE); LubAux->H = H; + LubAux->H_cap = H_cap; LubAux->dH_dp = D_H_DP; LubAux->dH_ddh = D_H_ddh; LubAux->gradP_mag = 0; LubAux->srate = srate; LubAux->mu_star = vis_w; + LubAux->visc_diss = mu_diss; + LubAux->dvisc_diss_dT = dmu_diss_dT; + LubAux->dvisc_diss_dpgrad = dmu_diss_dpgrad; for (i = 0; i < dim; i++) { LubAux->q[i] = q[i]; LubAux->v_avg[i] = v_avg[i]; @@ -5511,7 +5545,9 @@ void calculate_lub_q_v_old( } /* Curvature - numerical in planview direction */ - if (pd->e[pg->imtrx][SHELL_LUB_CURV]) { + if (mp->Lub_Curv_Combine) { + CURV = fv_old->sh_l_curv; + } else if (pd->e[pg->imtrx][SHELL_LUB_CURV]) { CURV += fv_old->sh_l_curv; } if (pd->e[pg->imtrx][SHELL_LUB_CURV_2]) { diff --git a/src/mm_std_models_shell.c b/src/mm_std_models_shell.c index 0825e350..70f58d88 100644 --- a/src/mm_std_models_shell.c +++ b/src/mm_std_models_shell.c @@ -661,14 +661,37 @@ double height_function_model(double *H_U, } exp_term2 = pow(MAX(exp_term, DBL_SEMI_SMALL), 1. / (2. * powerlaw + 1.)); if ((ls != NULL || pfd != NULL) && Fwall_model) { - H *= (1. - lsi->H) * exp_term2 + lsi->H; - dh_grad = (1. - lsi->H) * alpha / (2. * powerlaw + 1.) * pow(exp_term2, -2. * powerlaw); - for (j = 0; j < ei[pg->imtrx]->dof[FILL]; j++) { - dH_dF[j] = H_orig * lsi->d_H_dF[j] * (1. - exp_term2); + double factor = (mp->mp2nd->viscositymask[1] ? (1.0 - lsi->H) : lsi->H); + if (mp->Lub_LS_Interpolation == LOGARITHMIC) { + if (lsi->near || (fv->F > 0 && mp->mp2nd->viscositymask[1]) || + (fv->F < 0 && mp->mp2nd->viscositymask[0])) { + double H_log = (DOUBLE_NONZERO(exp_term2) ? log(1.0 / exp_term2) : 0.0); + if (fabs(fv->F) <= ls->Length_Scale) { + H = pow(H_orig * exp_term2, factor) * pow(H_orig, 1.0 - factor); + dh_grad = (H / (H_orig * exp_term2)) * factor * H_orig * alpha * + exp(-alpha * rel_dist) / (2. * powerlaw + 1.) * pow(H, -2. * powerlaw); + for (j = 0; j < ei[pg->imtrx]->dof[FILL]; j++) { + dH_dF[j] = H * H_log * lsi->d_H_dF[j]; + } + } + } else { + H *= exp_term2; + dh_grad = H_orig * alpha * exp(-alpha * rel_dist) / (2. * powerlaw + 1.) * + pow(H, -2. * powerlaw); + } + } else { + factor = (mp->mp2nd->viscositymask[1] ? (1.0 - lsi->H) : lsi->H); + H *= factor * exp_term2 + (1.0 - factor); + dh_grad = factor * H_orig * alpha * exp(-alpha * rel_dist) / (2. * powerlaw + 1.) * + pow(H, -2. * powerlaw); + for (j = 0; j < ei[pg->imtrx]->dof[FILL]; j++) { + dH_dF[j] = H_orig * lsi->d_H_dF[j] * (1. - exp_term2); + } } } else { H *= exp_term2; - dh_grad = alpha / (2. * powerlaw + 1.) * pow(exp_term2, -2. * powerlaw); + dh_grad = + H_orig * alpha * exp(-alpha * rel_dist) / (2. * powerlaw + 1.) * pow(H, -2. * powerlaw); } if (mp->HeightUFunctionModel == WALL_DISTMOD) { dH_U_dX[0] = dh_grad * fv->grad_ext_field[mp->heightU_ext_field_index][0]; @@ -694,7 +717,6 @@ double height_function_model(double *H_U, dH_L_dX[0] = dH_L_dX[1] = dH_L_dX[2] = 0.; } } - return (H); }