www.mooseframework.org
ComputeIsotropicLinearElasticPFFractureStress.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<ComputeStressBase>();
14  params.addClassDescription("Computes the stress and free energy derivatives for the phase field "
15  "fracture model, with linear isotropic elasticity");
16  params.addRequiredCoupledVar("c", "Order parameter for damage");
17  params.addParam<Real>("kdamage", 1e-6, "Stiffness of damaged matrix");
18  params.addParam<MaterialPropertyName>(
19  "F_name", "E_el", "Name of material property storing the elastic energy");
20  return params;
21 }
22 
24  const InputParameters & parameters)
25  : ComputeStressBase(parameters),
26  _c(coupledValue("c")),
27  _kdamage(getParam<Real>("kdamage")),
28  _l(getMaterialProperty<Real>("l")),
29  _gc(getMaterialProperty<Real>("gc_prop")),
30  _F(declareProperty<Real>(getParam<MaterialPropertyName>("F_name"))),
31  _dFdc(declarePropertyDerivative<Real>(getParam<MaterialPropertyName>("F_name"),
32  getVar("c", 0)->name())),
33  _d2Fdc2(declarePropertyDerivative<Real>(
34  getParam<MaterialPropertyName>("F_name"), getVar("c", 0)->name(), getVar("c", 0)->name())),
35  _d2Fdcdstrain(declareProperty<RankTwoTensor>("d2Fdcdstrain")),
36  _dstress_dc(
37  declarePropertyDerivative<RankTwoTensor>(_base_name + "stress", getVar("c", 0)->name())),
38  _hist(declareProperty<Real>("hist")),
39  _hist_old(getMaterialPropertyOld<Real>("hist"))
40 {
41 }
42 
43 void
45 {
46  _hist[_qp] = 0.0;
47 }
48 
49 void
51 {
52  const Real c = _c[_qp];
53 
54  // Compute cfactor to have c > 1 behave the same as c = 1
55  Real cfactor = 0.0;
56  if (c < 1.0)
57  cfactor = 1.0;
58 
59  // Isotropic elasticity is assumed and should be enforced
60  const Real lambda = _elasticity_tensor[_qp](0, 0, 1, 1);
61  const Real mu = _elasticity_tensor[_qp](0, 1, 0, 1);
62 
63  // Compute eigenvectors and eigenvalues of mechanical strain
64  RankTwoTensor eigvec;
65  std::vector<Real> eigval(LIBMESH_DIM);
66  _mechanical_strain[_qp].symmetricEigenvaluesEigenvectors(eigval, eigvec);
67 
68  // Calculate tensors of outerproduct of eigen vectors
69  std::vector<RankTwoTensor> etens(LIBMESH_DIM);
70 
71  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
72  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
73  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
74  etens[i](j, k) = eigvec(j, i) * eigvec(k, i);
75 
76  // Separate out positive and negative eigen values
77  std::vector<Real> epos(LIBMESH_DIM), eneg(LIBMESH_DIM);
78  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
79  {
80  epos[i] = (std::abs(eigval[i]) + eigval[i]) / 2.0;
81  eneg[i] = (std::abs(eigval[i]) - eigval[i]) / 2.0;
82  }
83 
84  // Seprate positive and negative sums of all eigenvalues
85  Real etr = 0.0;
86  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
87  etr += eigval[i];
88 
89  const Real etrpos = (std::abs(etr) + etr) / 2.0;
90  const Real etrneg = (std::abs(etr) - etr) / 2.0;
91 
92  // Calculate the tensile (postive) and compressive (negative) parts of stress
93  RankTwoTensor stress0pos, stress0neg;
94  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
95  {
96  stress0pos += etens[i] * (lambda * etrpos + 2.0 * mu * epos[i]);
97  stress0neg += etens[i] * (lambda * etrneg + 2.0 * mu * eneg[i]);
98  }
99 
100  // sum squares of epos and eneg
101  Real pval(0.0), nval(0.0);
102  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
103  {
104  pval += epos[i] * epos[i];
105  nval += eneg[i] * eneg[i];
106  }
107 
108  // Energy with positive principal strains
109  const Real G0_pos = lambda * etrpos * etrpos / 2.0 + mu * pval;
110  const Real G0_neg = lambda * etrneg * etrneg / 2.0 + mu * nval;
111 
112  // Assign history variable and derivative
113  if (G0_pos > _hist_old[_qp])
114  _hist[_qp] = G0_pos;
115  else
116  _hist[_qp] = _hist_old[_qp];
117 
118  // Damage associated with positive component of stress
119  _stress[_qp] =
120  stress0pos * cfactor * ((1.0 - c) * (1.0 - c) * (1 - _kdamage) + _kdamage) - stress0neg;
121 
122  // Elastic free energy density
123  _F[_qp] = _hist[_qp] * cfactor * ((1.0 - c) * (1.0 - c) * (1 - _kdamage) + _kdamage) - G0_neg +
124  _gc[_qp] / (2 * _l[_qp]) * c * c;
125 
126  // derivative of elastic free energy density wrt c
127  _dFdc[_qp] = -_hist[_qp] * 2.0 * cfactor * (1.0 - c) * (1 - _kdamage) + _gc[_qp] / _l[_qp] * c;
128 
129  // 2nd derivative of elastic free energy density wrt c
130  _d2Fdc2[_qp] = _hist[_qp] * 2.0 * cfactor * (1 - _kdamage) + _gc[_qp] / _l[_qp];
131 
132  // 2nd derivative wrt c and strain. Note that I am ignoring the history variable here, but this
133  // approach gets the fastest convergence
134  _d2Fdcdstrain[_qp] = -stress0pos * cfactor * (1.0 - c) * (1 - _kdamage);
135 
136  // Used in StressDivergencePFFracTensors off-diagonal Jacobian
137  _dstress_dc[_qp] = -stress0pos * 2.0 * cfactor * (1.0 - c) * (1 - _kdamage);
138 
139  // Estimate Jacobian mult, exact when c = 0 and VERY APPROXIMATE when c = 1
141 }
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
ComputeStressBase is the base class for stress tensors.
const VariableValue & _c
Coupled order parameter defining the crack.
const Real _kdamage
Small number to avoid non-positive definiteness at or near complete damage.
MaterialProperty< RankTwoTensor > & _stress
const MaterialProperty< Real > & _gc
Material property defining gc parameter, declared elsewhere.
virtual void initQpStatefulProperties()
Function required to initialize statefull material properties.
InputParameters validParams< ComputeIsotropicLinearElasticPFFractureStress >()
const MaterialProperty< RankTwoTensor > & _mechanical_strain
const MaterialProperty< Real > & _l
Material property defining crack width, declared elsewhere.
InputParameters validParams< ComputeStressBase >()
MaterialProperty< Real > & _hist
History variable that prevents crack healing, declared in this material.
const MaterialProperty< RankFourTensor > & _elasticity_tensor
virtual void computeQpStress()
Primary funciton of this material, computes stress and elastic energy.
MaterialProperty< Real > & _F
Elastic energy and derivatives, declared in this material.
const MaterialProperty< Real > & _hist_old
Old value of history variable.