www.mooseframework.org
ComputeLinearElasticPFFractureStress.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("Phase-field fracture model energy contribution to fracture for "
15  "elasticity and undamaged stress under compressive strain");
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  _F(declareProperty<Real>(getParam<MaterialPropertyName>("F_name"))),
29  _dFdc(declarePropertyDerivative<Real>(getParam<MaterialPropertyName>("F_name"),
30  getVar("c", 0)->name())),
31  _d2Fdc2(declarePropertyDerivative<Real>(
32  getParam<MaterialPropertyName>("F_name"), getVar("c", 0)->name(), getVar("c", 0)->name())),
33  _d2Fdcdstrain(declareProperty<RankTwoTensor>("d2Fdcdstrain")),
34  _dstress_dc(
35  declarePropertyDerivative<RankTwoTensor>(_base_name + "stress", getVar("c", 0)->name()))
36 {
37 }
38 
39 void
41 {
42  const Real c = _c[_qp];
43 
44  // Compute usual Jacobian mult
46 
47  // Zero out values when c > 1
48  Real cfactor = 1.0;
49  if (c > 1.0)
50  cfactor = 0.0;
51 
52  // Compute eigenvectors and eigenvalues of ustress
53  RankTwoTensor eigvec;
54  std::vector<Real> eigval(LIBMESH_DIM);
55  _mechanical_strain[_qp].symmetricEigenvaluesEigenvectors(eigval, eigvec);
56 
57  RankTwoTensor Ipos, Ineg, eigval_tensor;
58  RankTwoTensor I(RankTwoTensor::initIdentity);
59  eigval_tensor.fillFromInputVector(eigval);
60 
61  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
62  if (eigval[i] < 0.0)
63  Ineg(i, i) = 1.0;
64 
65  Ipos = I - Ineg;
66 
67  // Creation of some critical fourth order tensors
68  RankFourTensor QoQ = eigvec.mixedProductIkJl(eigvec);
69  RankFourTensor QToQT = (eigvec.transpose()).mixedProductIkJl(eigvec.transpose());
70  RankFourTensor IpoIp = Ipos.mixedProductIkJl(I);
71  RankFourTensor InoIn = Ineg.mixedProductIkJl(I);
72 
73  const RankFourTensor Jpos = _Jacobian_mult[_qp] * QoQ * IpoIp * QToQT;
74  const RankFourTensor Jneg = _Jacobian_mult[_qp] * QoQ * InoIn * QToQT;
75  _Jacobian_mult[_qp] = cfactor * Jpos * ((1.0 - c) * (1.0 - c) + _kdamage) + Jneg;
76 
77  RankTwoTensor stress0pos = Jpos * _mechanical_strain[_qp];
78  RankTwoTensor stress0neg = Jneg * _mechanical_strain[_qp];
79 
80  // Damage associated with positive component of stress
81  _stress[_qp] = cfactor * stress0pos * ((1.0 - c) * (1.0 - c) + _kdamage) + stress0neg;
82  // _stress[_qp] = _Jacobian_mult[_qp] * _mechanical_strain[_qp];
83 
84  // Energy with positive principal strains
85  const Real G0_pos = stress0pos.doubleContraction(_mechanical_strain[_qp]) / 2.0;
86  const Real G0_neg = stress0neg.doubleContraction(_mechanical_strain[_qp]) / 2.0;
87 
88  // Elastic free energy density
89  _F[_qp] = cfactor * G0_pos * ((1.0 - c) * (1.0 - c) + _kdamage) + G0_neg;
90 
91  // derivative of elastic free energy density wrt c
92  _dFdc[_qp] = -cfactor * G0_pos * 2.0 * (1.0 - c);
93 
94  // 2nd derivative of elastic free energy density wrt c
95  _d2Fdc2[_qp] = cfactor * G0_pos * 2.0;
96 
97  // 2nd derivative wrt c and strain
98  _d2Fdcdstrain[_qp] = -cfactor * stress0pos * (1.0 - c);
99 
100  // Used in StressDivergencePFFracTensors off-diagonal Jacobian
101  _dstress_dc[_qp] = -cfactor * stress0pos * 2.0 * (1.0 - c);
102 }
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
ComputeStressBase is the base class for stress tensors.
MaterialProperty< RankTwoTensor > & _stress
ComputeLinearElasticPFFractureStress(const InputParameters &parameters)
const MaterialProperty< RankTwoTensor > & _mechanical_strain
InputParameters validParams< ComputeLinearElasticPFFractureStress >()
InputParameters validParams< ComputeStressBase >()
Real _kdamage
Small number to avoid non-positive definiteness at or near complete damage.
const MaterialProperty< RankFourTensor > & _elasticity_tensor