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             : GenericReal<is_ad>
     114    48490076 : InertialForceTempl<is_ad>::computeQpResidual()
     115             : {
     116    48490076 :   if (_dt == 0)
     117           0 :     return 0;
     118    48490076 :   else if (_has_beta)
     119             :   {
     120    45394848 :     auto accel = 1.0 / _beta *
     121    45336672 :                  (((_u[_qp] - (*_u_old)[_qp]) / (_dt * _dt)) - (*_vel_old)[_qp] / _dt -
     122    45336672 :                   (*_accel_old)[_qp] * (0.5 - _beta));
     123    45336672 :     auto vel =
     124    45453024 :         (*_vel_old)[_qp] + (_dt * (1.0 - _gamma)) * (*_accel_old)[_qp] + _gamma * _dt * accel;
     125    45336672 :     return _test[_i][_qp] * _density[_qp] *
     126    45453024 :            (accel + vel * _eta[_qp] * (1.0 + _alpha) - _alpha * _eta[_qp] * (*_vel_old)[_qp]);
     127             :   }
     128             : 
     129             :   // Lumped mass option
     130             :   // Only lumping the masses here
     131             :   // will multiply by corresponding residual multiplier after lumping the matrix
     132             :   // Density scaling is a fictitious added density to increase the time step
     133     3153404 :   else if (_time_integrator.isLumped() && _time_integrator.isExplicit() && !is_ad)
     134      107016 :     return _test[_i][_qp] * (_density[_qp] + _density_scaling[_qp]);
     135             : 
     136             :   // Consistent mass option
     137             :   // Same for explicit, implicit, and implicit with HHT
     138             :   else
     139     3046388 :     return _test[_i][_qp] * _density[_qp] *
     140     3046388 :            ((*_u_dotdot_factor)[_qp] + (*_u_dot_factor)[_qp] * _eta[_qp] * (1.0 + _alpha) -
     141     3046388 :             _alpha * _eta[_qp] * (*_u_dot_old)[_qp]);
     142             : }
     143             : 
     144             : template <>
     145             : void
     146           0 : InertialForceTempl<true>::computeResidualAdditional()
     147             : {
     148           0 :   mooseError("Internal error");
     149             : }
     150             : 
     151             : template <>
     152             : void
     153     1229729 : InertialForceTempl<false>::computeResidualAdditional()
     154             : {
     155     1229729 :   if (_dt == 0)
     156             :     return;
     157             : 
     158             :   // Explicit lumped only
     159     1229729 :   if (!_time_integrator.isLumped() || !_time_integrator.isExplicit())
     160     1226402 :     return;
     161             : 
     162       20007 :   for (unsigned int i = 0; i < _test.size(); ++i)
     163       16680 :     this->_local_re(i) *= (*_u_dotdot_factor_dof)[i] + _eta[0] * (*_u_dot_factor_dof)[i];
     164             : }
     165             : 
     166             : template <>
     167             : Real
     168           0 : InertialForceTempl<true>::computeQpJacobian()
     169             : {
     170           0 :   mooseError("Internal error");
     171             : }
     172             : 
     173             : template <>
     174             : Real
     175    76093368 : InertialForceTempl<false>::computeQpJacobian()
     176             : {
     177    76093368 :   if (_dt == 0)
     178             :     return 0;
     179             :   else
     180             :   {
     181    76093368 :     if (_has_beta)
     182    70763392 :       return _test[_i][_qp] * _density[_qp] / (_beta * _dt * _dt) * _phi[this->_j][_qp] +
     183    70763392 :              _eta[_qp] * (1 + _alpha) * _test[_i][_qp] * _density[_qp] * _gamma / _beta / _dt *
     184    70763392 :                  _phi[this->_j][_qp];
     185             :     else
     186     5329976 :       return _test[_i][_qp] * _density[_qp] * (*_du_dotdot_du)[_qp] * _phi[this->_j][_qp] +
     187     5329976 :              _eta[_qp] * (1 + _alpha) * _test[_i][_qp] * _density[_qp] * (*_du_dot_du)[_qp] *
     188     5329976 :                  _phi[this->_j][_qp];
     189             :   }
     190             : }
     191             : 
     192             : template class InertialForceTempl<false>;
     193             : template class InertialForceTempl<true>;