Line data    Source code

       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             : 
      10             : #include "ComputeCrackedStress.h"
      11             : 
      12             : registerMooseObject("SolidMechanicsApp", ComputeCrackedStress);
      13             : 
      14             : InputParameters
      15           0 : ComputeCrackedStress::validParams()
      16             : {
      17           0 :   InputParameters params = Material::validParams();
      18           0 :   params.addClassDescription("Computes energy and modifies the stress for phase field fracture");
      19           0 :   params.addRequiredCoupledVar("c", "Order parameter for damage");
      20           0 :   params.addParam<Real>("kdamage", 1e-9, "Stiffness of damaged matrix");
      21           0 :   params.addParam<bool>("finite_strain_model", false, "The model is using finite strain");
      22           0 :   params.addParam<bool>(
      23           0 :       "use_current_history_variable", false, "Use the current value of the history variable.");
      24           0 :   params.addParam<MaterialPropertyName>(
      25             :       "F_name", "E_el", "Name of material property storing the elastic energy");
      26           0 :   params.addParam<MaterialPropertyName>(
      27             :       "kappa_name",
      28             :       "kappa_op",
      29             :       "Name of material property being created to store the interfacial parameter kappa");
      30           0 :   params.addParam<MaterialPropertyName>(
      31             :       "mobility_name", "L", "Name of material property being created to store the mobility L");
      32           0 :   params.addParam<std::string>("base_name", "The base name used to save the cracked stress");
      33           0 :   params.addRequiredParam<std::string>("uncracked_base_name",
      34             :                                        "The base name used to calculate the original stress");
      35           0 :   return params;
      36           0 : }
      37             : 
      38           0 : ComputeCrackedStress::ComputeCrackedStress(const InputParameters & parameters)
      39             :   : DerivativeMaterialInterface<Material>(parameters),
      40           0 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      41           0 :     _uncracked_base_name(getParam<std::string>("uncracked_base_name") + "_"),
      42           0 :     _finite_strain_model(getParam<bool>("finite_strain_model")),
      43           0 :     _use_current_hist(getParam<bool>("use_current_history_variable")),
      44           0 :     _strain(
      45           0 :         _finite_strain_model
      46           0 :             ? getMaterialPropertyByName<RankTwoTensor>(_uncracked_base_name + "elastic_strain")
      47           0 :             : getMaterialPropertyByName<RankTwoTensor>(_uncracked_base_name + "mechanical_strain")),
      48           0 :     _uncracked_stress(getMaterialPropertyByName<RankTwoTensor>(_uncracked_base_name + "stress")),
      49           0 :     _uncracked_Jacobian_mult(
      50           0 :         getMaterialPropertyByName<RankFourTensor>(_uncracked_base_name + "Jacobian_mult")),
      51           0 :     _c(coupledValue("c")),
      52           0 :     _gc_prop(getMaterialProperty<Real>("gc_prop")),
      53           0 :     _l(getMaterialProperty<Real>("l")),
      54           0 :     _visco(getMaterialProperty<Real>("visco")),
      55           0 :     _kdamage(getParam<Real>("kdamage")),
      56           0 :     _stress(declareProperty<RankTwoTensor>(_base_name + "stress")),
      57           0 :     _F(declareProperty<Real>(getParam<MaterialPropertyName>("F_name"))),
      58           0 :     _dFdc(declarePropertyDerivative<Real>(getParam<MaterialPropertyName>("F_name"),
      59           0 :                                           coupledName("c", 0))),
      60           0 :     _d2Fdc2(declarePropertyDerivative<Real>(
      61           0 :         getParam<MaterialPropertyName>("F_name"), coupledName("c", 0), coupledName("c", 0))),
      62           0 :     _d2Fdcdstrain(declareProperty<RankTwoTensor>("d2Fdcdstrain")),
      63           0 :     _dstress_dc(declarePropertyDerivative<RankTwoTensor>("stress", coupledName("c", 0))),
      64           0 :     _hist(declareProperty<Real>("hist")),
      65           0 :     _hist_old(getMaterialPropertyOld<Real>("hist")),
      66           0 :     _Jacobian_mult(declareProperty<RankFourTensor>(_base_name + "Jacobian_mult")),
      67           0 :     _kappa(declareProperty<Real>(getParam<MaterialPropertyName>("kappa_name"))),
      68           0 :     _L(declareProperty<Real>(getParam<MaterialPropertyName>("mobility_name")))
      69             : {
      70           0 : }
      71             : 
      72             : void
      73           0 : ComputeCrackedStress::initQpStatefulProperties()
      74             : {
      75           0 :   _stress[_qp].zero();
      76           0 :   _hist[_qp] = 0.0;
      77           0 : }
      78             : 
      79             : void
      80           0 : ComputeCrackedStress::computeQpProperties()
      81             : {
      82           0 :   const Real c = _c[_qp];
      83             : 
      84             :   // Zero out values when c > 1
      85             :   Real cfactor = 1.0;
      86           0 :   if (c > 1.0)
      87             :     cfactor = 0.0;
      88             : 
      89             :   // Create the positive and negative projection tensors
      90           0 :   RankFourTensor I4sym(RankFourTensor::initIdentitySymmetricFour);
      91             :   std::vector<Real> eigval;
      92           0 :   RankTwoTensor eigvec;
      93           0 :   RankFourTensor Ppos = _uncracked_stress[_qp].positiveProjectionEigenDecomposition(eigval, eigvec);
      94           0 :   RankFourTensor Pneg = I4sym - Ppos;
      95             : 
      96             :   // Project the positive and negative stresses
      97           0 :   RankTwoTensor stress0pos = Ppos * _uncracked_stress[_qp];
      98           0 :   RankTwoTensor stress0neg = Pneg * _uncracked_stress[_qp];
      99             : 
     100             :   // Compute the positive and negative elastic energies
     101           0 :   Real G0_pos = (stress0pos).doubleContraction(_strain[_qp]) / 2.0;
     102           0 :   Real G0_neg = (stress0neg).doubleContraction(_strain[_qp]) / 2.0;
     103             : 
     104             :   // Update the history variable
     105           0 :   if (G0_pos > _hist_old[_qp])
     106           0 :     _hist[_qp] = G0_pos;
     107             :   else
     108           0 :     _hist[_qp] = _hist_old[_qp];
     109             : 
     110           0 :   Real hist_variable = _hist_old[_qp];
     111           0 :   if (_use_current_hist)
     112           0 :     hist_variable = _hist[_qp];
     113             : 
     114             :   // Compute degredation function and derivatives
     115           0 :   Real h = cfactor * (1.0 - c) * (1.0 - c) * (1.0 - _kdamage) + _kdamage;
     116           0 :   Real dhdc = -2.0 * cfactor * (1.0 - c) * (1.0 - _kdamage);
     117           0 :   Real d2hdc2 = 2.0 * cfactor * (1.0 - _kdamage);
     118             : 
     119             :   // Compute stress and its derivatives
     120           0 :   _stress[_qp] = (Ppos * h + Pneg) * _uncracked_stress[_qp];
     121           0 :   _dstress_dc[_qp] = stress0pos * dhdc;
     122           0 :   _Jacobian_mult[_qp] = (Ppos * h + Pneg) * _uncracked_Jacobian_mult[_qp];
     123             : 
     124             :   // Compute energy and its derivatives
     125           0 :   _F[_qp] = hist_variable * h - G0_neg + _gc_prop[_qp] * c * c / (2 * _l[_qp]);
     126           0 :   _dFdc[_qp] = hist_variable * dhdc + _gc_prop[_qp] * c / _l[_qp];
     127           0 :   _d2Fdc2[_qp] = hist_variable * d2hdc2 + _gc_prop[_qp] / _l[_qp];
     128             : 
     129             :   // 2nd derivative wrt c and strain = 0.0 if we used the previous step's history varible
     130           0 :   if (_use_current_hist)
     131           0 :     _d2Fdcdstrain[_qp] = stress0pos * dhdc;
     132             : 
     133             :   // Assign L and kappa
     134           0 :   _kappa[_qp] = _gc_prop[_qp] * _l[_qp];
     135           0 :   _L[_qp] = 1.0 / (_gc_prop[_qp] * _visco[_qp]);
     136           0 : }