www.mooseframework.org
PorousFlowEnergyTimeDerivative.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 
12 #include "MooseVariable.h"
13 
15 
18 {
20  params.addParam<bool>("strain_at_nearest_qp",
21  false,
22  "When calculating nodal porosity that depends on strain, use the strain at "
23  "the nearest quadpoint. This adds a small extra computational burden, and "
24  "is not necessary for simulations involving only linear lagrange elements. "
25  " If you set this to true, you will also want to set the same parameter to "
26  "true for related Kernels and Materials");
27  params.addParam<std::string>(
28  "base_name",
29  "For mechanically-coupled systems, this Kernel will depend on the volumetric strain. "
30  "base_name should almost always be the same base_name as given to the TensorMechanics object "
31  "that computes strain. Supplying a base_name to this Kernel but not defining an associated "
32  "TensorMechanics strain calculator means that this Kernel will not depend on volumetric "
33  "strain. That could be useful when models contain solid mechanics that is not coupled to "
34  "porous flow, for example");
35  params.addRequiredParam<UserObjectName>(
36  "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
37  params.set<bool>("use_displaced_mesh") = false;
38  params.suppressParameter<bool>("use_displaced_mesh");
39  params.addClassDescription("Derivative of heat-energy-density wrt time");
40  return params;
41 }
42 
44  : TimeKernel(parameters),
45  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
46  _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
47  _num_phases(_dictator.numPhases()),
48  _fluid_present(_num_phases > 0),
49  _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
50  _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
51  _has_total_strain(hasMaterialProperty<RankTwoTensor>(_base_name + "total_strain")),
52  _total_strain_old(_has_total_strain
53  ? &getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")
54  : nullptr),
55  _porosity(getMaterialProperty<Real>("PorousFlow_porosity_nodal")),
56  _porosity_old(getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")),
57  _dporosity_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_nodal_dvar")),
58  _dporosity_dgradvar(
59  getMaterialProperty<std::vector<RealGradient>>("dPorousFlow_porosity_nodal_dgradvar")),
60  _nearest_qp(_strain_at_nearest_qp
61  ? &getMaterialProperty<unsigned int>("PorousFlow_nearestqp_nodal")
62  : nullptr),
63  _rock_energy_nodal(getMaterialProperty<Real>("PorousFlow_matrix_internal_energy_nodal")),
64  _rock_energy_nodal_old(getMaterialPropertyOld<Real>("PorousFlow_matrix_internal_energy_nodal")),
65  _drock_energy_nodal_dvar(
66  getMaterialProperty<std::vector<Real>>("dPorousFlow_matrix_internal_energy_nodal_dvar")),
67  _fluid_density(_fluid_present ? &getMaterialProperty<std::vector<Real>>(
68  "PorousFlow_fluid_phase_density_nodal")
69  : nullptr),
70  _fluid_density_old(_fluid_present ? &getMaterialPropertyOld<std::vector<Real>>(
71  "PorousFlow_fluid_phase_density_nodal")
72  : nullptr),
73  _dfluid_density_dvar(_fluid_present ? &getMaterialProperty<std::vector<std::vector<Real>>>(
74  "dPorousFlow_fluid_phase_density_nodal_dvar")
75  : nullptr),
76  _fluid_saturation_nodal(
77  _fluid_present ? &getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")
78  : nullptr),
79  _fluid_saturation_nodal_old(
80  _fluid_present ? &getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")
81  : nullptr),
82  _dfluid_saturation_nodal_dvar(_fluid_present
83  ? &getMaterialProperty<std::vector<std::vector<Real>>>(
84  "dPorousFlow_saturation_nodal_dvar")
85  : nullptr),
86  _energy_nodal(_fluid_present ? &getMaterialProperty<std::vector<Real>>(
87  "PorousFlow_fluid_phase_internal_energy_nodal")
88  : nullptr),
89  _energy_nodal_old(_fluid_present ? &getMaterialPropertyOld<std::vector<Real>>(
90  "PorousFlow_fluid_phase_internal_energy_nodal")
91  : nullptr),
92  _denergy_nodal_dvar(_fluid_present ? &getMaterialProperty<std::vector<std::vector<Real>>>(
93  "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
94  : nullptr)
95 {
96 }
97 
98 Real
100 {
102  Real energy = (1.0 - _porosity[_i]) * _rock_energy_nodal[_i];
103  Real energy_old = (1.0 - _porosity_old[_i]) * _rock_energy_nodal_old[_i];
104 
106  if (_fluid_present)
107  for (unsigned ph = 0; ph < _num_phases; ++ph)
108  {
109  energy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
110  (*_energy_nodal)[_i][ph] * _porosity[_i];
111  energy_old += (*_fluid_density_old)[_i][ph] * (*_fluid_saturation_nodal_old)[_i][ph] *
112  (*_energy_nodal_old)[_i][ph] * _porosity_old[_i];
113  }
114  const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
115 
116  return _test[_i][_qp] * (1.0 + strain) * (energy - energy_old) / _dt;
117 }
118 
119 Real
121 {
122  // If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
123  if (!_var_is_porflow_var)
124  return 0.0;
126 }
127 
128 Real
130 {
131  // If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
133  return 0.0;
135 }
136 
137 Real
139 {
140  const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
141 
142  const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
143 
144  // porosity is dependent on variables that are lumped to the nodes,
145  // but it can depend on the gradient
146  // of variables, which are NOT lumped to the nodes, hence:
147  Real denergy = -_dporosity_dgradvar[_i][pvar] * _grad_phi[_j][_i] * _rock_energy_nodal[_i];
148  for (unsigned ph = 0; ph < _num_phases; ++ph)
149  denergy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
150  (*_energy_nodal)[_i][ph] * _dporosity_dgradvar[_i][pvar] * _grad_phi[_j][nearest_qp];
151 
152  if (_i != _j)
153  return _test[_i][_qp] * (1.0 + strain) * denergy / _dt;
154 
155  // As the fluid energy is lumped to the nodes, only non-zero terms are for _i==_j
156  denergy += -_dporosity_dvar[_i][pvar] * _rock_energy_nodal[_i];
157  denergy += (1.0 - _porosity[_i]) * _drock_energy_nodal_dvar[_i][pvar];
158  for (unsigned ph = 0; ph < _num_phases; ++ph)
159  {
160  denergy += (*_dfluid_density_dvar)[_i][ph][pvar] * (*_fluid_saturation_nodal)[_i][ph] *
161  (*_energy_nodal)[_i][ph] * _porosity[_i];
162  denergy += (*_fluid_density)[_i][ph] * (*_dfluid_saturation_nodal_dvar)[_i][ph][pvar] *
163  (*_energy_nodal)[_i][ph] * _porosity[_i];
164  denergy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
165  (*_denergy_nodal_dvar)[_i][ph][pvar] * _porosity[_i];
166  denergy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
167  (*_energy_nodal)[_i][ph] * _dporosity_dvar[_i][pvar];
168  }
169  return _test[_i][_qp] * (1.0 + strain) * denergy / _dt;
170 }
Real computeQpJac(unsigned int pvar) const
Derivative of residual with respect to PorousFlow variable number pvar This is used by both computeQp...
bool notPorousFlowVariable(unsigned int moose_var_num) const
Returns true if moose_var_num is not a porous flow variabe.
MooseVariable & _var
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
unsigned int number() const
const VariablePhiGradient & _grad_phi
T & set(const std::string &name, bool quiet_mode=false)
const MaterialProperty< std::vector< Real > > *const _fluid_density
Nodal fluid density.
const MaterialProperty< std::vector< Real > > & _drock_energy_nodal_dvar
d(nodal rock energy density)/d(PorousFlow variable)
Real & _dt
const PorousFlowDictator & _dictator
PorousFlowDictator UserObject.
void addRequiredParam(const std::string &name, const std::string &doc_string)
PorousFlowEnergyTimeDerivative(const InputParameters &parameters)
void suppressParameter(const std::string &name)
const bool _fluid_present
Whether _num_phases > 0 (ie. there is a fluid present)
const bool _has_total_strain
Whether there is a Material called _base_name_total_strain.
const MaterialProperty< Real > & _porosity
Porosity at the nodes, but it can depend on grad(variables) which are actually evaluated at the qps...
const VariableTestValue & _test
const MaterialProperty< std::vector< RealGradient > > & _dporosity_dgradvar
d(porosity)/d(grad PorousFlow variable) - remember these derivatives will be wrt grad(vars) at qps ...
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
const MaterialProperty< Real > & _rock_energy_nodal
Nodal rock energy density.
const unsigned int _num_phases
Number of fluid phases.
Kernel = (heat_energy - heat_energy_old)/dt It is lumped to the nodes.
unsigned int _i
const MaterialProperty< Real > & _porosity_old
Old value of porosity.
unsigned int _j
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
void addClassDescription(const std::string &doc_string)
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
The PorousFlow variable number.
registerMooseObject("PorousFlowApp", PorousFlowEnergyTimeDerivative)
const MaterialProperty< Real > & _rock_energy_nodal_old
Old value of nodal rock energy density.
const MaterialProperty< std::vector< Real > > & _dporosity_dvar
d(porosity)/d(PorousFlow variable) - these derivatives will be wrt variables at the nodes ...
static InputParameters validParams()
const bool _var_is_porflow_var
Whether the Variable for this Kernel is a PorousFlow variable according to the Dictator.
void ErrorVector unsigned int
unsigned int _qp
const bool _strain_at_nearest_qp
Whether the porosity uses the volumetric strain at the closest quadpoint.