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 "ComputeLinearElasticPFFractureStress.h"
      11             : #include "MathUtils.h"
      12             : 
      13             : registerMooseObject("SolidMechanicsApp", ComputeLinearElasticPFFractureStress);
      14             : 
      15             : InputParameters
      16           0 : ComputeLinearElasticPFFractureStress::validParams()
      17             : {
      18           0 :   InputParameters params = ComputePFFractureStressBase::validParams();
      19           0 :   params.addClassDescription("Computes the stress and free energy derivatives for the phase field "
      20             :                              "fracture model, with small strain");
      21           0 :   MooseEnum Decomposition("strain_spectral strain_vol_dev stress_spectral none", "none");
      22           0 :   params.addParam<MooseEnum>("decomposition_type",
      23             :                              Decomposition,
      24           0 :                              "Decomposition approaches.  Choices are: " +
      25           0 :                                  Decomposition.getRawNames());
      26           0 :   return params;
      27           0 : }
      28             : 
      29           0 : ComputeLinearElasticPFFractureStress::ComputeLinearElasticPFFractureStress(
      30           0 :     const InputParameters & parameters)
      31             :   : ComputePFFractureStressBase(parameters),
      32             :     GuaranteeConsumer(this),
      33           0 :     _decomposition_type(getParam<MooseEnum>("decomposition_type").getEnum<Decomposition_type>())
      34             : {
      35           0 : }
      36             : 
      37             : void
      38           0 : ComputeLinearElasticPFFractureStress::initialSetup()
      39             : {
      40           0 :   if ((_decomposition_type == Decomposition_type::strain_vol_dev ||
      41           0 :        _decomposition_type == Decomposition_type::strain_spectral) &&
      42           0 :       !hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC))
      43           0 :     mooseError("Decomposition approach of strain_vol_dev and strain_spectral can only be used with "
      44             :                "isotropic elasticity tensor materials, use stress_spectral for anistropic "
      45             :                "elasticity tensor materials");
      46           0 : }
      47             : 
      48             : void
      49           0 : ComputeLinearElasticPFFractureStress::computeStrainSpectral(Real & F_pos, Real & F_neg)
      50             : {
      51             :   // Isotropic elasticity is assumed and should be enforced
      52           0 :   const Real lambda = _elasticity_tensor[_qp](0, 0, 1, 1);
      53           0 :   const Real mu = _elasticity_tensor[_qp](0, 1, 0, 1);
      54             : 
      55           0 :   RankTwoTensor I2(RankTwoTensor::initIdentity);
      56             : 
      57             :   // Compute eigenvectors and eigenvalues of mechanical strain and projection tensor
      58           0 :   RankTwoTensor eigvec;
      59           0 :   std::vector<Real> eigval(LIBMESH_DIM);
      60             :   RankFourTensor Ppos =
      61           0 :       _mechanical_strain[_qp].positiveProjectionEigenDecomposition(eigval, eigvec);
      62           0 :   RankFourTensor I4sym(RankFourTensor::initIdentitySymmetricFour);
      63             : 
      64             :   // Calculate tensors of outerproduct of eigen vectors
      65           0 :   std::vector<RankTwoTensor> etens(LIBMESH_DIM);
      66             : 
      67           0 :   for (const auto i : make_range(Moose::dim))
      68           0 :     etens[i] = RankTwoTensor::selfOuterProduct(eigvec.column(i));
      69             : 
      70             :   // Separate out positive and negative eigen values
      71           0 :   std::vector<Real> epos(LIBMESH_DIM), eneg(LIBMESH_DIM);
      72           0 :   for (const auto i : make_range(Moose::dim))
      73             :   {
      74           0 :     epos[i] = (std::abs(eigval[i]) + eigval[i]) / 2.0;
      75           0 :     eneg[i] = -(std::abs(eigval[i]) - eigval[i]) / 2.0;
      76             :   }
      77             : 
      78             :   // Seprate positive and negative sums of all eigenvalues
      79             :   Real etr = 0.0;
      80           0 :   for (const auto i : make_range(Moose::dim))
      81           0 :     etr += eigval[i];
      82             : 
      83           0 :   const Real etrpos = (std::abs(etr) + etr) / 2.0;
      84           0 :   const Real etrneg = -(std::abs(etr) - etr) / 2.0;
      85             : 
      86             :   // Calculate the tensile (postive) and compressive (negative) parts of stress
      87           0 :   RankTwoTensor stress0pos, stress0neg;
      88           0 :   for (const auto i : make_range(Moose::dim))
      89             :   {
      90           0 :     stress0pos += etens[i] * (lambda * etrpos + 2.0 * mu * epos[i]);
      91           0 :     stress0neg += etens[i] * (lambda * etrneg + 2.0 * mu * eneg[i]);
      92             :   }
      93             : 
      94             :   // sum squares of epos and eneg
      95             :   Real pval(0.0), nval(0.0);
      96           0 :   for (const auto i : make_range(Moose::dim))
      97             :   {
      98           0 :     pval += epos[i] * epos[i];
      99           0 :     nval += eneg[i] * eneg[i];
     100             :   }
     101             : 
     102           0 :   _stress[_qp] = stress0pos * _D[_qp] - _pressure[_qp] * I2 * _I[_qp] + stress0neg;
     103             : 
     104             :   // Energy with positive principal strains
     105           0 :   F_pos = lambda * etrpos * etrpos / 2.0 + mu * pval;
     106           0 :   F_neg = -lambda * etrneg * etrneg / 2.0 + mu * nval;
     107             : 
     108             :   // 2nd derivative wrt c and strain = 0.0 if we used the previous step's history varible
     109           0 :   if (_use_current_hist)
     110           0 :     _d2Fdcdstrain[_qp] = stress0pos * _dDdc[_qp];
     111             : 
     112             :   // Used in StressDivergencePFFracTensors off-diagonal Jacobian
     113           0 :   _dstress_dc[_qp] = stress0pos * _dDdc[_qp] - _pressure[_qp] * I2 * _dIdc[_qp];
     114             : 
     115           0 :   _Jacobian_mult[_qp] = (I4sym - (1 - _D[_qp]) * Ppos) * _elasticity_tensor[_qp];
     116           0 : }
     117             : 
     118             : void
     119           0 : ComputeLinearElasticPFFractureStress::computeStressSpectral(Real & F_pos, Real & F_neg)
     120             : {
     121             :   // Compute Uncracked stress
     122           0 :   RankTwoTensor stress = _elasticity_tensor[_qp] * _mechanical_strain[_qp];
     123             : 
     124           0 :   RankTwoTensor I2(RankTwoTensor::initIdentity);
     125             : 
     126             :   // Create the positive and negative projection tensors
     127           0 :   RankFourTensor I4sym(RankFourTensor::initIdentitySymmetricFour);
     128             :   std::vector<Real> eigval;
     129           0 :   RankTwoTensor eigvec;
     130           0 :   RankFourTensor Ppos = stress.positiveProjectionEigenDecomposition(eigval, eigvec);
     131             : 
     132             :   // Project the positive and negative stresses
     133           0 :   RankTwoTensor stress0pos = Ppos * stress;
     134           0 :   RankTwoTensor stress0neg = stress - stress0pos;
     135             : 
     136             :   // Compute the positive and negative elastic energies
     137           0 :   F_pos = (stress0pos).doubleContraction(_mechanical_strain[_qp]) / 2.0;
     138           0 :   F_neg = (stress0neg).doubleContraction(_mechanical_strain[_qp]) / 2.0;
     139             : 
     140           0 :   _stress[_qp] = stress0pos * _D[_qp] - _pressure[_qp] * I2 * _I[_qp] + stress0neg;
     141             : 
     142             :   // 2nd derivative wrt c and strain = 0.0 if we used the previous step's history varible
     143           0 :   if (_use_current_hist)
     144           0 :     _d2Fdcdstrain[_qp] = stress0pos * _dDdc[_qp];
     145             : 
     146             :   // Used in StressDivergencePFFracTensors off-diagonal Jacobian
     147           0 :   _dstress_dc[_qp] = stress0pos * _dDdc[_qp] - _pressure[_qp] * I2 * _dIdc[_qp];
     148             : 
     149           0 :   _Jacobian_mult[_qp] = (I4sym - (1 - _D[_qp]) * Ppos) * _elasticity_tensor[_qp];
     150           0 : }
     151             : 
     152             : void
     153           0 : ComputeLinearElasticPFFractureStress::computeStrainVolDev(Real & F_pos, Real & F_neg)
     154             : {
     155             :   // Isotropic elasticity is assumed and should be enforced
     156           0 :   const Real lambda = _elasticity_tensor[_qp](0, 0, 1, 1);
     157           0 :   const Real mu = _elasticity_tensor[_qp](0, 1, 0, 1);
     158           0 :   const Real k = lambda + 2.0 * mu / LIBMESH_DIM;
     159             : 
     160           0 :   RankTwoTensor I2(RankTwoTensor::initIdentity);
     161             :   RankFourTensor I2I2 = I2.outerProduct(I2);
     162             : 
     163           0 :   RankFourTensor Jacobian_pos, Jacobian_neg;
     164           0 :   RankTwoTensor strain0vol, strain0dev;
     165           0 :   RankTwoTensor stress0pos, stress0neg;
     166             :   Real strain0tr, strain0tr_neg, strain0tr_pos;
     167             : 
     168           0 :   strain0dev = _mechanical_strain[_qp].deviatoric();
     169           0 :   strain0vol = _mechanical_strain[_qp] - strain0dev;
     170           0 :   strain0tr = _mechanical_strain[_qp].trace();
     171           0 :   strain0tr_neg = std::min(strain0tr, 0.0);
     172           0 :   strain0tr_pos = strain0tr - strain0tr_neg;
     173           0 :   stress0neg = k * strain0tr_neg * I2;
     174           0 :   stress0pos = _elasticity_tensor[_qp] * _mechanical_strain[_qp] - stress0neg;
     175             :   // Energy with positive principal strains
     176           0 :   RankTwoTensor strain0dev2 = strain0dev * strain0dev;
     177           0 :   F_pos = 0.5 * k * strain0tr_pos * strain0tr_pos + mu * strain0dev2.trace();
     178           0 :   F_neg = 0.5 * k * strain0tr_neg * strain0tr_neg;
     179             : 
     180           0 :   _stress[_qp] = stress0pos * _D[_qp] - _pressure[_qp] * I2 * _I[_qp] + stress0neg;
     181             : 
     182             :   // 2nd derivative wrt c and strain = 0.0 if we used the previous step's history varible
     183           0 :   if (_use_current_hist)
     184           0 :     _d2Fdcdstrain[_qp] = stress0pos * _dDdc[_qp];
     185             : 
     186             :   // Used in StressDivergencePFFracTensors off-diagonal Jacobian
     187           0 :   _dstress_dc[_qp] = stress0pos * _dDdc[_qp] - _pressure[_qp] * I2 * _dIdc[_qp];
     188             : 
     189           0 :   if (strain0tr < 0)
     190           0 :     Jacobian_neg = k * I2I2;
     191           0 :   Jacobian_pos = _elasticity_tensor[_qp] - Jacobian_neg;
     192           0 :   _Jacobian_mult[_qp] = _D[_qp] * Jacobian_pos + Jacobian_neg;
     193           0 : }
     194             : 
     195             : void
     196           0 : ComputeLinearElasticPFFractureStress::computeQpStress()
     197             : {
     198             :   Real F_pos, F_neg;
     199           0 :   RankTwoTensor I2(RankTwoTensor::initIdentity);
     200             : 
     201           0 :   switch (_decomposition_type)
     202             :   {
     203           0 :     case Decomposition_type::strain_spectral:
     204           0 :       computeStrainSpectral(F_pos, F_neg);
     205             :       break;
     206           0 :     case Decomposition_type::strain_vol_dev:
     207           0 :       computeStrainVolDev(F_pos, F_neg);
     208             :       break;
     209           0 :     case Decomposition_type::stress_spectral:
     210           0 :       computeStressSpectral(F_pos, F_neg);
     211             :       break;
     212           0 :     default:
     213             :     {
     214           0 :       RankTwoTensor stress = _elasticity_tensor[_qp] * _mechanical_strain[_qp];
     215           0 :       F_pos = stress.doubleContraction(_mechanical_strain[_qp]) / 2.0;
     216           0 :       F_neg = 0.0;
     217           0 :       if (_use_current_hist)
     218           0 :         _d2Fdcdstrain[_qp] = stress * _dDdc[_qp];
     219             : 
     220           0 :       _stress[_qp] = _D[_qp] * stress - _pressure[_qp] * I2 * _I[_qp];
     221           0 :       _dstress_dc[_qp] = stress * _dDdc[_qp] - _pressure[_qp] * I2 * _dIdc[_qp];
     222           0 :       _Jacobian_mult[_qp] = _D[_qp] * _elasticity_tensor[_qp];
     223             :     }
     224             :   }
     225             : 
     226             :   // // Assign history variable
     227           0 :   Real hist_variable = _H_old[_qp];
     228           0 :   if (_use_snes_vi_solver)
     229             :   {
     230           0 :     _H[_qp] = F_pos;
     231             : 
     232           0 :     if (_use_current_hist)
     233             :       hist_variable = _H[_qp];
     234             :   }
     235             :   else
     236             :   {
     237           0 :     if (F_pos > _H_old[_qp])
     238           0 :       _H[_qp] = F_pos;
     239             :     else
     240           0 :       _H[_qp] = _H_old[_qp];
     241             : 
     242           0 :     if (_use_current_hist)
     243           0 :       hist_variable = _H[_qp];
     244             : 
     245           0 :     if (hist_variable < _barrier[_qp])
     246             :       hist_variable = _barrier[_qp];
     247             :   }
     248             : 
     249             :   // Elastic free energy density
     250           0 :   _E[_qp] =
     251           0 :       hist_variable * _D[_qp] + F_neg - _pressure[_qp] * _mechanical_strain[_qp].trace() * _I[_qp];
     252           0 :   _dEdc[_qp] =
     253           0 :       hist_variable * _dDdc[_qp] - _pressure[_qp] * _mechanical_strain[_qp].trace() * _dIdc[_qp];
     254           0 :   _d2Ed2c[_qp] = hist_variable * _d2Dd2c[_qp] -
     255           0 :                  _pressure[_qp] * _mechanical_strain[_qp].trace() * _d2Id2c[_qp];
     256           0 : }