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 "InertialForce.h"
      11             : #include "SubProblem.h"
      12             : #include "TimeIntegrator.h"
      13             : #include "NonlinearSystemBase.h"
      14             : 
      15             : registerMooseObject("SolidMechanicsApp", InertialForce);
      16             : registerMooseObject("SolidMechanicsApp", ADInertialForce);
      17             : 
      18             : template <bool is_ad>
      19             : InputParameters
      20        1658 : InertialForceTempl<is_ad>::validParams()
      21             : {
      22        1658 :   InputParameters params = InertialForceParent<is_ad>::validParams();
      23        1658 :   params.addClassDescription("Calculates the residual for the inertial force "
      24             :                              "($M \\cdot acceleration$) and the contribution of mass"
      25             :                              " dependent Rayleigh damping and HHT time "
      26             :                              " integration scheme ($\\eta \\cdot M \\cdot"
      27             :                              " ((1+\\alpha)velq2-\\alpha \\cdot vel-old) $)");
      28        1658 :   params.set<bool>("use_displaced_mesh") = true;
      29        3316 :   params.addCoupledVar("velocity", "velocity variable");
      30        3316 :   params.addCoupledVar("acceleration", "acceleration variable");
      31        3316 :   params.addParam<Real>("beta", "beta parameter for Newmark Time integration");
      32        3316 :   params.addParam<Real>("gamma", "gamma parameter for Newmark Time integration");
      33        3316 :   params.addParam<MaterialPropertyName>("eta",
      34        3316 :                                         0.0,
      35             :                                         "Name of material property or a constant real "
      36             :                                         "number defining the eta parameter for the "
      37             :                                         "Rayleigh damping.");
      38        3316 :   params.addParam<MaterialPropertyName>(
      39             :       "density_scaling",
      40        3316 :       0.0,
      41             :       "Name of material property to add mass scaling in explicit simulations.");
      42        3316 :   params.addParam<Real>("alpha",
      43        3316 :                         0,
      44             :                         "alpha parameter for mass dependent numerical damping induced "
      45             :                         "by HHT time integration scheme");
      46        3316 :   params.addParam<MaterialPropertyName>(
      47             :       "density", "density", "Name of Material Property that provides the density");
      48        1658 :   return params;
      49           0 : }
      50             : 
      51             : template <bool is_ad>
      52         828 : InertialForceTempl<is_ad>::InertialForceTempl(const InputParameters & parameters)
      53             :   : InertialForceParent<is_ad>(parameters),
      54         828 :     _density(this->template getGenericMaterialProperty<Real, is_ad>("density")),
      55        1656 :     _has_beta(this->isParamValid("beta")),
      56        1656 :     _has_gamma(this->isParamValid("gamma")),
      57        1364 :     _beta(_has_beta ? this->template getParam<Real>("beta") : 0.1),
      58        1364 :     _gamma(_has_gamma ? this->template getParam<Real>("gamma") : 0.1),
      59        1656 :     _has_velocity(this->isParamValid("velocity")),
      60        1656 :     _has_acceleration(this->isParamValid("acceleration")),
      61        1656 :     _eta(this->template getGenericMaterialProperty<Real, is_ad>("eta")),
      62        1656 :     _density_scaling(this->template getMaterialProperty<Real>("density_scaling")),
      63        1656 :     _alpha(this->template getParam<Real>("alpha")),
      64        1656 :     _time_integrator(_sys.getTimeIntegrator(this->_var.number()))
      65             : {
      66         828 :   if (_has_beta && _has_gamma && _has_velocity && _has_acceleration)
      67             :   {
      68         536 :     _vel_old = &this->coupledValueOld("velocity");
      69         536 :     _accel_old = &this->coupledValueOld("acceleration");
      70         536 :     _u_old = &this->valueOld();
      71             :   }
      72         292 :   else if (!_has_beta && !_has_gamma && !_has_velocity && !_has_acceleration)
      73             :   {
      74         292 :     _u_dot_old = &(_var.uDotOld());
      75         292 :     _du_dot_du = &(_var.duDotDu());
      76         292 :     _du_dotdot_du = &(_var.duDotDotDu());
      77             : 
      78         292 :     this->addFEVariableCoupleableVectorTag(_time_integrator.uDotFactorTag());
      79         292 :     this->addFEVariableCoupleableVectorTag(_time_integrator.uDotDotFactorTag());
      80             : 
      81         292 :     _u_dot_factor = &_var.vectorTagValue(_time_integrator.uDotFactorTag());
      82         292 :     _u_dotdot_factor = &_var.vectorTagValue(_time_integrator.uDotDotFactorTag());
      83             : 
      84         292 :     if (_time_integrator.isLumped())
      85             :     {
      86         124 :       _u_dot_factor_dof = &_var.vectorTagDofValue(_time_integrator.uDotFactorTag());
      87         124 :       _u_dotdot_factor_dof = &_var.vectorTagDofValue(_time_integrator.uDotDotFactorTag());
      88             :     }
      89             :   }
      90             :   else
      91           0 :     mooseError("InertialForce: Either all or none of `beta`, `gamma`, `velocity`and `acceleration` "
      92             :                "should be provided as input.");
      93             : 
      94             :   // Check if HHT and explicit are being used simultaneously
      95         828 :   if (_alpha != 0 && _time_integrator.isExplicit())
      96           0 :     mooseError("InertialForce: HHT time integration parameter can only be used with Newmark-Beta "
      97             :                "time integration.");
      98             : 
      99             :   // Check if beta and explicit are being used simultaneously
     100         828 :   if (_has_beta && _time_integrator.isExplicit())
     101           0 :     mooseError("InertialForce: Newmark-beta integration parameter, beta, cannot be provided along "
     102             :                "with an explicit time "
     103             :                "integrator.");
     104             : 
     105         828 :   if (!(_time_integrator.isLumped() && _time_integrator.isExplicit()))
     106        1408 :     if (parameters.isParamSetByUser("density_scaling"))
     107           0 :       this->paramError(
     108             :           "density_scaling",
     109             :           "Density (mass) scaling can only be used in lumped mass, explicit simulations");
     110         828 : }
     111             : 
     112             : template <bool is_ad>
     113             : std::pair<GenericReal<is_ad>, GenericReal<is_ad>>
     114    45336672 : InertialForceTempl<is_ad>::computeNewmarkBetaVelAccel(const GenericReal<is_ad> & u,
     115             :                                                       const Real u_old,
     116             :                                                       const Real vel_old,
     117             :                                                       const Real accel_old,
     118             :                                                       const Real beta,
     119             :                                                       const Real gamma,
     120             :                                                       const Real dt)
     121             : {
     122       58176 :   const auto accel =
     123    45453024 :       1.0 / beta * (((u - u_old) / (dt * dt)) - vel_old / dt - accel_old * (0.5 - beta));
     124    45453024 :   const auto vel = vel_old + (dt * (1.0 - gamma)) * accel_old + gamma * dt * accel;
     125    45336672 :   return {vel, accel};
     126             : }
     127             : 
     128             : template <bool is_ad>
     129             : GenericReal<is_ad>
     130    48490076 : InertialForceTempl<is_ad>::computeQpResidual()
     131             : {
     132    48490076 :   if (_dt == 0)
     133           0 :     return 0;
     134    48490076 :   else if (_has_beta)
     135             :   {
     136    45336672 :     const auto [vel, accel] = computeNewmarkBetaVelAccel(
     137    45336672 :         _u[_qp], (*_u_old)[_qp], (*_vel_old)[_qp], (*_accel_old)[_qp], _beta, _gamma, _dt);
     138    45336672 :     return _test[_i][_qp] * _density[_qp] *
     139    45453024 :            (accel + vel * _eta[_qp] * (1.0 + _alpha) - _alpha * _eta[_qp] * (*_vel_old)[_qp]);
     140             :   }
     141             : 
     142             :   // Lumped mass option
     143             :   // Only lumping the masses here
     144             :   // will multiply by corresponding residual multiplier after lumping the matrix
     145             :   // Density scaling is a fictitious added density to increase the time step
     146     3153404 :   else if (_time_integrator.isLumped() && _time_integrator.isExplicit() && !is_ad)
     147      107016 :     return _test[_i][_qp] * (_density[_qp] + _density_scaling[_qp]);
     148             : 
     149             :   // Consistent mass option
     150             :   // Same for explicit, implicit, and implicit with HHT
     151             :   else
     152     3046388 :     return _test[_i][_qp] * _density[_qp] *
     153     3046388 :            ((*_u_dotdot_factor)[_qp] + (*_u_dot_factor)[_qp] * _eta[_qp] * (1.0 + _alpha) -
     154     3046388 :             _alpha * _eta[_qp] * (*_u_dot_old)[_qp]);
     155             : }
     156             : 
     157             : template <>
     158             : void
     159           0 : InertialForceTempl<true>::computeResidualAdditional()
     160             : {
     161           0 :   mooseError("Internal error");
     162             : }
     163             : 
     164             : template <>
     165             : void
     166     1229729 : InertialForceTempl<false>::computeResidualAdditional()
     167             : {
     168     1229729 :   if (_dt == 0)
     169             :     return;
     170             : 
     171             :   // Explicit lumped only
     172     1229729 :   if (!_time_integrator.isLumped() || !_time_integrator.isExplicit())
     173     1226402 :     return;
     174             : 
     175       20007 :   for (unsigned int i = 0; i < _test.size(); ++i)
     176       16680 :     this->_local_re(i) *= (*_u_dotdot_factor_dof)[i] + _eta[0] * (*_u_dot_factor_dof)[i];
     177             : }
     178             : 
     179             : template <>
     180             : Real
     181           0 : InertialForceTempl<true>::computeQpJacobian()
     182             : {
     183           0 :   mooseError("Internal error");
     184             : }
     185             : 
     186             : template <>
     187             : Real
     188    76093368 : InertialForceTempl<false>::computeQpJacobian()
     189             : {
     190    76093368 :   if (_dt == 0)
     191             :     return 0;
     192             :   else
     193             :   {
     194    76093368 :     if (_has_beta)
     195    70763392 :       return _test[_i][_qp] * _density[_qp] / (_beta * _dt * _dt) * _phi[this->_j][_qp] +
     196    70763392 :              _eta[_qp] * (1 + _alpha) * _test[_i][_qp] * _density[_qp] * _gamma / _beta / _dt *
     197    70763392 :                  _phi[this->_j][_qp];
     198             :     else
     199     5329976 :       return _test[_i][_qp] * _density[_qp] * (*_du_dotdot_du)[_qp] * _phi[this->_j][_qp] +
     200     5329976 :              _eta[_qp] * (1 + _alpha) * _test[_i][_qp] * _density[_qp] * (*_du_dot_du)[_qp] *
     201     5329976 :                  _phi[this->_j][_qp];
     202             :   }
     203             : }
     204             : 
     205             : template class InertialForceTempl<false>;
     206             : template class InertialForceTempl<true>;