www.mooseframework.org
PorousFlow1PhaseMD_Gaussian.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
9 
10 template <>
11 InputParameters
13 {
14  InputParameters params = validParams<PorousFlowVariableBase>();
15  params.addRequiredCoupledVar("mass_density",
16  "Variable that represents log(mass-density) of the single phase");
17  params.addRequiredRangeCheckedParam<Real>(
18  "al",
19  "al>0",
20  "For this class, the capillary function is assumed to be saturation = "
21  "exp(-(al*porepressure)^2) for porepressure<0.");
22  params.addRequiredRangeCheckedParam<Real>(
23  "density_P0", "density_P0>0", "The density of the fluid phase at zero porepressure");
24  params.addRequiredRangeCheckedParam<Real>(
25  "bulk_modulus", "bulk_modulus>0", "The constant bulk modulus of the fluid phase");
26  params.addClassDescription("This Material is used for the single-phase situation where "
27  "log(mass-density) is the primary variable. calculates the 1 "
28  "porepressure and the 1 saturation in a 1-phase isothermal situation, "
29  "and derivatives of these with respect to the PorousFlowVariables. A "
30  "gaussian capillary function is assumed");
31  return params;
32 }
33 
35  : PorousFlowVariableBase(parameters),
36 
37  _al(getParam<Real>("al")),
38  _al2(std::pow(_al, 2)),
39  _logdens0(std::log(getParam<Real>("density_P0"))),
40  _bulk(getParam<Real>("bulk_modulus")),
41  _recip_bulk(1.0 / _al / _bulk),
42  _recip_bulk2(std::pow(_recip_bulk, 2)),
43 
44  _md_var(_nodal_material ? coupledNodalValue("mass_density") : coupledValue("mass_density")),
45  _gradmd_qp_var(coupledGradient("mass_density")),
46  _md_varnum(coupled("mass_density")),
47  _pvar(_dictator.isPorousFlowVariable(_md_varnum) ? _dictator.porousFlowVariableNum(_md_varnum)
48  : 0)
49 {
50  if (_num_phases != 1)
51  mooseError("The Dictator proclaims that the number of phases is ",
52  _dictator.numPhases(),
53  " whereas PorousFlow1PhaseMD_Gaussian can only be used for 1-phase simulations. Be "
54  "aware that the Dictator has noted your mistake.");
55 }
56 
57 void
59 {
61  buildPS();
62 }
63 
64 void
66 {
67  // size stuff correctly and prepare the derivative matrices with zeroes
69 
70  buildPS();
71 
72  if (!_nodal_material)
73  {
74  if (_md_var[_qp] >= _logdens0)
75  {
76  (*_gradp_qp)[_qp][0] = _gradmd_qp_var[_qp] * _bulk;
77  (*_grads_qp)[_qp][0] = 0.0;
78  }
79  else
80  {
81  (*_gradp_qp)[_qp][0] =
82  _gradmd_qp_var[_qp] / (_recip_bulk - 2.0 * _al * _porepressure[_qp][0]) / _al;
83  (*_grads_qp)[_qp][0] =
84  -2.0 * _al2 * _porepressure[_qp][0] * _saturation[_qp][0] * (*_gradp_qp)[_qp][0];
85  }
86  }
87 
88  if (_dictator.notPorousFlowVariable(_md_varnum))
89  return;
90 
91  if (_md_var[_qp] >= _logdens0)
92  {
93  // fully saturated at the node or quadpoint
94  _dporepressure_dvar[_qp][0][_pvar] = _bulk;
95  _dsaturation_dvar[_qp][0][_pvar] = 0.0;
96  }
97  else
98  {
99  const Real pp = _porepressure[_qp][0];
100  _dporepressure_dvar[_qp][0][_pvar] = 1.0 / (_recip_bulk - 2.0 * _al * pp) / _al;
101  const Real sat = _saturation[_qp][0];
102  _dsaturation_dvar[_qp][0][_pvar] = -2.0 * _al2 * pp * sat * _dporepressure_dvar[_qp][0][_pvar];
103  }
104 
105  if (!_nodal_material)
106  {
107  if (_md_var[_qp] >= _logdens0)
108  {
109  // fully saturated at the quadpoint
110  (*_dgradp_qp_dgradv)[_qp][0][_pvar] = _bulk;
111  (*_dgradp_qp_dv)[_qp][0][_pvar] = 0.0;
112  (*_dgrads_qp_dgradv)[_qp][0][_pvar] = 0.0;
113  (*_dgrads_qp_dv)[_qp][0][_pvar] = 0.0;
114  }
115  else
116  {
117  const Real pp = _porepressure[_qp][0];
118  (*_dgradp_qp_dgradv)[_qp][0][_pvar] = 1.0 / (_recip_bulk - 2.0 * _al * pp) / _al;
119  (*_dgradp_qp_dv)[_qp][0][_pvar] = _gradmd_qp_var[_qp] * 2.0 * _al *
120  _dporepressure_dvar[_qp][0][_pvar] /
121  std::pow(_recip_bulk - 2.0 * _al * pp, 2.0) / _al;
122  const Real sat = _saturation[_qp][0];
123  (*_dgrads_qp_dgradv)[_qp][0][_pvar] =
124  -2.0 * _al2 * pp * sat * (*_dgradp_qp_dgradv)[_qp][0][_pvar];
125  (*_dgrads_qp_dv)[_qp][0][_pvar] =
126  -2.0 * _al2 * _dporepressure_dvar[_qp][0][_pvar] * sat * (*_gradp_qp)[_qp][0];
127  (*_dgrads_qp_dv)[_qp][0][_pvar] +=
128  -2.0 * _al2 * pp * _dsaturation_dvar[_qp][0][_pvar] * (*_gradp_qp)[_qp][0];
129  (*_dgrads_qp_dv)[_qp][0][_pvar] += -2.0 * _al2 * pp * sat * (*_dgradp_qp_dv)[_qp][0][_pvar];
130  }
131  }
132 }
133 
134 void
136 {
137  if (_md_var[_qp] >= _logdens0)
138  {
139  // full saturation
140  _porepressure[_qp][0] = (_md_var[_qp] - _logdens0) * _bulk;
141  _saturation[_qp][0] = 1.0;
142  }
143  else
144  {
145  // v = logdens0 + p/bulk - (al p)^2
146  // 0 = (v-logdens0) - p/bulk + (al p)^2
147  // 2 al p = (1/al/bulk) +/- sqrt((1/al/bulk)^2 - 4(v-logdens0)) (the "minus" sign is chosen)
148  // s = exp(-(al*p)^2)
149  _porepressure[_qp][0] =
150  (_recip_bulk - std::sqrt(_recip_bulk2 + 4.0 * (_logdens0 - _md_var[_qp]))) / (2.0 * _al);
151  _saturation[_qp][0] = std::exp(-std::pow(_al * _porepressure[_qp][0], 2.0));
152  }
153 }
const unsigned int _md_varnum
Moose variable number of the mass-density.
virtual void computeQpProperties() override
const Real _al
Gaussian parameter: saturation = exp(-(al*p)^2)
const Real _logdens0
fluid density = _dens0*exp(P/_bulk)
const unsigned int _pvar
PorousFlow variable number of the mass-density.
MaterialProperty< std::vector< std::vector< Real > > > & _dporepressure_dvar
d(porepressure)/d(PorousFlow variable)
MaterialProperty< std::vector< Real > > & _saturation
Computed nodal or qp saturation of the phases.
InputParameters validParams< PorousFlow1PhaseMD_Gaussian >()
MaterialProperty< std::vector< std::vector< Real > > > & _dsaturation_dvar
d(saturation)/d(PorousFlow variable)
InputParameters validParams< PorousFlowVariableBase >()
virtual void initQpStatefulProperties() override
PorousFlow1PhaseMD_Gaussian(const InputParameters &parameters)
virtual void initQpStatefulProperties() override
virtual void computeQpProperties() override
const VariableValue & _md_var
Nodal or quadpoint value of mass-density of the fluid phase.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
const unsigned int _num_phases
Number of phases.
const Real _bulk
fluid density = _dens0*exp(P/_bulk)
const VariableGradient & _gradmd_qp_var
Gradient(_mass-density at quadpoints)
Base class for thermophysical variable materials, which assemble materials for primary variables such...
MaterialProperty< std::vector< Real > > & _porepressure
Computed nodal or quadpoint values of porepressure of the phases.