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 "ComputeFiniteBeamStrain.h"
      11             : #include "Assembly.h"
      12             : #include "NonlinearSystem.h"
      13             : #include "MooseVariable.h"
      14             : 
      15             : #include "libmesh/quadrature.h"
      16             : #include "libmesh/utility.h"
      17             : 
      18             : registerMooseObject("SolidMechanicsApp", ComputeFiniteBeamStrain);
      19             : 
      20             : InputParameters
      21          72 : ComputeFiniteBeamStrain::validParams()
      22             : {
      23          72 :   InputParameters params = ComputeIncrementalBeamStrain::validParams();
      24          72 :   params.addClassDescription("Compute a rotation increment for finite rotations of the beam and "
      25             :                              "computes the small/large strain increments in the current rotated "
      26             :                              "configuration of the beam.");
      27          72 :   return params;
      28           0 : }
      29             : 
      30          54 : ComputeFiniteBeamStrain::ComputeFiniteBeamStrain(const InputParameters & parameters)
      31             :   : ComputeIncrementalBeamStrain(parameters),
      32         108 :     _total_rotation_old(getMaterialPropertyOld<RankTwoTensor>("total_rotation"))
      33             : {
      34          54 : }
      35             : 
      36             : void
      37       88088 : ComputeFiniteBeamStrain::computeRotation()
      38             : {
      39             :   // First - convert the differential beam displacement to the beam configuration at previous time
      40             :   // step
      41       88088 :   const RealVectorValue delta_disp_local(_total_rotation_old[0] * (_disp1 - _disp0));
      42             : 
      43             :   // Second - calculate the incremental rotation matrix using Euler angles
      44             :   const Real intermediate_length_1 =
      45       88088 :       std::sqrt(Utility::pow<2>(_original_length[0] + delta_disp_local(0)) +
      46             :                 Utility::pow<2>(delta_disp_local(2)));
      47       88088 :   const Real cos_alpha = (_original_length[0] + delta_disp_local(0)) / intermediate_length_1;
      48       88088 :   const Real sin_alpha = std::sqrt(1.0 - Utility::pow<2>(cos_alpha));
      49             : 
      50             :   const Real intermediate_length_2 =
      51       88088 :       std::sqrt(Utility::pow<2>(intermediate_length_1) + Utility::pow<2>(delta_disp_local(1)));
      52       88088 :   const Real sin_beta = delta_disp_local(1) / intermediate_length_2;
      53       88088 :   const Real cos_beta = std::sqrt(1.0 - Utility::pow<2>(sin_beta));
      54             : 
      55       88088 :   const RealVectorValue rotation_d_1(cos_alpha * cos_beta, sin_beta, sin_alpha * cos_beta);
      56       88088 :   const RealVectorValue rotation_d_2(-cos_alpha * sin_beta, cos_beta, -sin_alpha * sin_beta);
      57             :   const RealVectorValue rotation_d_3(-sin_alpha, 0.0, cos_alpha);
      58             : 
      59             :   const auto rotation_d =
      60       88088 :       RankTwoTensor ::initializeFromRows(rotation_d_1, rotation_d_2, rotation_d_3);
      61             : 
      62             :   // Convert average rotational displacement to the beam configuration at previous time step
      63       88088 :   const RealVectorValue avg_rot_local(_total_rotation_old[0] * (_rot0 + _rot1));
      64             : 
      65             :   const Real gamma_increment =
      66       88088 :       0.5 * (rotation_d_1(0) * avg_rot_local(0) + rotation_d_1(1) * avg_rot_local(1) +
      67       88088 :              rotation_d_1(2) * avg_rot_local(2));
      68             : 
      69       88088 :   RankTwoTensor rotation_a;
      70       88088 :   rotation_a(0, 0) = 1.0;
      71       88088 :   rotation_a(1, 1) = std::cos(gamma_increment);
      72       88088 :   rotation_a(1, 2) = std::sin(gamma_increment);
      73       88088 :   rotation_a(2, 1) = -rotation_a(1, 2);
      74       88088 :   rotation_a(2, 2) = rotation_a(1, 1);
      75             : 
      76             :   // Total rotation matrix from global configuration to beam local config at current time
      77       88088 :   _total_rotation[0] = rotation_a * rotation_d * _total_rotation_old[0];
      78       88088 : }