www.mooseframework.org
ComputeLinearElasticPFFractureStress.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 #include "MathUtils.h"
12 
14 
15 template <>
16 InputParameters
18 {
19  InputParameters params = validParams<ComputePFFractureStressBase>();
20  params.addClassDescription("Computes the stress and free energy derivatives for the phase field "
21  "fracture model, with small strain");
22  MooseEnum Decomposition("strain_spectral strain_vol_dev stress_spectral none", "none");
23  params.addParam<MooseEnum>("decomposition_type",
24  Decomposition,
25  "Decomposition approaches. Choices are: " +
26  Decomposition.getRawNames());
27  return params;
28 }
29 
31  const InputParameters & parameters)
32  : ComputePFFractureStressBase(parameters),
33  GuaranteeConsumer(this),
34  _decomposition_type(getParam<MooseEnum>("decomposition_type").getEnum<Decomposition_type>())
35 {
36 }
37 
38 void
40 {
44  mooseError("Decomposition approach of strain_vol_dev and strain_spectral can only be used with "
45  "isotropic elasticity tensor materials, use stress_spectral for anistropic "
46  "elasticity tensor materials");
47 }
48 
49 void
51 {
52  // Isotropic elasticity is assumed and should be enforced
53  const Real lambda = _elasticity_tensor[_qp](0, 0, 1, 1);
54  const Real mu = _elasticity_tensor[_qp](0, 1, 0, 1);
55 
56  // Compute eigenvectors and eigenvalues of mechanical strain and projection tensor
57  RankTwoTensor eigvec;
58  std::vector<Real> eigval(LIBMESH_DIM);
59  RankFourTensor Ppos =
60  _mechanical_strain[_qp].positiveProjectionEigenDecomposition(eigval, eigvec);
61  RankFourTensor I4sym(RankFourTensor::initIdentitySymmetricFour);
62 
63  // Calculate tensors of outerproduct of eigen vectors
64  std::vector<RankTwoTensor> etens(LIBMESH_DIM);
65 
66  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
67  etens[i].vectorOuterProduct(eigvec.column(i), eigvec.column(i));
68 
69  // Separate out positive and negative eigen values
70  std::vector<Real> epos(LIBMESH_DIM), eneg(LIBMESH_DIM);
71  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
72  {
73  epos[i] = (std::abs(eigval[i]) + eigval[i]) / 2.0;
74  eneg[i] = -(std::abs(eigval[i]) - eigval[i]) / 2.0;
75  }
76 
77  // Seprate positive and negative sums of all eigenvalues
78  Real etr = 0.0;
79  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
80  etr += eigval[i];
81 
82  const Real etrpos = (std::abs(etr) + etr) / 2.0;
83  const Real etrneg = -(std::abs(etr) - etr) / 2.0;
84 
85  // Calculate the tensile (postive) and compressive (negative) parts of stress
86  RankTwoTensor stress0pos, stress0neg;
87  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
88  {
89  stress0pos += etens[i] * (lambda * etrpos + 2.0 * mu * epos[i]);
90  stress0neg += etens[i] * (lambda * etrneg + 2.0 * mu * eneg[i]);
91  }
92 
93  // sum squares of epos and eneg
94  Real pval(0.0), nval(0.0);
95  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
96  {
97  pval += epos[i] * epos[i];
98  nval += eneg[i] * eneg[i];
99  }
100 
101  _stress[_qp] = stress0pos * _D[_qp] + stress0neg;
102 
103  // Energy with positive principal strains
104  F_pos = lambda * etrpos * etrpos / 2.0 + mu * pval;
105  F_neg = -lambda * etrneg * etrneg / 2.0 + mu * nval;
106 
107  // 2nd derivative wrt c and strain = 0.0 if we used the previous step's history varible
108  if (_use_current_hist)
109  _d2Fdcdstrain[_qp] = stress0pos * _dDdc[_qp];
110 
111  // Used in StressDivergencePFFracTensors off-diagonal Jacobian
112  _dstress_dc[_qp] = stress0pos * _dDdc[_qp];
113 
114  _Jacobian_mult[_qp] = (I4sym - (1 - _D[_qp]) * Ppos) * _elasticity_tensor[_qp];
115 }
116 
117 void
119 {
120  // Compute Uncracked stress
122 
123  // Create the positive and negative projection tensors
124  RankFourTensor I4sym(RankFourTensor::initIdentitySymmetricFour);
125  std::vector<Real> eigval;
126  RankTwoTensor eigvec;
127  RankFourTensor Ppos = stress.positiveProjectionEigenDecomposition(eigval, eigvec);
128 
129  // Project the positive and negative stresses
130  RankTwoTensor stress0pos = Ppos * stress;
131  RankTwoTensor stress0neg = stress - stress0pos;
132 
133  // Compute the positive and negative elastic energies
134  F_pos = (stress0pos).doubleContraction(_mechanical_strain[_qp]) / 2.0;
135  F_neg = (stress0neg).doubleContraction(_mechanical_strain[_qp]) / 2.0;
136 
137  _stress[_qp] = stress0pos * _D[_qp] + stress0neg;
138 
139  // 2nd derivative wrt c and strain = 0.0 if we used the previous step's history varible
140  if (_use_current_hist)
141  _d2Fdcdstrain[_qp] = stress0pos * _dDdc[_qp];
142 
143  // Used in StressDivergencePFFracTensors off-diagonal Jacobian
144  _dstress_dc[_qp] = stress0pos * _dDdc[_qp];
145 
146  _Jacobian_mult[_qp] = (I4sym - (1 - _D[_qp]) * Ppos) * _elasticity_tensor[_qp];
147 }
148 
149 void
151 {
152  // Isotropic elasticity is assumed and should be enforced
153  const Real lambda = _elasticity_tensor[_qp](0, 0, 1, 1);
154  const Real mu = _elasticity_tensor[_qp](0, 1, 0, 1);
155  const Real k = lambda + 2.0 * mu / LIBMESH_DIM;
156 
157  RankTwoTensor I2(RankTwoTensor::initIdentity);
158  RankFourTensor I2I2 = I2.outerProduct(I2);
159 
160  RankFourTensor Jacobian_pos, Jacobian_neg;
161  RankTwoTensor strain0vol, strain0dev;
162  RankTwoTensor stress0pos, stress0neg;
163  Real strain0tr, strain0tr_neg, strain0tr_pos;
164 
165  strain0dev = _mechanical_strain[_qp].deviatoric();
166  strain0vol = _mechanical_strain[_qp] - strain0dev;
167  strain0tr = _mechanical_strain[_qp].trace();
168  strain0tr_neg = std::min(strain0tr, 0.0);
169  strain0tr_pos = strain0tr - strain0tr_neg;
170  stress0neg = k * strain0tr_neg * I2;
171  stress0pos = _elasticity_tensor[_qp] * _mechanical_strain[_qp] - stress0neg;
172  // Energy with positive principal strains
173  RankTwoTensor strain0dev2 = strain0dev * strain0dev;
174  F_pos = 0.5 * k * strain0tr_pos * strain0tr_pos + mu * strain0dev2.trace();
175  F_neg = 0.5 * k * strain0tr_neg * strain0tr_neg;
176 
177  _stress[_qp] = stress0pos * _D[_qp] + stress0neg;
178 
179  // 2nd derivative wrt c and strain = 0.0 if we used the previous step's history varible
180  if (_use_current_hist)
181  _d2Fdcdstrain[_qp] = stress0pos * _dDdc[_qp];
182 
183  // Used in StressDivergencePFFracTensors off-diagonal Jacobian
184  _dstress_dc[_qp] = stress0pos * _dDdc[_qp];
185 
186  if (strain0tr < 0)
187  Jacobian_neg = k * I2I2;
188  Jacobian_pos = _elasticity_tensor[_qp] - Jacobian_neg;
189  _Jacobian_mult[_qp] = _D[_qp] * Jacobian_pos + Jacobian_neg;
190 }
191 
192 void
194 {
195  Real F_pos, F_neg;
196 
197  switch (_decomposition_type)
198  {
200  computeStrainSpectral(F_pos, F_neg);
201  break;
203  computeStrainVolDev(F_pos, F_neg);
204  break;
206  computeStressSpectral(F_pos, F_neg);
207  break;
208  default:
209  {
210  _stress[_qp] = _D[_qp] * _elasticity_tensor[_qp] * _mechanical_strain[_qp];
211  F_pos = (_stress[_qp]).doubleContraction(_mechanical_strain[_qp]) / 2.0;
212  F_neg = 0.0;
213  if (_use_current_hist)
214  _d2Fdcdstrain[_qp] = _stress[_qp] * _dDdc[_qp];
215 
216  _dstress_dc[_qp] = _stress[_qp] * _dDdc[_qp];
217  _Jacobian_mult[_qp] = _D[_qp] * _elasticity_tensor[_qp];
218  }
219  }
220 
221  // // Assign history variable
222  if (F_pos > _H_old[_qp])
223  _H[_qp] = F_pos;
224  else
225  _H[_qp] = _H_old[_qp];
226 
227  Real hist_variable = _H_old[_qp];
228  if (_use_current_hist)
229  hist_variable = _H[_qp];
230 
231  if (hist_variable < _barrier[_qp])
232  hist_variable = _barrier[_qp];
233 
234  // Elastic free energy density
235  _E[_qp] = hist_variable * _D[_qp] + F_neg;
236  _dEdc[_qp] = hist_variable * _dDdc[_qp];
237  _d2Ed2c[_qp] = hist_variable * _d2Dd2c[_qp];
238 }
MaterialProperty< Real > & _H
History variable that prevents crack healing, declared in this material.
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
const MaterialProperty< Real > & _d2Dd2c
Second-order derivative of degradation w.r.t damage variable.
const MaterialProperty< RankFourTensor > & _elasticity_tensor
Elasticity tensor material property.
MaterialProperty< Real > & _E
Material property for elastic energy.
MaterialProperty< Real > & _d2Ed2c
Second-order derivative of elastic energy w.r.t damage variable.
virtual void computeQpStress() override
Compute the stress and store it in the _stress material property for the current quadrature point...
ComputePFFractureStressBase is the base class for stress in phase field fracture model.
MaterialProperty< RankTwoTensor > & _dstress_dc
Derivative of stress w.r.t damage variable.
MaterialProperty< RankTwoTensor > & _d2Fdcdstrain
Second-order derivative of elastic energy w.r.t damage variable and strain.
InputParameters validParams< ComputePFFractureStressBase >()
MaterialProperty< RankTwoTensor > & _stress
Stress material property.
void computeStrainSpectral(Real &F_pos, Real &F_neg)
Method to split elastic energy based on strain spectral decomposition.
const MaterialProperty< Real > & _barrier
material property for fracture energy barrier
registerMooseObject("TensorMechanicsApp", ComputeLinearElasticPFFractureStress)
void computeStressSpectral(Real &F_pos, Real &F_neg)
Method to split elastic energy based on stress spectral decomposition.
enum ComputeLinearElasticPFFractureStress::Decomposition_type _decomposition_type
const MaterialProperty< Real > & _dDdc
Derivative of degradation function w.r.t damage variable.
MaterialProperty< Real > & _dEdc
Derivative of elastic energy w.r.t damage variable.
Phase-field fracture This class computes the stress and energy contribution for the small strain Line...
ComputeLinearElasticPFFractureStress(const InputParameters &parameters)
bool _use_current_hist
Use current value of history variable.
const MaterialProperty< RankTwoTensor > & _mechanical_strain
Mechanical strain material property.
const MaterialProperty< Real > & _H_old
Old value of history variable.
InputParameters validParams< ComputeLinearElasticPFFractureStress >()
const std::string _elasticity_tensor_name
Name of the elasticity tensor material property.
bool hasGuaranteedMaterialProperty(const MaterialPropertyName &prop, Guarantee guarantee)
const MaterialProperty< Real > & _D
Material property for energetic degradation function.
void computeStrainVolDev(Real &F_pos, Real &F_neg)
Method to split elastic energy based on strain volumetric/deviatoric decomposition.
Add-on class that provides the functionality to check if guarantees for material properties are provi...