Skip to content

Commit

Permalink
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'main' into dc_s
Browse files Browse the repository at this point in the history
wortiz authored Aug 28, 2024
2 parents 3f3d449 + 8d85936 commit 7d42484
Showing 11 changed files with 285 additions and 116 deletions.
7 changes: 7 additions & 0 deletions include/mm_as_structs.h
Original file line number Diff line number Diff line change
@@ -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;
1 change: 1 addition & 0 deletions include/mm_mp_structs.h
Original file line number Diff line number Diff line change
@@ -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;
6 changes: 5 additions & 1 deletion include/mm_post_def.h
Original file line number Diff line number Diff line change
@@ -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);

2 changes: 1 addition & 1 deletion include/rf_fem_const.h
Original file line number Diff line number Diff line change
@@ -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 */
1 change: 1 addition & 0 deletions src/dp_vif.c
Original file line number Diff line number Diff line change
@@ -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);
85 changes: 58 additions & 27 deletions src/mm_fill_shell.c
Original file line number Diff line number Diff line change
@@ -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)];
}
}
Loading

0 comments on commit 7d42484

Please sign in to comment.