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 "ComputeFiniteStrain.h"
      11             : #include "Assembly.h"
      12             : 
      13             : #include "libmesh/quadrature.h"
      14             : #include "libmesh/utility.h"
      15             : 
      16             : MooseEnum
      17       24776 : ComputeFiniteStrain::decompositionType()
      18             : {
      19       49552 :   return MooseEnum("TaylorExpansion EigenSolution HughesWinget", "TaylorExpansion");
      20             : }
      21             : 
      22             : registerMooseObject("SolidMechanicsApp", ComputeFiniteStrain);
      23             : 
      24             : InputParameters
      25       11896 : ComputeFiniteStrain::validParams()
      26             : {
      27       11896 :   InputParameters params = ComputeIncrementalStrainBase::validParams();
      28       11896 :   params.addClassDescription(
      29             :       "Compute a strain increment and rotation increment for finite strains.");
      30       23792 :   params.addParam<MooseEnum>("decomposition_method",
      31       23792 :                              ComputeFiniteStrain::decompositionType(),
      32             :                              "Methods to calculate the strain and rotation increments");
      33       11896 :   return params;
      34           0 : }
      35             : 
      36        8916 : ComputeFiniteStrain::ComputeFiniteStrain(const InputParameters & parameters)
      37             :   : ComputeIncrementalStrainBase(parameters),
      38        8916 :     _Fhat(_fe_problem.getMaxQps()),
      39       17832 :     _decomposition_method(getParam<MooseEnum>("decomposition_method").getEnum<DecompMethod>()),
      40        8916 :     _use_hw(_decomposition_method == DecompMethod::HughesWinget),
      41        8916 :     _def_grad_mid(_use_hw ? &declareProperty<RankTwoTensor>(_base_name + "def_grad_mid") : nullptr),
      42       17832 :     _f_bar(_use_hw ? &declareProperty<RankTwoTensor>(_base_name + "f_bar") : nullptr)
      43             : {
      44        8916 : }
      45             : 
      46             : void
      47     4922260 : ComputeFiniteStrain::computeProperties()
      48             : {
      49     4922260 :   RankTwoTensor ave_Fhat;
      50             :   Real ave_dfgrd_det = 0.0;
      51    40690014 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      52             :   {
      53             :     // Deformation gradient
      54             :     auto A = RankTwoTensor::initializeFromRows(
      55    35767754 :         (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
      56             : 
      57             :     // Old Deformation gradient
      58             :     auto Fbar = RankTwoTensor::initializeFromRows(
      59    35767754 :         (*_grad_disp_old[0])[_qp], (*_grad_disp_old[1])[_qp], (*_grad_disp_old[2])[_qp]);
      60             : 
      61             :     // Gauss point deformation gradient
      62    35767754 :     _deformation_gradient[_qp] = A;
      63    35767754 :     _deformation_gradient[_qp].addIa(1.0);
      64             : 
      65             :     // deformation gradient midpoint (for Hughes-Winget kinematics)
      66    35767754 :     if (_use_hw)
      67             :     {
      68       40288 :       (*_def_grad_mid)[_qp].setToIdentity();
      69       40288 :       (*_def_grad_mid)[_qp] += 0.5 * (A + Fbar);
      70             :     }
      71             : 
      72             :     // A = gradU - gradUold
      73    35767754 :     A -= Fbar;
      74             : 
      75             :     //_f_bar = dDu/Dx_o (for Hughes-Winget kinematics)
      76    35767754 :     if (_use_hw)
      77       40288 :       (*_f_bar)[_qp] = A;
      78             : 
      79             :     // Fbar = ( I + gradUold)
      80    35767754 :     Fbar.addIa(1.0);
      81             : 
      82             :     // Incremental deformation gradient _Fhat = I + A Fbar^-1
      83    35767754 :     _Fhat[_qp] = A * Fbar.inverse();
      84    35767754 :     _Fhat[_qp].addIa(1.0);
      85             : 
      86    35767754 :     if (_volumetric_locking_correction)
      87             :     {
      88             :       // Calculate average _Fhat for volumetric locking correction
      89    23506144 :       ave_Fhat += _Fhat[_qp] * _JxW[_qp] * _coord[_qp];
      90             : 
      91             :       // Average deformation gradient
      92    23506144 :       ave_dfgrd_det += _deformation_gradient[_qp].det() * _JxW[_qp] * _coord[_qp];
      93             :     }
      94             :   }
      95             : 
      96     4922260 :   if (_volumetric_locking_correction)
      97             :   {
      98             :     // needed for volumetric locking correction
      99     2959548 :     ave_Fhat /= _current_elem_volume;
     100             :     // average deformation gradient
     101     2959548 :     ave_dfgrd_det /= _current_elem_volume;
     102             :   }
     103             : 
     104    40689358 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     105             :   {
     106    35767180 :     if (_volumetric_locking_correction)
     107             :     {
     108             :       // Finalize volumetric locking correction
     109    23506144 :       _Fhat[_qp] *= std::cbrt(ave_Fhat.det() / _Fhat[_qp].det());
     110    23506144 :       _deformation_gradient[_qp] *= std::cbrt(ave_dfgrd_det / _deformation_gradient[_qp].det());
     111             :     }
     112             : 
     113    35767180 :     computeQpStrain();
     114             :   }
     115     4922178 : }
     116             : 
     117             : void
     118    44147832 : ComputeFiniteStrain::computeQpStrain()
     119             : {
     120    44147832 :   RankTwoTensor total_strain_increment;
     121             : 
     122             :   // three ways to calculate these increments: TaylorExpansion(default), EigenSolution, or
     123             :   // HughesWinget
     124    44147832 :   computeQpIncrements(total_strain_increment, _rotation_increment[_qp]);
     125             : 
     126    44147750 :   _strain_increment[_qp] = total_strain_increment;
     127             : 
     128             :   // Remove the eigenstrain increment
     129    44147750 :   subtractEigenstrainIncrementFromStrain(_strain_increment[_qp]);
     130             : 
     131    44147750 :   if (_dt > 0)
     132    44088326 :     _strain_rate[_qp] = _strain_increment[_qp] / _dt;
     133             :   else
     134       59424 :     _strain_rate[_qp].zero();
     135             : 
     136             :   // if HughesWinget, rotate old strains here
     137    44147750 :   RankTwoTensor mechanical_strain_old = _mechanical_strain_old[_qp];
     138    44147750 :   RankTwoTensor total_strain_old = _total_strain_old[_qp];
     139    44147750 :   if (_use_hw)
     140             :   {
     141       40288 :     mechanical_strain_old = _rotation_increment[_qp] * _mechanical_strain_old[_qp] *
     142       40288 :                             _rotation_increment[_qp].transpose();
     143       40288 :     total_strain_old =
     144       40288 :         _rotation_increment[_qp] * _total_strain_old[_qp] * _rotation_increment[_qp].transpose();
     145             :   }
     146             : 
     147             :   // Update strain in intermediate configuration
     148    44147750 :   _mechanical_strain[_qp] = mechanical_strain_old + _strain_increment[_qp];
     149    44147750 :   _total_strain[_qp] = total_strain_old + total_strain_increment;
     150             : 
     151             :   // Rotate strain to current configuration, unless HughesWinget
     152    44147750 :   if (!_use_hw)
     153             :   {
     154    44107462 :     _mechanical_strain[_qp] =
     155    44107462 :         _rotation_increment[_qp] * _mechanical_strain[_qp] * _rotation_increment[_qp].transpose();
     156    44107462 :     _total_strain[_qp] =
     157    44107462 :         _rotation_increment[_qp] * _total_strain[_qp] * _rotation_increment[_qp].transpose();
     158             :   }
     159             : 
     160    44147750 :   if (_global_strain)
     161           0 :     _total_strain[_qp] += (*_global_strain)[_qp];
     162    44147750 : }
     163             : 
     164             : void
     165    44145444 : ComputeFiniteStrain::computeQpIncrements(RankTwoTensor & total_strain_increment,
     166             :                                          RankTwoTensor & rotation_increment)
     167             : {
     168    44145444 :   switch (_decomposition_method)
     169             :   {
     170    43527060 :     case DecompMethod::TaylorExpansion:
     171             :     {
     172             :       // inverse of _Fhat
     173    43527060 :       const RankTwoTensor invFhat = _Fhat[_qp].inverse();
     174             : 
     175             :       // A = I - _Fhat^-1
     176    43527060 :       RankTwoTensor A(RankTwoTensor::initIdentity);
     177    43527060 :       A -= invFhat;
     178             : 
     179             :       // Cinv - I = A A^T - A - A^T;
     180    43527060 :       RankTwoTensor Cinv_I = A * A.transpose() - A - A.transpose();
     181             : 
     182             :       // strain rate D from Taylor expansion, Chat = (-1/2(Chat^-1 - I) + 1/4*(Chat^-1 - I)^2 + ...
     183    43527060 :       total_strain_increment = -Cinv_I * 0.5 + Cinv_I * Cinv_I * 0.25;
     184             : 
     185    43527060 :       const Real a[3] = {invFhat(1, 2) - invFhat(2, 1),
     186    43527060 :                          invFhat(2, 0) - invFhat(0, 2),
     187    43527060 :                          invFhat(0, 1) - invFhat(1, 0)};
     188             : 
     189    43527060 :       Real q = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) / 4.0;
     190    43527060 :       Real trFhatinv_1 = invFhat.trace() - 1.0;
     191    43527060 :       const Real p = trFhatinv_1 * trFhatinv_1 / 4.0;
     192             : 
     193             :       // cos theta_a
     194    43527060 :       const Real C1_squared = p +
     195    43527060 :                               3.0 * Utility::pow<2>(p) * (1.0 - (p + q)) / Utility::pow<2>(p + q) -
     196    43527060 :                               2.0 * Utility::pow<3>(p) * (1.0 - (p + q)) / Utility::pow<3>(p + q);
     197    43527060 :       if (C1_squared <= 0.0)
     198          82 :         mooseException(
     199             :             "Cannot take square root of a number less than or equal to zero in the calculation of "
     200             :             "C1 for the Rashid approximation for the rotation tensor. This zero or negative number "
     201             :             "may occur when elements become heavily distorted.");
     202             : 
     203    43526978 :       const Real C1 = std::sqrt(C1_squared);
     204             : 
     205             :       Real C2;
     206    43526978 :       if (q > 0.01)
     207             :         // (1-cos theta_a)/4q
     208      248680 :         C2 = (1.0 - C1) / (4.0 * q);
     209             :       else
     210             :         // alternate form for small q
     211    43278298 :         C2 = 0.125 + q * 0.03125 * (Utility::pow<2>(p) - 12.0 * (p - 1.0)) / Utility::pow<2>(p) +
     212    43278298 :              Utility::pow<2>(q) * (p - 2.0) * (Utility::pow<2>(p) - 10.0 * p + 32.0) /
     213             :                  Utility::pow<3>(p) +
     214    43278298 :              Utility::pow<3>(q) *
     215    43278298 :                  (1104.0 - 992.0 * p + 376.0 * Utility::pow<2>(p) - 72.0 * Utility::pow<3>(p) +
     216    43278298 :                   5.0 * Utility::pow<4>(p)) /
     217    43278298 :                  (512.0 * Utility::pow<4>(p));
     218             : 
     219             :       const Real C3_test =
     220    43526978 :           (p * q * (3.0 - q) + Utility::pow<3>(p) + Utility::pow<2>(q)) / Utility::pow<3>(p + q);
     221             : 
     222    43526978 :       if (C3_test <= 0.0)
     223           0 :         mooseException(
     224             :             "Cannot take square root of a number less than or equal to zero in the calculation of "
     225             :             "C3_test for the Rashid approximation for the rotation tensor. This zero or negative "
     226             :             "number may occur when elements become heavily distorted.");
     227             : 
     228    43526978 :       const Real C3 = 0.5 * std::sqrt(C3_test); // sin theta_a/(2 sqrt(q))
     229             : 
     230             :       // Calculate incremental rotation. Note that this value is the transpose of that from Rashid,
     231             :       // 93, so we transpose it before storing
     232    43526978 :       RankTwoTensor R_incr;
     233    43526978 :       R_incr.addIa(C1);
     234   174107912 :       for (unsigned int i = 0; i < 3; ++i)
     235   522323736 :         for (unsigned int j = 0; j < 3; ++j)
     236   391742802 :           R_incr(i, j) += C2 * a[i] * a[j];
     237             : 
     238    43526978 :       R_incr(0, 1) += C3 * a[2];
     239    43526978 :       R_incr(0, 2) -= C3 * a[1];
     240    43526978 :       R_incr(1, 0) -= C3 * a[2];
     241    43526978 :       R_incr(1, 2) += C3 * a[0];
     242    43526978 :       R_incr(2, 0) += C3 * a[1];
     243    43526978 :       R_incr(2, 1) -= C3 * a[0];
     244             : 
     245    43526978 :       rotation_increment = R_incr.transpose();
     246             :       break;
     247             :     }
     248             : 
     249      578096 :     case DecompMethod::EigenSolution:
     250             :     {
     251      578096 :       FactorizedRankTwoTensor Chat = RankTwoTensor::transposeTimes(_Fhat[_qp]);
     252      578096 :       FactorizedRankTwoTensor Uhat = MathUtils::sqrt(Chat);
     253     1156192 :       rotation_increment = _Fhat[_qp] * Uhat.inverse().get();
     254     1156192 :       total_strain_increment = MathUtils::log(Uhat).get();
     255             :       break;
     256             :     }
     257             : 
     258       40288 :     case DecompMethod::HughesWinget:
     259             :     {
     260       40288 :       const RankTwoTensor G = (*_f_bar)[_qp] * (*_def_grad_mid)[_qp].inverse();
     261             : 
     262       40288 :       total_strain_increment = 0.5 * (G + G.transpose());
     263       40288 :       const RankTwoTensor W = 0.5 * (G - G.transpose());
     264             : 
     265       40288 :       RankTwoTensor Q_1(RankTwoTensor::initIdentity);
     266       40288 :       RankTwoTensor Q_2(RankTwoTensor::initIdentity);
     267             : 
     268       40288 :       Q_1 -= 0.5 * W;
     269       40288 :       Q_2 += 0.5 * W;
     270             : 
     271       40288 :       rotation_increment = Q_1.inverse() * Q_2;
     272             : 
     273             :       break;
     274             :     }
     275             : 
     276           0 :     default:
     277           0 :       mooseError("ComputeFiniteStrain Error: Pass valid decomposition type: TaylorExpansion, "
     278             :                  "EigenSolution, or HughesWinget.");
     279             :   }
     280    44145362 : }