LynxDensityThermal.C 4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
/******************************************************************************/
/*                       LYNX, a MOOSE-based application                      */
/*                                                                            */
/*          Copyright (C) 2017 by Antoine B. Jacquey and Mauro Cacace         */
/*             GFZ Potsdam, German Research Centre for Geosciences            */
/*                                                                            */
/*    This program is free software: you can redistribute it and/or modify    */
/*    it under the terms of the GNU General Public License as published by    */
/*      the Free Software Foundation, either version 3 of the License, or     */
/*                     (at your option) any later version.                    */
/*                                                                            */
/*       This program is distributed in the hope that it will be useful,      */
/*       but WITHOUT ANY WARRANTY; without even the implied warranty of       */
/*        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       */
/*                GNU General Public License for more details.                */
/*                                                                            */
/*      You should have received a copy of the GNU General Public License     */
/*    along with this program. If not, see <http://www.gnu.org/licenses/>     */
/******************************************************************************/

#include "LynxDensityThermal.h"

registerMooseObject("LynxApp", LynxDensityThermal);

template <>
InputParameters
validParams<LynxDensityThermal>()
{
  InputParameters params = validParams<LynxDensityBase>();
  params.addClassDescription(
      "Material calculating densities as a simple linear function of temperature.");
Antoine Jacquey's avatar
Antoine Jacquey committed
32
  params.addParam<std::vector<Real>>("beta_fluid", "The fluid thermal expansion coefficient.");
33 34 35 36 37 38 39 40 41 42 43 44
  params.addParam<std::vector<Real>>("beta_solid", "The solid thermal expansion coefficient.");
  params.addParam<Real>("reference_temperature", 0.0, "The reference temperature.");
  params.addParam<FunctionName>("reference_temperature_fct",
                                "The reference temperature given by a function.");
  return params;
}

LynxDensityThermal::LynxDensityThermal(const InputParameters & parameters)
  : LynxDensityBase(parameters),
    _temperature(coupledValue("temperature")),
    _drho_dtemp(declareProperty<Real>("drho_dtemp")),
    _dinvrho_dtemp(declareProperty<Real>("dinvrho_dtemp")),
Antoine Jacquey's avatar
Antoine Jacquey committed
45 46
    _beta_fluid(isParamValid("beta_fluid") ? getLynxParam<Real>("beta_fluid")
                                           : std::vector<Real>(_n_composition, 0.0)),
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
    _beta_solid(isParamValid("beta_solid") ? getLynxParam<Real>("beta_solid")
                                           : std::vector<Real>(_n_composition, 0.0)),
    _temp_ref(getParam<Real>("reference_temperature")),
    _temp_ref_fct(isParamValid("reference_temperature_fct")
                      ? &getFunction("reference_temperature_fct")
                      : NULL)
{
}

void
LynxDensityThermal::computeQpProperties()
{
  computeQpGravity();

  Real temp_ref = _temp_ref;
  if (_temp_ref_fct)
    temp_ref = _temp_ref_fct->value(_t, _q_point[_qp]);

Antoine Jacquey's avatar
Antoine Jacquey committed
65 66
  _rho_f[_qp] = averageProperty(_fluid_density) *
                (1.0 - averageProperty(_beta_fluid) * (_temperature[_qp] - temp_ref));
67 68 69
  _rho_s[_qp] = averageProperty(_solid_density) *
                (1.0 - averageProperty(_beta_solid) * (_temperature[_qp] - temp_ref));

Antoine Jacquey's avatar
Antoine Jacquey committed
70
  _rho_b[_qp] = _porosity[_qp] * _rho_f[_qp] + (1.0 - _porosity[_qp]) * _rho_s[_qp];
71 72
  _reference_rho_b[_qp] = _rho_b[_qp];

Antoine Jacquey's avatar
Antoine Jacquey committed
73
  Real drho_f_dtemp = -1.0 * averageProperty(_beta_fluid) * averageProperty(_fluid_density);
74 75
  Real drho_s_dtemp = -1.0 * averageProperty(_beta_solid) * averageProperty(_solid_density);

Antoine Jacquey's avatar
Antoine Jacquey committed
76
  _drho_dtemp[_qp] = _porosity[_qp] * drho_f_dtemp + (1.0 - _porosity[_qp]) * drho_s_dtemp;
77 78 79 80

  _dinvrho_dtemp[_qp] =
      (_rho_b[_qp] > 0.0) ? -1.0 / std::pow(_rho_b[_qp], 2.0) * _drho_dtemp[_qp] : 0.0;
}