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 "MooseMesh.h"
      16             : #include "SeismicSource.h"
      17             : 
      18             : /**
      19             :  *  This class applies a force at the given point. The direction of the force is
      20             :  *  determined by the strike and dip angle of the fault and the slip direction
      21             :  *(rake).
      22             :  *  The magnitude of the force is dependent on the seismic moment
      23             :  *(shear_modulus*area*slip).
      24             :  *  Strike = 0 implies x axis is aligned with geographic North. If a 3-D model
      25             :  *is not being used,
      26             :  *  please adjust the strike, rake and dip angles accordingly.
      27             :  *  Reference: Quantitative Seismology by Keiiti Aki and Paul G. Richards,
      28             :  *Chapter 4.
      29             :  **/
      30             : 
      31             : registerMooseObject("MastodonApp", SeismicSource);
      32             : 
      33             : InputParameters
      34          40 : SeismicSource::validParams()
      35             : {
      36          40 :   InputParameters params = DiracKernel::validParams();
      37          40 :   params.addClassDescription("This class applies a seismic source at a given point.");
      38          40 :   params += SeismicSource::commonParameters();
      39             : 
      40          80 :   params.addRequiredParam<unsigned int>("component",
      41             :                                         "The direction in which the force is applied.");
      42          40 :   return params;
      43           0 : }
      44             : 
      45             : InputParameters
      46          60 : SeismicSource::commonParameters()
      47             : {
      48          60 :   InputParameters params = emptyInputParameters();
      49         120 :   params.addRequiredParam<Real>("strike", "The strike angle of the fault in degrees.");
      50         120 :   params.addRequiredParam<Real>("dip", "The dip angle of the fault in degrees.");
      51         120 :   params.addRequiredParam<Real>("rake", "The rake angle of the fault in degrees.");
      52         120 :   params.addRequiredParam<Real>("area", "The area over which slip is distributed.");
      53         120 :   params.addRequiredParam<FunctionName>("slip",
      54             :                                         "The function describing slip as a function of time.");
      55         120 :   params.addRequiredParam<Real>("shear_modulus",
      56             :                                 "The shear modulus of the soil around the seismic source.");
      57         120 :   params.addParam<Real>("alpha", 0, "The Hilber Hughes Taylor (HHT) time integration parameter.");
      58         120 :   params.addParam<std::vector<Real>>("point",
      59             :                                      "The x, y and z coordinate of a single point source.");
      60         120 :   params.addParam<std::vector<FunctionName>>(
      61             :       "position_function",
      62             :       "The vector of function names that describes the x, "
      63             :       "y and z coordinates of the source. Only the x position"
      64             :       "function is required for 1 dimensional problems, x and y"
      65             :       "positions are required for 2 dimensional problems and x, "
      66             :       "y and z positions are required for 3D problems. The first"
      67             :       "column in these functions gives source number (starting"
      68             :       "with 1) and the second column contains the coordinate.");
      69         120 :   params.addParam<unsigned int>("number",
      70             :                                 "The number of point sources. This number should be same as "
      71             :                                 "number of entries in the function files describing the "
      72             :                                 "position of the point source.");
      73         120 :   params.addParam<Real>("rupture_speed",
      74             :                         "The speed at which the earthquake rupture propogates "
      75             :                         "through the fault (usually around 0.8 * shear wave "
      76             :                         "speed).");
      77         120 :   params.addParam<std::vector<Real>>("epicenter", "The x, y and z coordinates of the epicenter.");
      78          60 :   return params;
      79           0 : }
      80             : 
      81          31 : SeismicSource::SeismicSource(const InputParameters & parameters)
      82             :   : DiracKernel(parameters),
      83          31 :     _component(getParam<unsigned int>("component")),
      84          62 :     _alpha(getParam<Real>("alpha")),
      85          62 :     _shear_modulus(getParam<Real>("shear_modulus")),
      86          62 :     _area(getParam<Real>("area")),
      87          31 :     _slip_function(&getFunction("slip")),
      88          31 :     _moment_magnitude(0.0),
      89          31 :     _force(0.0),
      90          31 :     _rupture_speed(std::numeric_limits<unsigned int>::max()),
      91          62 :     _epicenter(_mesh.dimension(), 0.0)
      92             : {
      93         118 :   if (!isParamValid("point") && !isParamValid("position_function"))
      94           2 :     mooseError("Error in " + name() +
      95             :                ". Either a point of a set of points should be given as input.");
      96             : 
      97         113 :   if (isParamValid("position_function") && !isParamValid("number"))
      98           2 :     mooseError("Error in " + name() +
      99             :                ". Please specify number of source points defined in the function.");
     100             : 
     101          58 :   if (isParamValid("position_function"))
     102             :   {
     103             :     std::vector<FunctionName> pos_function =
     104          51 :         getParam<std::vector<FunctionName>>("position_function");
     105          17 :     if (pos_function.size() != _mesh.dimension())
     106           2 :       mooseError("Error in " + name() + ". The number of position functions "
     107             :                                         "should be equal to mesh dimension.");
     108          16 :   }
     109             : 
     110          28 :   if (_component >= _mesh.dimension())
     111           0 :     mooseError("Error in " + name() + ". component cannot exceed mesh dimensions.");
     112             : 
     113             :   // calulating seismic moment tensor
     114             :   Real pi = libMesh::pi;
     115          56 :   Real dip = getParam<Real>("dip") / 180.0 * pi;
     116          56 :   Real strike = getParam<Real>("strike") / 180.0 * pi;
     117          56 :   Real rake = getParam<Real>("rake") / 180.0 * pi;
     118             : 
     119          28 :   if (_mesh.dimension() == 3)
     120             :   {
     121          14 :     _moment.resize(3, std::vector<Real>(3, 0.0));
     122          14 :     _moment[0][0] = -(sin(dip) * cos(rake) * sin(2.0 * strike) +
     123          14 :                       sin(2.0 * dip) * sin(rake) * sin(strike) * sin(strike));
     124          14 :     _moment[0][1] = (sin(dip) * cos(rake) * cos(2.0 * strike) +
     125          14 :                      0.5 * sin(2.0 * dip) * sin(rake) * sin(2.0 * strike));
     126          14 :     _moment[1][0] = _moment[0][1];
     127          14 :     _moment[0][2] =
     128          14 :         -(cos(dip) * cos(rake) * cos(strike) + cos(2.0 * dip) * sin(rake) * sin(strike));
     129          14 :     _moment[2][0] = _moment[0][2];
     130          14 :     _moment[1][1] = (sin(dip) * cos(rake) * sin(2.0 * strike) -
     131          14 :                      sin(2.0 * dip) * sin(rake) * cos(strike) * cos(strike));
     132          14 :     _moment[1][2] =
     133          14 :         -(cos(dip) * cos(rake) * sin(strike) - cos(2.0 * dip) * sin(rake) * cos(strike));
     134          14 :     _moment[2][1] = _moment[1][2];
     135          14 :     _moment[2][2] = sin(2.0 * dip) * sin(rake);
     136             :   }
     137          14 :   else if (_mesh.dimension() == 2)
     138             :   {
     139             :     // x direction is aligned with north and y direction is aligned with vertical.
     140             :     // This will result in an in-plane earthquake wave.
     141          13 :     if (strike != 0.0)
     142           2 :       mooseError("Error in " + name() +
     143             :                  ". A non-zero strike angle for 2D models will create an "
     144             :                  "out-of-plane earthquake wave. This is currently not "
     145             :                  "supported.");
     146          12 :     _moment.resize(2, std::vector<Real>(2, 0.0));
     147          12 :     _moment[0][0] = (sin(dip) * cos(rake) * sin(2.0 * strike) -
     148          12 :                      sin(2.0 * dip) * sin(rake) * cos(strike) * cos(strike));
     149          12 :     _moment[0][1] =
     150          12 :         -(cos(dip) * cos(rake) * sin(strike) - cos(2.0 * dip) * sin(rake) * cos(strike));
     151          12 :     _moment[1][0] = _moment[0][1];
     152          12 :     _moment[1][1] = sin(2.0 * dip) * sin(rake);
     153             :   }
     154             :   else
     155           2 :     mooseError("Error in " + name() + ". Only mesh dimensions of 2 and 3 are currently supported.");
     156             : 
     157          72 :   if (isParamValid("rupture_speed") && isParamValid("epicenter"))
     158             :   {
     159          18 :     if (isParamValid("point"))
     160           2 :       mooseError("Error in " + name() +
     161             :                  ". Rupture speed and epicenter should only be provided "
     162             :                  "when multiple point sources are specified.");
     163             : 
     164          24 :     std::vector<Real> epi = getParam<std::vector<Real>>("epicenter");
     165           8 :     if (epi.size() != _mesh.dimension())
     166           2 :       mooseError("Error in " + name() + ". Epicenter should be same size as mesh dimension.");
     167             : 
     168          14 :     _rupture_speed = getParam<Real>("rupture_speed");
     169          14 :     _epicenter = getParam<std::vector<Real>>("epicenter");
     170             : 
     171           7 :     if (_rupture_speed <= 0.0)
     172           2 :       mooseError("Error in " + name() + ". Rupture speed has to be positive.");
     173           6 :   }
     174          68 :   else if ((isParamValid("rupture_speed") & !isParamValid("epicenter")) ||
     175          81 :            (!isParamValid("rupture_speed") & isParamValid("epicenter")))
     176           4 :     mooseError("Error in " + name() +
     177             :                ". Either both rupture speed and epicenter should be "
     178             :                "provided or neither should be provided.");
     179          21 : }
     180             : 
     181             : void
     182        2132 : SeismicSource::addPoints()
     183             : {
     184        4264 :   if (isParamValid("point"))
     185             :   {
     186         270 :     std::vector<Real> point_param = getParam<std::vector<Real>>("point");
     187          90 :     _source_location(0) = point_param[0];
     188          90 :     if (point_param.size() > 1)
     189             :     {
     190          90 :       _source_location(1) = point_param[1];
     191          90 :       if (point_param.size() > 2)
     192          90 :         _source_location(2) = point_param[2];
     193             :     }
     194          90 :     addPoint(_source_location);
     195          90 :   }
     196        4084 :   else if (isParamValid("position_function"))
     197             :   {
     198        4084 :     unsigned int number = getParam<unsigned int>("number");
     199             :     std::vector<FunctionName> position_function =
     200        4084 :         getParam<std::vector<FunctionName>>("position_function");
     201        2042 :     std::vector<const Function *> pos_function(_mesh.dimension());
     202             : 
     203        6126 :     for (unsigned int i = 0; i < _mesh.dimension(); ++i)
     204        4084 :       pos_function[i] = &getFunctionByName(position_function[i]);
     205             : 
     206       16336 :     for (unsigned int i = 0; i < number; ++i)
     207             :     {
     208       42882 :       for (unsigned int j = 0; j < _mesh.dimension(); ++j)
     209       28588 :         _source_location(j) = (pos_function[j])->value(i + 1, _qp);
     210             : 
     211       14294 :       addPoint(_source_location, i);
     212             :     }
     213        2042 :   }
     214        2132 : }
     215             : 
     216             : Real
     217       35424 : SeismicSource::computeQpResidual()
     218             : {
     219             :   /**
     220             :    * The time at which each point source ruptures is given by time_shift.
     221             :    * The time_shift parameter is calculated based on rupture speed and the distance
     222             :    * of the source from the epicenter.
     223             :    **/
     224             : 
     225             :   Real time_shift = 0.0;
     226       35424 :   if (_rupture_speed != std::numeric_limits<unsigned int>::max())
     227             :   {
     228             :     Real distance = 0.0;
     229       22960 :     if (_mesh.dimension() == 1)
     230           0 :       distance = std::abs(_physical_point[_qp](0) - _epicenter[0]);
     231       22960 :     else if (_mesh.dimension() == 2)
     232       22960 :       distance = std::sqrt(pow(_physical_point[_qp](0) - _epicenter[0], 2.0) +
     233       22960 :                            pow(_physical_point[_qp](1) - _epicenter[1], 2.0));
     234           0 :     else if (_mesh.dimension() == 3)
     235           0 :       distance = std::sqrt(pow(_physical_point[_qp](0) - _epicenter[0], 2.0) +
     236           0 :                            pow(_physical_point[_qp](1) - _epicenter[1], 2.0) +
     237           0 :                            pow(_physical_point[_qp](2) - _epicenter[2], 2.0));
     238             : 
     239       22960 :     time_shift = distance / _rupture_speed;
     240             :   }
     241             : 
     242             :   /**
     243             :    *  f_i(x,t) = - d M_ij(x,t)/ d x_j (summation over index j)
     244             :    *  For a point source applied at p, M_ij(x,t) = M_ij(t)*delta(x-p)
     245             :    *  Therefore f_i(x,t) = - M_ij(t)*d delta(x-p)/ d x_j
     246             :    *  This in weak form simplifies to the volume integral of M_ij(t)* delta(x-p)
     247             :    ** d test_i / d x_j
     248             :    **/
     249       35424 :   if (_t > time_shift)
     250       27192 :     _moment_magnitude =
     251       27192 :         _shear_modulus * _area * _slip_function->value(_t + _alpha * _dt - time_shift, _qp);
     252             :   else
     253        8232 :     _moment_magnitude = 0.0;
     254             : 
     255       35424 :   _force = 0.0;
     256      106752 :   for (unsigned int i = 0; i < _mesh.dimension(); i++)
     257       71328 :     _force += _moment[_component][i] * _grad_test[_i][_qp](i);
     258             : 
     259       35424 :   return -_force * _moment_magnitude;
     260             : }