https://mooseframework.inl.gov
HyperElasticPhaseFieldIsoDamage.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
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 #include "libmesh/utility.h"
12 
14 
17 {
19  params.addParam<bool>("numerical_stiffness", false, "Flag for numerical stiffness");
20  params.addParam<Real>("damage_stiffness", 1e-8, "Avoid zero after complete damage");
21  params.addParam<Real>("zero_tol", 1e-12, "Tolerance for numerical zero");
22  params.addParam<Real>(
23  "zero_perturb", 1e-8, "Perturbation value when strain value less than numerical zero");
24  params.addParam<Real>("perturbation_scale_factor", 1e-5, "Perturbation scale factor");
25  params.addRequiredCoupledVar("c", "Damage variable");
26  params.addParam<bool>(
27  "use_current_history_variable", false, "Use the current value of the history variable.");
28  params.addParam<MaterialPropertyName>(
29  "F_name", "E_el", "Name of material property storing the elastic energy");
30  params.addClassDescription(
31  "Computes damaged stress and energy in the intermediate configuration assuming isotropy");
32 
33  return params;
34 }
35 
38  _num_stiffness(getParam<bool>("numerical_stiffness")),
39  _kdamage(getParam<Real>("damage_stiffness")),
40  _use_current_hist(getParam<bool>("use_current_history_variable")),
41  _l(getMaterialProperty<Real>("l")),
42  _gc(getMaterialProperty<Real>("gc_prop")),
43  _zero_tol(getParam<Real>("zero_tol")),
44  _zero_pert(getParam<Real>("zero_perturb")),
45  _pert_val(getParam<Real>("perturbation_scale_factor")),
46  _c(coupledValue("c")),
47  _save_state(false),
48  _dstress_dc(
49  declarePropertyDerivative<RankTwoTensor>(_base_name + "stress", coupledName("c", 0))),
50  _etens(LIBMESH_DIM),
51  _F(declareProperty<Real>(getParam<MaterialPropertyName>("F_name"))),
52  _dFdc(declarePropertyDerivative<Real>(getParam<MaterialPropertyName>("F_name"),
53  coupledName("c", 0))),
54  _d2Fdc2(declarePropertyDerivative<Real>(
55  getParam<MaterialPropertyName>("F_name"), coupledName("c", 0), coupledName("c", 0))),
56  _d2Fdcdstrain(declareProperty<RankTwoTensor>("d2Fdcdstrain")),
57  _hist(declareProperty<Real>("hist")),
58  _hist_old(getMaterialPropertyOld<Real>("hist"))
59 {
60 }
61 
62 void
64 {
66 
67  _save_state = true;
69  _pk2[_qp] = _pk2_tmp;
70 
71  _save_state = false;
72  if (_num_stiffness)
74 
75  if (_num_stiffness)
77 
78  _dce_dfe.zero();
79  for (const auto i : make_range(Moose::dim))
80  for (const auto j : make_range(Moose::dim))
81  for (const auto k : make_range(Moose::dim))
82  {
83  _dce_dfe(i, j, k, i) = _dce_dfe(i, j, k, i) + _fe(k, j);
84  _dce_dfe(i, j, k, j) = _dce_dfe(i, j, k, j) + _fe(k, i);
85  }
86 
88 }
89 
90 void
92 {
93  Real lambda = _elasticity_tensor[_qp](0, 0, 1, 1);
94  Real mu = _elasticity_tensor[_qp](0, 1, 0, 1);
95 
96  Real c = _c[_qp];
97  Real xfac = Utility::pow<2>(1.0 - c) + _kdamage;
98 
99  std::vector<Real> eigval;
100  RankTwoTensor evec;
102 
103  for (const auto i : make_range(Moose::dim))
105 
106  Real etr = 0.0;
107  for (const auto i : make_range(Moose::dim))
108  etr += eigval[i];
109 
110  Real etrpos = (std::abs(etr) + etr) / 2.0;
111  Real etrneg = (std::abs(etr) - etr) / 2.0;
112 
113  RankTwoTensor pk2pos, pk2neg;
114 
115  for (const auto i : make_range(Moose::dim))
116  {
117  pk2pos += _etens[i] * (lambda * etrpos + 2.0 * mu * (std::abs(eigval[i]) + eigval[i]) / 2.0);
118  pk2neg += _etens[i] * (lambda * etrneg + 2.0 * mu * (std::abs(eigval[i]) - eigval[i]) / 2.0);
119  }
120 
121  _pk2_tmp = pk2pos * xfac - pk2neg;
122 
123  if (_save_state)
124  {
125  std::vector<Real> epos(LIBMESH_DIM), eneg(LIBMESH_DIM);
126  for (const auto i : make_range(Moose::dim))
127  {
128  epos[i] = (std::abs(eigval[i]) + eigval[i]) / 2.0;
129  eneg[i] = (std::abs(eigval[i]) - eigval[i]) / 2.0;
130  }
131 
132  // sum squares of epos and eneg
133  Real pval(0.0), nval(0.0);
134  for (const auto i : make_range(Moose::dim))
135  {
136  pval += epos[i] * epos[i];
137  nval += eneg[i] * eneg[i];
138  }
139 
140  // Energy with positive principal strains
141  const Real G0_pos = lambda * etrpos * etrpos / 2.0 + mu * pval;
142  const Real G0_neg = lambda * etrneg * etrneg / 2.0 + mu * nval;
143 
144  // Assign history variable and derivative
145  if (G0_pos > _hist_old[_qp])
146  _hist[_qp] = G0_pos;
147  else
148  _hist[_qp] = _hist_old[_qp];
149 
150  Real hist_variable = _hist_old[_qp];
151  if (_use_current_hist)
152  hist_variable = _hist[_qp];
153 
154  // Elastic free energy density
155  _F[_qp] = hist_variable * xfac - G0_neg + _gc[_qp] / (2 * _l[_qp]) * c * c;
156 
157  // derivative of elastic free energy density wrt c
158  _dFdc[_qp] = -hist_variable * 2.0 * (1.0 - c) * (1 - _kdamage) + _gc[_qp] / _l[_qp] * c;
159 
160  // 2nd derivative of elastic free energy density wrt c
161  _d2Fdc2[_qp] = hist_variable * 2.0 * (1 - _kdamage) + _gc[_qp] / _l[_qp];
162 
163  _dG0_dee = pk2pos;
164 
165  _dpk2_dc = -pk2pos * 2.0 * (1.0 - c);
166  }
167 }
168 
169 void
171 {
172  RankTwoTensor ee_tmp;
173 
174  for (const auto i : make_range(Moose::dim))
175  for (unsigned int j = i; j < LIBMESH_DIM; ++j)
176  {
177  Real ee_pert = _zero_pert;
178  if (std::abs(_ee(i, j)) > _zero_tol)
179  ee_pert = _pert_val * std::abs(_ee(i, j));
180 
181  ee_tmp = _ee;
182  _ee(i, j) += ee_pert;
183 
185 
186  for (const auto k : make_range(Moose::dim))
187  for (const auto l : make_range(Moose::dim))
188  {
189  _dpk2_dee(k, l, i, j) = (_pk2_tmp(k, l) - _pk2[_qp](k, l)) / ee_pert;
190  _dpk2_dee(k, l, j, i) = (_pk2_tmp(k, l) - _pk2[_qp](k, l)) / ee_pert;
191  }
192  _ee = ee_tmp;
193  }
194 }
195 
196 void
198 {
200 
202  RankTwoTensor dG0_dfe = _dce_dfe.innerProductTranspose(dG0_dce);
203  RankTwoTensor dG0_df = _dfe_df.innerProductTranspose(dG0_dfe);
204 
205  // 2nd derivative wrt c and strain = 0.0 if we used the previous step's history varible
206  if (_use_current_hist)
207  _d2Fdcdstrain[_qp] =
208  -_df_dstretch_inc.innerProductTranspose(dG0_df) * 2.0 * (1.0 - _c[_qp]) * (1 - _kdamage);
209 
210  usingTensorIndices(i_, j_, k_, l_);
211  _dstress_dc[_qp] = _fe.times<i_, k_, j_, l_>(_fe) * _dpk2_dc;
212 }
This class solves the viscoplastic flow rate equations in the total form Involves 4 different types o...
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
virtual void computeElasticStrain()
Computes elastic Lagrangian strain.
virtual void computePK2StressAndDerivative()
This function computes PK2 stress.
RankTwoTensorTempl< T > innerProductTranspose(const RankTwoTensorTempl< T > &) const
const MaterialProperty< Real > & _hist_old
Old value of history variable.
bool _use_current_hist
Use current value of history variable.
virtual void computeDamageStress()
This function computes PK2 stress modified to account for damage Computes numerical stiffness if flag...
static constexpr std::size_t dim
HyperElasticPhaseFieldIsoDamage(const InputParameters &parameters)
virtual void computeQpJacobian()
This function computes the Jacobian.
Real _zero_pert
Perturbation value for near zero or zero strain components.
static RankTwoTensorTempl< Real > selfOuterProduct(const libMesh::TypeVector< Real > &)
virtual void computeNumStiffness()
This function computes numerical stiffness.
static const std::string mu
Definition: NS.h:123
libMesh::VectorValue< Real > column(const unsigned int i) const
bool _save_state
Flag to save couple material properties.
const VariableValue & _c
Compupled damage variable.
Real _pert_val
Perturbation value for strain components.
Real _zero_tol
Used in numerical stiffness calculation to check near zero values.
const MaterialProperty< Real > & _l
Material property defining crack width, declared elsewhere.
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
This class solves visco plastic model based on isotropically damaged stress The damage parameter is o...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
MaterialProperty< RankTwoTensor > & _dstress_dc
void symmetricEigenvaluesEigenvectors(std::vector< Real > &eigvals, RankTwoTensorTempl< Real > &eigvecs) const
RankFourTensorTempl< Real > times(const RankTwoTensorTempl< Real > &b) const
registerMooseObject("SolidMechanicsApp", HyperElasticPhaseFieldIsoDamage)
MaterialProperty< RankTwoTensor > & _d2Fdcdstrain
IntRange< T > make_range(T beg, T end)
void addClassDescription(const std::string &doc_string)
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
bool _num_stiffness
Flag to compute numerical stiffness.
MaterialProperty< Real > & _F
Elastic energy and derivatives, declared in this material.
virtual void computeQpJacobian()
This function computes tensors used to construct diagonal and off-diagonal Jacobian.
Real _kdamage
Small stiffness of completely damaged material point.
const MaterialProperty< RankFourTensor > & _elasticity_tensor
Elasticity tensor material property.
static const std::string k
Definition: NS.h:130
MaterialProperty< Real > & _hist
History variable that prevents crack healing, declared in this material.
const MaterialProperty< Real > & _gc
Material property defining gc parameter, declared elsewhere.