www.mooseframework.org
HEVPFlowRatePowerLawJ2.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 /****************************************************************/
8 
9 template <>
10 InputParameters
12 {
13  InputParameters params = validParams<HEVPFlowRateUOBase>();
14  params.addParam<Real>(
15  "reference_flow_rate", 0.001, "Reference flow rate for rate dependent flow");
16  params.addParam<Real>("flow_rate_exponent", 10.0, "Power law exponent in flow rate equation");
17  params.addParam<Real>("flow_rate_tol", 1e3, "Tolerance for flow rate");
18  params.addClassDescription(
19  "User object to evaluate power law flow rate and flow direction based on J2");
20 
21  return params;
22 }
23 
24 HEVPFlowRatePowerLawJ2::HEVPFlowRatePowerLawJ2(const InputParameters & parameters)
25  : HEVPFlowRateUOBase(parameters),
26  _ref_flow_rate(getParam<Real>("reference_flow_rate")),
27  _flow_rate_exponent(getParam<Real>("flow_rate_exponent")),
28  _flow_rate_tol(getParam<Real>("flow_rate_tol"))
29 {
30 }
31 
32 bool
33 HEVPFlowRatePowerLawJ2::computeValue(unsigned int qp, Real & val) const
34 {
35  RankTwoTensor pk2_dev = computePK2Deviatoric(_pk2[qp], _ce[qp]);
36  Real eqv_stress = computeEqvStress(pk2_dev, _ce[qp]);
37  val = std::pow(eqv_stress / _strength[qp], _flow_rate_exponent) * _ref_flow_rate;
38 
39  if (val > _flow_rate_tol)
40  {
41 #ifdef DEBUG
42  mooseWarning(
43  "Flow rate greater than ", _flow_rate_tol, " ", val, " ", eqv_stress, " ", _strength[qp]);
44 #endif
45  return false;
46  }
47  return true;
48 }
49 
50 bool
51 HEVPFlowRatePowerLawJ2::computeDirection(unsigned int qp, RankTwoTensor & val) const
52 {
53  RankTwoTensor pk2_dev = computePK2Deviatoric(_pk2[qp], _ce[qp]);
54  Real eqv_stress = computeEqvStress(pk2_dev, _ce[qp]);
55 
56  val.zero();
57  if (eqv_stress > 0.0)
58  val = 1.5 / eqv_stress * _ce[qp] * pk2_dev * _ce[qp];
59 
60  return true;
61 }
62 
63 bool
65  const std::string & coupled_var_name,
66  Real & val) const
67 {
68  val = 0.0;
69 
70  if (_strength_prop_name == coupled_var_name)
71  {
72  RankTwoTensor pk2_dev = computePK2Deviatoric(_pk2[qp], _ce[qp]);
73  Real eqv_stress = computeEqvStress(pk2_dev, _ce[qp]);
75  std::pow(eqv_stress / _strength[qp], _flow_rate_exponent) / _strength[qp];
76  }
77 
78  return true;
79 }
80 
81 bool
83  const std::string & coupled_var_name,
84  RankTwoTensor & val) const
85 {
86  val.zero();
87 
88  if (_pk2_prop_name == coupled_var_name)
89  {
90  RankTwoTensor pk2_dev = computePK2Deviatoric(_pk2[qp], _ce[qp]);
91  Real eqv_stress = computeEqvStress(pk2_dev, _ce[qp]);
92  Real dflowrate_dseqv = _ref_flow_rate * _flow_rate_exponent *
93  std::pow(eqv_stress / _strength[qp], _flow_rate_exponent - 1.0) /
94  _strength[qp];
95 
96  RankTwoTensor tau = pk2_dev * _ce[qp];
97  RankTwoTensor dseqv_dpk2dev;
98 
99  dseqv_dpk2dev.zero();
100  if (eqv_stress > 0.0)
101  dseqv_dpk2dev = 1.5 / eqv_stress * tau * _ce[qp];
102 
103  RankTwoTensor ce_inv = _ce[qp].inverse();
104 
105  RankFourTensor dpk2dev_dpk2;
106  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
107  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
108  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
109  for (unsigned int l = 0; l < LIBMESH_DIM; ++l)
110  {
111  dpk2dev_dpk2(i, j, k, l) = 0.0;
112  if (i == k && j == l)
113  dpk2dev_dpk2(i, j, k, l) = 1.0;
114  dpk2dev_dpk2(i, j, k, l) -= ce_inv(i, j) * _ce[qp](k, l) / 3.0;
115  }
116  val = dflowrate_dseqv * dpk2dev_dpk2.transposeMajor() * dseqv_dpk2dev;
117  }
118  return true;
119 }
120 
121 RankTwoTensor
123  const RankTwoTensor & ce) const
124 {
125  return pk2 - (pk2.doubleContraction(ce) * ce.inverse()) / 3.0;
126 }
127 
128 Real
129 HEVPFlowRatePowerLawJ2::computeEqvStress(const RankTwoTensor & pk2_dev,
130  const RankTwoTensor & ce) const
131 {
132  RankTwoTensor sdev = pk2_dev * ce;
133  Real val = sdev.doubleContraction(sdev.transpose());
134  return std::sqrt(1.5 * val);
135 }
Real computeEqvStress(const RankTwoTensor &, const RankTwoTensor &) const
std::string _pk2_prop_name
InputParameters validParams< HEVPFlowRatePowerLawJ2 >()
RankTwoTensor computePK2Deviatoric(const RankTwoTensor &, const RankTwoTensor &) const
virtual bool computeDerivative(unsigned int, const std::string &, Real &) const
virtual bool computeTensorDerivative(unsigned int, const std::string &, RankTwoTensor &) const
virtual bool computeDirection(unsigned int, RankTwoTensor &) const
This user object is a pure virtual base classs Derived classes computes flow rate, direction and derivatives.
const MaterialProperty< Real > & _strength
const MaterialProperty< RankTwoTensor > & _pk2
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
InputParameters validParams< HEVPFlowRateUOBase >()
virtual bool computeValue(unsigned int, Real &) const
HEVPFlowRatePowerLawJ2(const InputParameters &parameters)
const MaterialProperty< RankTwoTensor > & _ce
std::string _strength_prop_name