LinearFVPhaseChangeSource.C

Go to the documentation of this file.

#include "LinearFVPhaseChangeSource.h"
 
#include "MooseEnum.h"
 
#include <algorithm>
 
registerMooseObject("NavierStokesApp", LinearFVPhaseChangeSource);
 
InputParameters
LinearFVPhaseChangeSource::validParams()
{
  InputParameters params = LinearFVElementalKernel::validParams();
 
  params.addClassDescription(
      "Linear FV elemental kernel that adds the apparent heat-capacity "
      "phase-change source term: rho * L * (df/dT) * T_dot, with f the liquid "
      "fraction.");
 
  params.addRequiredParam<MooseFunctorName>("L", "Latent heat.");
  params.addRequiredParam<MooseFunctorName>(NS::density, "The mixture density.");
  params.addRequiredParam<MooseFunctorName>("T_solidus", "The solidus temperature.");
  params.addRequiredParam<MooseFunctorName>("T_liquidus", "The liquidus temperature.");
 
  MooseEnum smoothing("smooth sharp", "smooth");
  params.addParam<MooseEnum>(
      "smoothing",
      smoothing,
      "Shape of the liquid fraction over the mushy interval. 'smooth' (default) uses a "
      "smoothstep cubic function and 'sharp' uses a linear liquid fraction.");
 
  return params;
}
 
LinearFVPhaseChangeSource::LinearFVPhaseChangeSource(const InputParameters & params)
  : LinearFVElementalKernel(params),
    _L(getFunctor<Real>("L")),
    _rho(getFunctor<Real>(NS::density)),
    _T_solidus(getFunctor<Real>("T_solidus")),
    _T_liquidus(getFunctor<Real>("T_liquidus")),
    _smooth(getParam<MooseEnum>("smoothing") == "smooth"),
    _time_integrator(_sys.getTimeIntegrator(_var_num)),
    _factor_history(_time_integrator.numStatesRequired(), 0.0),
    _state_args(_time_integrator.numStatesRequired(), determineState())
{
}
 
Real
LinearFVPhaseChangeSource::computeDfDT(const Real T, const Real T_sol, const Real dT_pc) const
{
  if (_smooth)
  {
    // f(s) = 3 s^2 - 2 s^3  =>  df/dT = 6 s (1 - s) / dT_pc
    const Real s = std::clamp((T - T_sol) / dT_pc, 0.0, 1.0);
    return 6.0 * s * (1.0 - s) / dT_pc;
  }
 
  // Sharp: linear liquid fraction, i.e. a top-hat derivative over the open mushy interval
  if (T > T_sol && T < T_sol + dT_pc)
    return 1.0 / dT_pc;
  return 0.0;
}
 
Real
LinearFVPhaseChangeSource::computeMatrixContribution()
{
  // Element context
  const auto state = determineState();
  const auto elem_arg = makeElemArg(_current_elem_info->elem());
 
  const Real T_sol = _T_solidus(elem_arg, state);
  const Real T_liq = _T_liquidus(elem_arg, state);
  const Real dT_pc = T_liq - T_sol;
 
  // Guard against degenerate or inverted mushy interval
  if (dT_pc <= 0.0)
    return 0.0;
 
  const Real T = _var.getElemValue(*_current_elem_info, state);
  const Real dfdT = computeDfDT(T, T_sol, dT_pc);
 
  // Apparent heat capacity term rho * L * (df/dT) * T_dot
  const Real rhoL = _rho(elem_arg, state) * _L(elem_arg, state);
  const Real Tdot_coeff = _time_integrator.timeDerivativeMatrixContribution(1.0);
 
  return rhoL * dfdT * Tdot_coeff * _current_elem_volume;
}
 
Real
LinearFVPhaseChangeSource::computeRightHandSideContribution()
{
  // Element context
  const auto state = determineState();
  const auto elem_arg = makeElemArg(_current_elem_info->elem());
 
  const Real T_sol = _T_solidus(elem_arg, state);
  const Real T_liq = _T_liquidus(elem_arg, state);
  const Real dT_pc = T_liq - T_sol;
 
  // Guard against degenerate or inverted mushy interval
  if (dT_pc <= 0.0)
    return 0.0;
 
  const Real T = _var.getElemValue(*_current_elem_info, state);
  const Real dfdT = computeDfDT(T, T_sol, dT_pc);
 
  // Apparent heat capacity term rho * L * (df/dT) * T_dot (explicit RHS part)
  const Real rhoL = _rho(elem_arg, state) * _L(elem_arg, state);
  const Real Told_dt = _time_integrator.timeDerivativeRHSContribution(_dof_id, _factor_history);
 
  return rhoL * dfdT * Told_dt * _current_elem_volume;
}
 
void
LinearFVPhaseChangeSource::setCurrentElemInfo(const ElemInfo * elem_info)
{
  LinearFVElementalKernel::setCurrentElemInfo(elem_info);
  for (const auto i : index_range(_factor_history))
    _factor_history[i] = 1.0;
}