Line data    Source code

       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             : 
      10             : #include "ADComputeFiniteStrain.h"
      11             : #include "RankTwoTensor.h"
      12             : #include "RankFourTensor.h"
      13             : #include "SymmetricRankTwoTensor.h"
      14             : #include "SymmetricRankFourTensor.h"
      15             : 
      16             : #include "libmesh/quadrature.h"
      17             : #include "libmesh/utility.h"
      18             : 
      19             : registerMooseObject("TensorMechanicsApp", ADComputeFiniteStrain);
      20             : registerMooseObject("TensorMechanicsApp", ADSymmetricFiniteStrain);
      21             : 
      22             : template <typename R2, typename R4>
      23             : MooseEnum
      24        1440 : ADComputeFiniteStrainTempl<R2, R4>::decompositionType()
      25             : {
      26        2880 :   return MooseEnum("TaylorExpansion EigenSolution", "TaylorExpansion");
      27             : }
      28             : 
      29             : template <typename R2, typename R4>
      30             : InputParameters
      31        1440 : ADComputeFiniteStrainTempl<R2, R4>::validParams()
      32             : {
      33             :   InputParameters params = ADComputeIncrementalStrainBase::validParams();
      34        1440 :   params.addClassDescription(
      35             :       "Compute a strain increment and rotation increment for finite strains.");
      36        2880 :   params.addParam<MooseEnum>("decomposition_method",
      37             :                              ADComputeFiniteStrainTempl<R2, R4>::decompositionType(),
      38             :                              "Methods to calculate the strain and rotation increments");
      39        1440 :   return params;
      40           0 : }
      41             : 
      42             : template <typename R2, typename R4>
      43        1080 : ADComputeFiniteStrainTempl<R2, R4>::ADComputeFiniteStrainTempl(const InputParameters & parameters)
      44             :   : ADComputeIncrementalStrainBaseTempl<R2>(parameters),
      45        1080 :     _Fhat(this->_fe_problem.getMaxQps()),
      46        1080 :     _decomposition_method(
      47        2160 :         this->template getParam<MooseEnum>("decomposition_method").template getEnum<DecompMethod>())
      48             : {
      49        1080 : }
      50             : 
      51             : template <typename R2, typename R4>
      52             : void
      53      758500 : ADComputeFiniteStrainTempl<R2, R4>::computeProperties()
      54             : {
      55      758500 :   ADRankTwoTensor ave_Fhat;
      56     6355260 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      57             :   {
      58             :     // Deformation gradient
      59     5596760 :     auto A = ADRankTwoTensor::initializeFromRows(
      60     5596760 :         (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
      61             : 
      62             :     // Old Deformation gradient
      63     5596760 :     auto Fbar = ADRankTwoTensor::initializeFromRows(
      64     5596760 :         (*_grad_disp_old[0])[_qp], (*_grad_disp_old[1])[_qp], (*_grad_disp_old[2])[_qp]);
      65             : 
      66             :     // A = gradU - gradUold
      67     5596760 :     A -= Fbar;
      68             : 
      69             :     // Fbar = ( I + gradUold)
      70     5596760 :     Fbar.addIa(1.0);
      71             : 
      72             :     // Incremental deformation gradient _Fhat = I + A Fbar^-1
      73     5596760 :     _Fhat[_qp] = A * Fbar.inverse();
      74     5596760 :     _Fhat[_qp].addIa(1.0);
      75             : 
      76             :     // Calculate average _Fhat for volumetric locking correction
      77     5596760 :     if (_volumetric_locking_correction)
      78     3124176 :       ave_Fhat += _Fhat[_qp] * _JxW[_qp] * _coord[_qp];
      79             :   }
      80             : 
      81      758500 :   if (_volumetric_locking_correction)
      82      390522 :     ave_Fhat /= _current_elem_volume;
      83             : 
      84      758500 :   const auto ave_Fhat_det = ave_Fhat.det();
      85     6355257 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      86             :   {
      87             :     // Finalize volumetric locking correction
      88     5596758 :     if (_volumetric_locking_correction)
      89     6248352 :       _Fhat[_qp] *= std::cbrt(ave_Fhat_det / _Fhat[_qp].det());
      90             : 
      91     5596758 :     computeQpStrain();
      92             :   }
      93      758499 : }
      94             : 
      95             : template <typename R2, typename R4>
      96             : void
      97     6776044 : ADComputeFiniteStrainTempl<R2, R4>::computeQpStrain()
      98             : {
      99     6776044 :   ADR2 total_strain_increment;
     100             : 
     101             :   // two ways to calculate these increments: TaylorExpansion(default) or EigenSolution
     102     6776044 :   computeQpIncrements(total_strain_increment, _rotation_increment[_qp]);
     103             : 
     104     6776043 :   _strain_increment[_qp] = total_strain_increment;
     105             : 
     106             :   // Remove the eigenstrain increment
     107     6776043 :   this->subtractEigenstrainIncrementFromStrain(_strain_increment[_qp]);
     108             : 
     109     6776043 :   if (_dt > 0)
     110     6729967 :     _strain_rate[_qp] = _strain_increment[_qp] / _dt;
     111             :   else
     112       46076 :     _strain_rate[_qp].zero();
     113             : 
     114             :   // Update strain in intermediate configuration
     115     6776043 :   _mechanical_strain[_qp] = _mechanical_strain_old[_qp] + _strain_increment[_qp];
     116     6776043 :   _total_strain[_qp] = _total_strain_old[_qp] + total_strain_increment;
     117             : 
     118             :   // Rotate strain to current configuration
     119     6776043 :   _mechanical_strain[_qp].rotate(_rotation_increment[_qp]);
     120     6776043 :   _total_strain[_qp].rotate(_rotation_increment[_qp]);
     121             : 
     122     6776043 :   if (_global_strain)
     123           0 :     _total_strain[_qp] += (*_global_strain)[_qp];
     124     6776043 : }
     125             : 
     126             : template <typename R2, typename R4>
     127             : void
     128     6776044 : ADComputeFiniteStrainTempl<R2, R4>::computeQpIncrements(ADR2 & total_strain_increment,
     129             :                                                         ADRankTwoTensor & rotation_increment)
     130             : {
     131     6776044 :   switch (_decomposition_method)
     132             :   {
     133     6776044 :     case DecompMethod::TaylorExpansion:
     134             :     {
     135             :       // inverse of _Fhat
     136     6776044 :       const ADRankTwoTensor invFhat = _Fhat[_qp].inverse();
     137             : 
     138             :       // A = I - _Fhat^-1
     139     6776044 :       ADRankTwoTensor A(ADRankTwoTensor::initIdentity);
     140     6776044 :       A -= invFhat;
     141             : 
     142             :       // Cinv - I = A A^T - (A + A^T);
     143     6776044 :       ADR2 Cinv_I = ADR2::timesTranspose(A) - ADR2::plusTranspose(A);
     144             : 
     145             :       // strain rate D from Taylor expansion, Chat = (-1/2(Chat^-1 - I) + 1/4*(Chat^-1 - I)^2 + ...
     146     6780364 :       total_strain_increment = -Cinv_I * 0.5 + Cinv_I.square() * 0.25;
     147             : 
     148             :       const ADReal a[3] = {invFhat(1, 2) - invFhat(2, 1),
     149             :                            invFhat(2, 0) - invFhat(0, 2),
     150             :                            invFhat(0, 1) - invFhat(1, 0)};
     151             : 
     152     6776044 :       const auto q = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) / 4.0;
     153     6776044 :       const auto trFhatinv_1 = invFhat.trace() - 1.0;
     154    13552088 :       const auto p = trFhatinv_1 * trFhatinv_1 / 4.0;
     155             : 
     156             :       // cos theta_a
     157     6776044 :       const ADReal C1_squared =
     158    13552088 :           p + 3.0 * Utility::pow<2>(p) * (1.0 - (p + q)) / Utility::pow<2>(p + q) -
     159    20328132 :           2.0 * Utility::pow<3>(p) * (1.0 - (p + q)) / Utility::pow<3>(p + q);
     160     6776044 :       if (C1_squared <= 0.0)
     161           0 :         mooseException(
     162             :             "Cannot take square root of a number less than or equal to zero in the calculation of "
     163             :             "C1 for the Rashid approximation for the rotation tensor. This zero or negative number "
     164             :             "may occur when elements become heavily distorted.");
     165             : 
     166     6776044 :       const ADReal C1 = std::sqrt(C1_squared);
     167             : 
     168             :       ADReal C2;
     169     6776044 :       if (q > 0.01)
     170             :         // (1-cos theta_a)/4q
     171       34482 :         C2 = (1.0 - C1) / (4.0 * q);
     172             :       else
     173             :         // alternate form for small q
     174    27058200 :         C2 = 0.125 + q * 0.03125 * (Utility::pow<2>(p) - 12.0 * (p - 1.0)) / Utility::pow<2>(p) +
     175     6764550 :              Utility::pow<2>(q) * (p - 2.0) * (Utility::pow<2>(p) - 10.0 * p + 32.0) /
     176             :                  Utility::pow<3>(p) +
     177             :              Utility::pow<3>(q) *
     178    27058200 :                  (1104.0 - 992.0 * p + 376.0 * Utility::pow<2>(p) - 72.0 * Utility::pow<3>(p) +
     179     6764550 :                   5.0 * Utility::pow<4>(p)) /
     180    13529100 :                  (512.0 * Utility::pow<4>(p));
     181             : 
     182     6776045 :       const ADReal C3_test =
     183    13552088 :           (p * q * (3.0 - q) + Utility::pow<3>(p) + Utility::pow<2>(q)) / Utility::pow<3>(p + q);
     184     6776044 :       if (C3_test <= 0.0)
     185           1 :         mooseException(
     186             :             "Cannot take square root of a number less than or equal to zero in the calculation of "
     187             :             "C3_test for the Rashid approximation for the rotation tensor. This zero or negative "
     188             :             "number may occur when elements become heavily distorted.");
     189    13552086 :       const ADReal C3 = 0.5 * std::sqrt(C3_test); // sin theta_a/(2 sqrt(q))
     190             : 
     191             :       // Calculate incremental rotation. Note that this value is the transpose of that from Rashid,
     192             :       // 93, so we transpose it before storing
     193     6776043 :       ADRankTwoTensor R_incr;
     194     6776043 :       R_incr.addIa(C1);
     195    27104172 :       for (unsigned int i = 0; i < 3; ++i)
     196    81312516 :         for (unsigned int j = 0; j < 3; ++j)
     197   121968774 :           R_incr(i, j) += C2 * a[i] * a[j];
     198             : 
     199     6776043 :       R_incr(0, 1) += C3 * a[2];
     200     6776043 :       R_incr(0, 2) -= C3 * a[1];
     201     6776043 :       R_incr(1, 0) -= C3 * a[2];
     202     6776043 :       R_incr(1, 2) += C3 * a[0];
     203     6776043 :       R_incr(2, 0) += C3 * a[1];
     204     6776044 :       R_incr(2, 1) -= C3 * a[0];
     205             : 
     206     6776044 :       rotation_increment = R_incr.transpose();
     207             :       break;
     208             :     }
     209             : 
     210           0 :     case DecompMethod::EigenSolution:
     211             :     {
     212           0 :       FADR2 Chat = ADR2::transposeTimes(_Fhat[_qp]);
     213           0 :       FADR2 Uhat = MathUtils::sqrt(Chat);
     214           0 :       rotation_increment = _Fhat[_qp] * Uhat.inverse().template get<ADRankTwoTensor>();
     215           0 :       total_strain_increment = MathUtils::log(Uhat).template get<ADR2>();
     216             :       break;
     217             :     }
     218             : 
     219           0 :     default:
     220           0 :       mooseError("ADComputeFiniteStrain Error: Pass valid decomposition type: TaylorExpansion or "
     221             :                  "EigenSolution.");
     222             :   }
     223     6776043 : }
     224             : 
     225             : template class ADComputeFiniteStrainTempl<RankTwoTensor, RankFourTensor>;
     226             : template class ADComputeFiniteStrainTempl<SymmetricRankTwoTensor, SymmetricRankFourTensor>;