Line data    Source code

       1             : /*************************************************/
       2             : /*           DO NOT MODIFY THIS HEADER           */
       3             : /*                                               */
       4             : /*                     MASTODON                  */
       5             : /*                                               */
       6             : /*    (c) 2015 Battelle Energy Alliance, LLC     */
       7             : /*            ALL RIGHTS RESERVED                */
       8             : /*                                               */
       9             : /*   Prepared by Battelle Energy Alliance, LLC   */
      10             : /*     With the U. S. Department of Energy       */
      11             : /*                                               */
      12             : /*     See COPYRIGHT for full restrictions       */
      13             : /*************************************************/
      14             : #include "Function.h"
      15             : #include "MooseError.h"
      16             : #include "MooseMesh.h"
      17             : #include "NonReflectingBC.h"
      18             : 
      19             : registerMooseObject("MastodonApp", NonReflectingBC);
      20             : 
      21             : InputParameters
      22          12 : NonReflectingBC::validParams()
      23             : {
      24          12 :   InputParameters params = IntegratedBC::validParams();
      25          12 :   params.addClassDescription("Applies Lysmer damper in the normal and "
      26             :                              "tangential directions to soil boundary.");
      27          12 :   params += NonReflectingBC::commonParameters();
      28          24 :   params.addRequiredParam<unsigned int>("component",
      29             :                                         "The direction in which the Lysmer damper is applied.");
      30          12 :   params.set<bool>("use_displaced_mesh") = true;
      31          12 :   return params;
      32           0 : }
      33             : 
      34             : InputParameters
      35          14 : NonReflectingBC::commonParameters()
      36             : {
      37          14 :   InputParameters params = emptyInputParameters();
      38          28 :   params.addCoupledVar("displacements",
      39             :                        "The vector of displacement variables. "
      40             :                        "The size of this vector must be same "
      41             :                        "as the number of dimensions.");
      42          28 :   params.addCoupledVar("velocities",
      43             :                        "The vector of velocity variables that "
      44             :                        "are coupled to the displacement "
      45             :                        "variables. The size of this vector must "
      46             :                        "be same as that of displacements.");
      47          28 :   params.addCoupledVar("accelerations",
      48             :                        "The vector of acceleration variables that are coupled "
      49             :                        "to the displacement variables. The size of this vector "
      50             :                        "must be same as that of displacements.");
      51          28 :   params.addRequiredParam<Real>("beta", "The beta parameter for Newmark time integration.");
      52          28 :   params.addRequiredParam<Real>("gamma", "The gamma parameter for Newmark time integration.");
      53          28 :   params.addParam<Real>("alpha", 0.0, "The alpha parameter for HHT time integration.");
      54          28 :   params.addRequiredRangeCheckedParam<Real>("density", "density>0.0", "Density of the material.");
      55          28 :   params.addRequiredRangeCheckedParam<Real>(
      56             :       "p_wave_speed", "p_wave_speed>0.0", "P-wave speed of the material.");
      57          28 :   params.addRequiredRangeCheckedParam<Real>(
      58             :       "shear_wave_speed", "shear_wave_speed>0.0", "shear wave speed of the material.");
      59          14 :   return params;
      60           0 : }
      61             : 
      62          10 : NonReflectingBC::NonReflectingBC(const InputParameters & parameters)
      63             :   : IntegratedBC(parameters),
      64          10 :     _component(getParam<unsigned int>("component")),
      65          10 :     _ndisp(coupledComponents("displacements")),
      66          10 :     _disp(3),
      67          10 :     _disp_var(3),
      68          10 :     _disp_old(3),
      69          10 :     _vel_old(3),
      70          10 :     _accel_old(3),
      71          20 :     _beta(getParam<Real>("beta")),
      72          20 :     _gamma(getParam<Real>("gamma")),
      73          20 :     _alpha(getParam<Real>("alpha")),
      74          20 :     _density(getParam<Real>("density")),
      75          20 :     _p_wave_speed(getParam<Real>("p_wave_speed")),
      76          30 :     _shear_wave_speed(getParam<Real>("shear_wave_speed"))
      77             : {
      78             : 
      79             :   // Error checking on variable vectors
      80          10 :   if (_ndisp != _mesh.dimension())
      81           1 :     mooseError("The number of variables listed in the 'displacements' parameter in \"",
      82             :                name(),
      83             :                "\" block must match the mesh dimension.");
      84             : 
      85          18 :   if (coupledComponents("velocities") != _mesh.dimension())
      86           1 :     mooseError("The number of variables listed in the 'velocities' parameter in \"",
      87             :                name(),
      88             :                "\" block must match the mesh dimension.");
      89             : 
      90          16 :   if (coupledComponents("accelerations") != _mesh.dimension())
      91           1 :     mooseError("The number of variables listed in the 'accelerations' parameter in \"",
      92             :                name(),
      93             :                "\" block must match the mesh dimension.");
      94             : 
      95           7 :   if (_component >= _mesh.dimension())
      96           1 :     mooseError(
      97             :         "The 'component' parameter in \"", name(), "\" block should be less than mesh dimension.");
      98             : 
      99             :   // Populate coupled variable information
     100          24 :   for (unsigned int i = 0; i < _ndisp; ++i)
     101             :   {
     102          18 :     _disp[i] = &coupledValue("displacements", i);
     103          18 :     _disp_var[i] = coupled("displacements", i);
     104          18 :     _disp_old[i] = &coupledValueOld("displacements", i);
     105          18 :     _vel_old[i] = &coupledValueOld("velocities", i);
     106          18 :     _accel_old[i] = &coupledValueOld("accelerations", i);
     107             :   }
     108           6 : }
     109             : 
     110             : Real
     111     1038336 : NonReflectingBC::computeQpResidual()
     112             : {
     113     1038336 :   std::vector<Real> vel(3, 0.0);
     114             :   Real accel(0.0);
     115             : 
     116             :   Real normal_vel = 0.0;
     117     4153344 :   for (unsigned int i = 0; i < _ndisp; ++i)
     118             :   {
     119     3115008 :     accel = 1. / _beta * (((*_disp[i])[_qp] - (*_disp_old[i])[_qp]) / (_dt * _dt) -
     120     3115008 :                           (*_vel_old[i])[_qp] / _dt - (*_accel_old[i])[_qp] * (0.5 - _beta));
     121     3115008 :     vel[i] =
     122     3115008 :         (*_vel_old[i])[_qp] + (_dt * (1. - _gamma)) * (*_accel_old[i])[_qp] + _gamma * _dt * accel;
     123     3115008 :     vel[i] = (1. + _alpha) * vel[i] - _alpha * (*_vel_old[i])[_qp]; // HHT time integration
     124     3115008 :     normal_vel += vel[i] * _normals[_qp](i);
     125             :   }
     126             :   // residual is test[i][_qp] *( density* V_p * normal component of velocity +
     127             :   // density * V_s* tangential component of velocity)
     128     1038336 :   return _test[_i][_qp] * _density *
     129     1038336 :          (_p_wave_speed * normal_vel * _normals[_qp](_component) +
     130     1038336 :           _shear_wave_speed * (vel[_component] - normal_vel * _normals[_qp](_component)));
     131     1038336 : }
     132             : 
     133             : Real
     134     2555904 : NonReflectingBC::computeQpJacobian()
     135             : {
     136     2555904 :   return _test[_i][_qp] * _density *
     137     2555904 :          ((_p_wave_speed - _shear_wave_speed) * _normals[_qp](_component) *
     138     2555904 :               _normals[_qp](_component) +
     139     2555904 :           _shear_wave_speed) *
     140     2555904 :          (1. + _alpha) * _gamma / _beta / _dt * _phi[_j][_qp];
     141             : }
     142             : 
     143             : Real
     144           0 : NonReflectingBC::computeQpOffDiagJacobian(unsigned int jvar)
     145             : {
     146             :   unsigned int coupled_component = 0;
     147             :   bool active(false);
     148             : 
     149           0 :   for (unsigned int i = 0; i < _ndisp; ++i)
     150           0 :     if (jvar == _disp_var[i])
     151             :     {
     152             :       coupled_component = i;
     153             :       active = true;
     154             :     }
     155             : 
     156           0 :   if (active)
     157           0 :     return _test[_i][_qp] * _density * (_p_wave_speed - _shear_wave_speed) *
     158           0 :            _normals[_qp](_component) * _normals[_qp](coupled_component) * (1. + _alpha) * _gamma /
     159           0 :            _beta / _dt * _phi[_j][_qp];
     160             : 
     161             :   return 0.0;
     162             : }