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 "ComputeFiniteStrainNOSPD.h"
      11             : 
      12             : #include "libmesh/utility.h"
      13             : 
      14             : registerMooseObject("PeridynamicsApp", ComputeFiniteStrainNOSPD);
      15             : 
      16             : MooseEnum
      17          23 : ComputeFiniteStrainNOSPD::decompositionType()
      18             : {
      19          46 :   return MooseEnum("TaylorExpansion EigenSolution", "TaylorExpansion");
      20             : }
      21             : 
      22             : InputParameters
      23          23 : ComputeFiniteStrainNOSPD::validParams()
      24             : {
      25          23 :   InputParameters params = ComputeStrainBaseNOSPD::validParams();
      26          23 :   params.addClassDescription(
      27             :       "Class for computing nodal quantities for residual and jacobian calculation "
      28             :       "for peridynamic correspondence models under finite strain assumptions");
      29             : 
      30          46 :   params.addParam<MooseEnum>("decomposition_method",
      31          46 :                              ComputeFiniteStrainNOSPD::decompositionType(),
      32             :                              "Methods to calculate the strain and rotation increments");
      33             : 
      34          23 :   return params;
      35           0 : }
      36             : 
      37          18 : ComputeFiniteStrainNOSPD::ComputeFiniteStrainNOSPD(const InputParameters & parameters)
      38             :   : ComputeStrainBaseNOSPD(parameters),
      39          18 :     _strain_rate(declareProperty<RankTwoTensor>("strain_rate")),
      40          18 :     _strain_increment(declareProperty<RankTwoTensor>("strain_increment")),
      41          18 :     _rotation_increment(declareProperty<RankTwoTensor>("rotation_increment")),
      42          36 :     _deformation_gradient_old(getMaterialPropertyOld<RankTwoTensor>("deformation_gradient")),
      43          36 :     _mechanical_strain_old(getMaterialPropertyOld<RankTwoTensor>("mechanical_strain")),
      44          36 :     _total_strain_old(getMaterialPropertyOld<RankTwoTensor>("total_strain")),
      45          18 :     _eigenstrains_old(_eigenstrain_names.size()),
      46          18 :     _Fhat(2),
      47          54 :     _decomposition_method(getParam<MooseEnum>("decomposition_method").getEnum<DecompMethod>())
      48             : {
      49          18 :   for (unsigned int i = 0; i < _eigenstrains_old.size(); ++i)
      50           0 :     _eigenstrains_old[i] = &getMaterialPropertyOld<RankTwoTensor>(_eigenstrain_names[i]);
      51          18 : }
      52             : 
      53             : void
      54      155448 : ComputeFiniteStrainNOSPD::computeQpStrain()
      55             : {
      56      155448 :   computeQpDeformationGradient();
      57             : 
      58      155448 :   computeQpFhat();
      59             : 
      60      155448 :   RankTwoTensor total_strain_increment;
      61             : 
      62             :   // Two ways to calculate these increments: TaylorExpansion(default) or EigenSolution
      63      155448 :   computeQpStrainRotationIncrements(total_strain_increment, _rotation_increment[_qp]);
      64             : 
      65      155448 :   _strain_increment[_qp] = total_strain_increment;
      66             : 
      67             :   // Remove the eigenstrain increment
      68      155448 :   subtractEigenstrainIncrementFromStrain(_strain_increment[_qp]);
      69             : 
      70      155448 :   if (_dt > 0)
      71      155448 :     _strain_rate[_qp] = _strain_increment[_qp] / _dt;
      72             :   else
      73           0 :     _strain_rate[_qp].zero();
      74             : 
      75             :   // Update strain in intermediate configuration
      76      155448 :   _mechanical_strain[_qp] = _mechanical_strain_old[_qp] + _strain_increment[_qp];
      77      155448 :   _total_strain[_qp] = _total_strain_old[_qp] + total_strain_increment;
      78             : 
      79             :   // Rotate strain to current configuration
      80      155448 :   _mechanical_strain[_qp] =
      81      155448 :       _rotation_increment[_qp] * _mechanical_strain[_qp] * _rotation_increment[_qp].transpose();
      82      155448 :   _total_strain[_qp] =
      83      155448 :       _rotation_increment[_qp] * _total_strain[_qp] * _rotation_increment[_qp].transpose();
      84             : 
      85             :   // zero out all strain measures for broken bond
      86      155448 :   if (_bond_status_var->getElementalValue(_current_elem) < 0.5)
      87             :   {
      88           0 :     _strain_rate[_qp].zero();
      89           0 :     _strain_increment[_qp].zero();
      90           0 :     _rotation_increment[_qp].zero();
      91           0 :     _mechanical_strain[_qp].zero();
      92           0 :     _total_strain[_qp].zero();
      93             :   }
      94      155448 : }
      95             : 
      96             : void
      97           0 : ComputeFiniteStrainNOSPD::computeQpFhat()
      98             : {
      99             :   // Incremental deformation gradient of current step w.r.t previous step:
     100             :   // _Fhat = deformation_gradient * inv(deformation_gradient_old)
     101           0 :   _Fhat[_qp] = _deformation_gradient[_qp] * _deformation_gradient_old[_qp].inverse();
     102           0 : }
     103             : 
     104             : void
     105      155448 : ComputeFiniteStrainNOSPD::computeQpStrainRotationIncrements(RankTwoTensor & total_strain_increment,
     106             :                                                             RankTwoTensor & rotation_increment)
     107             : {
     108      155448 :   switch (_decomposition_method)
     109             :   {
     110      155448 :     case DecompMethod::TaylorExpansion:
     111             :     {
     112             :       // inverse of _Fhat
     113      155448 :       RankTwoTensor invFhat(_Fhat[_qp].inverse());
     114             : 
     115             :       // A = I - _Fhat^-1
     116      155448 :       RankTwoTensor A(RankTwoTensor::initIdentity);
     117      155448 :       A -= invFhat;
     118             : 
     119             :       // Cinv - I = A A^T - A - A^T;
     120      155448 :       RankTwoTensor Cinv_I = A * A.transpose() - A - A.transpose();
     121             : 
     122             :       // strain rate D from Taylor expansion, Chat = (-1/2(Chat^-1 - I) + 1/4*(Chat^-1 - I)^2 + ...
     123      155448 :       total_strain_increment = -Cinv_I * 0.5 + Cinv_I * Cinv_I * 0.25;
     124             : 
     125      155448 :       const Real a[3] = {invFhat(1, 2) - invFhat(2, 1),
     126      155448 :                          invFhat(2, 0) - invFhat(0, 2),
     127      155448 :                          invFhat(0, 1) - invFhat(1, 0)};
     128             : 
     129      155448 :       Real q = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) / 4.0;
     130      155448 :       Real trFhatinv_1 = invFhat.trace() - 1.0;
     131      155448 :       const Real p = trFhatinv_1 * trFhatinv_1 / 4.0;
     132             : 
     133             :       // cos theta_a
     134      155448 :       const Real C1_squared = p +
     135      155448 :                               3.0 * Utility::pow<2>(p) * (1.0 - (p + q)) / Utility::pow<2>(p + q) -
     136      155448 :                               2.0 * Utility::pow<3>(p) * (1.0 - (p + q)) / Utility::pow<3>(p + q);
     137      155448 :       if (C1_squared <= 0.0)
     138           0 :         mooseException(
     139             :             "Cannot take square root of a number less than or equal to zero in the calculation of "
     140             :             "C1 for the Rashid approximation for the rotation tensor.");
     141             : 
     142      155448 :       const Real C1 = std::sqrt(C1_squared);
     143             : 
     144             :       Real C2;
     145      155448 :       if (q > 0.01)
     146             :         // (1-cos theta_a)/4q
     147           0 :         C2 = (1.0 - C1) / (4.0 * q);
     148             :       else
     149             :         // alternate form for small q
     150      155448 :         C2 = 0.125 + q * 0.03125 * (Utility::pow<2>(p) - 12.0 * (p - 1.0)) / Utility::pow<2>(p) +
     151      155448 :              Utility::pow<2>(q) * (p - 2.0) * (Utility::pow<2>(p) - 10.0 * p + 32.0) /
     152             :                  Utility::pow<3>(p) +
     153      155448 :              Utility::pow<3>(q) *
     154      155448 :                  (1104.0 - 992.0 * p + 376.0 * Utility::pow<2>(p) - 72.0 * Utility::pow<3>(p) +
     155      155448 :                   5.0 * Utility::pow<4>(p)) /
     156      155448 :                  (512.0 * Utility::pow<4>(p));
     157             : 
     158             :       const Real C3_test =
     159      155448 :           (p * q * (3.0 - q) + Utility::pow<3>(p) + Utility::pow<2>(q)) / Utility::pow<3>(p + q);
     160             : 
     161      155448 :       if (C3_test <= 0.0)
     162           0 :         mooseException(
     163             :             "Cannot take square root of a number less than or equal to zero in the calculation of "
     164             :             "C3_test for the Rashid approximation for the rotation tensor.");
     165             : 
     166      155448 :       const Real C3 = 0.5 * std::sqrt(C3_test); // sin theta_a/(2 sqrt(q))
     167             : 
     168             :       // Calculate incremental rotation. Note that this value is the transpose of that from Rashid,
     169             :       // 93, so we transpose it before storing
     170      155448 :       RankTwoTensor R_incr;
     171      155448 :       R_incr.addIa(C1);
     172      621792 :       for (unsigned int i = 0; i < 3; ++i)
     173     1865376 :         for (unsigned int j = 0; j < 3; ++j)
     174     1399032 :           R_incr(i, j) += C2 * a[i] * a[j];
     175             : 
     176      155448 :       R_incr(0, 1) += C3 * a[2];
     177      155448 :       R_incr(0, 2) -= C3 * a[1];
     178      155448 :       R_incr(1, 0) -= C3 * a[2];
     179      155448 :       R_incr(1, 2) += C3 * a[0];
     180      155448 :       R_incr(2, 0) += C3 * a[1];
     181      155448 :       R_incr(2, 1) -= C3 * a[0];
     182             : 
     183      155448 :       rotation_increment = R_incr.transpose();
     184             :       break;
     185             :     }
     186           0 :     case DecompMethod::EigenSolution:
     187             :     {
     188           0 :       FactorizedRankTwoTensor Chat = RankTwoTensor::transposeTimes(_Fhat[_qp]);
     189           0 :       FactorizedRankTwoTensor Uhat = MathUtils::sqrt(Chat);
     190           0 :       rotation_increment = _Fhat[_qp] * Uhat.inverse().get();
     191           0 :       total_strain_increment = MathUtils::log(Uhat).get();
     192             :       break;
     193             :     }
     194             : 
     195           0 :     default:
     196           0 :       mooseError("ComputeFiniteStrainNOSPD Error: Invalid decomposition type! Please specify : "
     197             :                  "TaylorExpansion or "
     198             :                  "EigenSolution.");
     199             :   }
     200      155448 : }
     201             : 
     202             : void
     203      155448 : ComputeFiniteStrainNOSPD::subtractEigenstrainIncrementFromStrain(RankTwoTensor & strain)
     204             : {
     205      155448 :   for (unsigned int i = 0; i < _eigenstrains.size(); ++i)
     206             :   {
     207           0 :     strain -= (*_eigenstrains[i])[_qp];
     208           0 :     strain += (*_eigenstrains_old[i])[_qp];
     209             :   }
     210      155448 : }