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 "UpdatedLagrangianStressDivergence.h"
      11             : 
      12             : registerMooseObject("SolidMechanicsApp", UpdatedLagrangianStressDivergence);
      13             : 
      14             : template <class G>
      15         728 : UpdatedLagrangianStressDivergenceBase<G>::UpdatedLagrangianStressDivergenceBase(
      16             :     const InputParameters & parameters)
      17             :   : LagrangianStressDivergenceBase(parameters),
      18         728 :     _stress(getMaterialPropertyByName<RankTwoTensor>(_base_name + "cauchy_stress")),
      19        1456 :     _material_jacobian(getMaterialPropertyByName<RankFourTensor>(_base_name + "cauchy_jacobian")),
      20             : 
      21             :     // Assembly quantities in the reference frame for stabilization
      22         728 :     _assembly_undisplaced(_fe_problem.assembly(_tid, this->_sys.number())),
      23         728 :     _grad_phi_undisplaced(_assembly_undisplaced.gradPhi()),
      24         728 :     _JxW_undisplaced(_assembly_undisplaced.JxW()),
      25         728 :     _coord_undisplaced(_assembly_undisplaced.coordTransformation()),
      26         728 :     _q_point_undisplaced(_assembly_undisplaced.qPoints())
      27             : {
      28             :   // This kernel requires used_displaced_mesh to be true if large kinematics
      29             :   // is on
      30        1742 :   if (_large_kinematics && (!getParam<bool>("use_displaced_mesh")))
      31           0 :     mooseError("The UpdatedLagrangianStressDivergence kernels requires "
      32             :                "used_displaced_mesh = true for large_kinematics = true");
      33             : 
      34             :   // Similarly, if large kinematics is off so should use_displaced_mesh
      35        1508 :   if (!_large_kinematics && (getParam<bool>("use_displaced_mesh")))
      36           0 :     mooseError("The UpdatedLagrangianStressDivergence kernels requires "
      37             :                "used_displaced_mesh = false for large_kinematics = false");
      38             : 
      39             :   // TODO: add weak plane stress support
      40         728 :   if (_out_of_plane_strain)
      41           0 :     mooseError("The UpdatedLagrangianStressDivergence kernels do not yet support the weak plane "
      42             :                "stress formulation. Please use the TotalLagrangianStressDivergecen kernels.");
      43         728 : }
      44             : 
      45             : template <class G>
      46             : RankTwoTensor
      47   498518272 : UpdatedLagrangianStressDivergenceBase<G>::gradTest(unsigned int component)
      48             : {
      49             :   // F-bar doesn't modify the test function
      50   498518272 :   return G::gradOp(component, _grad_test[_i][_qp], _test[_i][_qp], _q_point[_qp]);
      51             : }
      52             : 
      53             : template <class G>
      54             : RankTwoTensor
      55   256109760 : UpdatedLagrangianStressDivergenceBase<G>::gradTrial(unsigned int component)
      56             : {
      57   256109760 :   return _stabilize_strain ? gradTrialStabilized(component) : gradTrialUnstabilized(component);
      58             : }
      59             : 
      60             : template <class G>
      61             : RankTwoTensor
      62   398456704 : UpdatedLagrangianStressDivergenceBase<G>::gradTrialUnstabilized(unsigned int component)
      63             : {
      64             :   // Without F-bar stabilization, simply return the gradient of the trial functions
      65   398456704 :   return G::gradOp(component, _grad_phi[_j][_qp], _phi[_j][_qp], _q_point[_qp]);
      66             : }
      67             : 
      68             : template <class G>
      69             : RankTwoTensor
      70    29799680 : UpdatedLagrangianStressDivergenceBase<G>::gradTrialStabilized(unsigned int component)
      71             : {
      72             :   // The base unstabilized trial function gradient
      73    29799680 :   const auto Gb = G::gradOp(component, _grad_phi[_j][_qp], _phi[_j][_qp], _q_point[_qp]);
      74             :   // The average trial function gradient
      75    29799680 :   const auto Ga = _avg_grad_trial[component][_j];
      76             :   // For updated Lagrangian, the modification is always linear
      77    29799680 :   return Gb + (Ga.trace() - Gb.trace()) / 3.0 * RankTwoTensor::Identity();
      78             : }
      79             : 
      80             : template <class G>
      81             : void
      82       62560 : UpdatedLagrangianStressDivergenceBase<G>::precalculateJacobianDisplacement(unsigned int component)
      83             : {
      84             :   // For updated Lagrangian, the averaging is taken on the reference frame. If large kinematics is
      85             :   // used, the averaged gradients should be pushed forward to the current frame.
      86      542880 :   for (auto j : make_range(_phi.size()))
      87             :   {
      88      480320 :     _avg_grad_trial[component][j] = StabilizationUtils::elementAverage(
      89     4244160 :         [this, component, j](unsigned int qp)
      90             :         {
      91             :           return G::gradOp(
      92     3763840 :               component, _grad_phi_undisplaced[j][qp], _phi[j][qp], _q_point_undisplaced[qp]);
      93             :         },
      94             :         _JxW_undisplaced,
      95             :         _coord_undisplaced);
      96      480320 :     if (_large_kinematics)
      97             :       // Push forward to the current frame.
      98             :       // The average deformation gradient is the same at all qps.
      99       51104 :       _avg_grad_trial[component][j] *= _F_avg[0].inverse();
     100             :   }
     101       62560 : }
     102             : 
     103             : template <class G>
     104             : Real
     105   242211904 : UpdatedLagrangianStressDivergenceBase<G>::computeQpResidual()
     106             : {
     107   242211904 :   return gradTest(_alpha).doubleContraction(_stress[_qp]);
     108             : }
     109             : 
     110             : template <class G>
     111             : Real
     112   256109760 : UpdatedLagrangianStressDivergenceBase<G>::computeQpJacobianDisplacement(unsigned int alpha,
     113             :                                                                         unsigned int beta)
     114             : {
     115   256109760 :   const auto grad_test = gradTest(alpha);
     116   256109760 :   const auto grad_trial = gradTrial(beta);
     117             : 
     118             :   //           J^{alpha beta} = J^{alpha beta}_material + J^{alpha beta}_geometric
     119             :   //  J^{alpha beta}_material = phi^alpha_{i, j} T_{ijkl} f^{-1}_{km} g^beta_{ml}
     120             :   // J^{alpha beta}_geometric = sigma_{ij} (phi^alpha_{k, k} psi^beta_{i, j} -
     121             :   //                                        phi^alpha_{k, j} psi^beta_{i, k})
     122             : 
     123             :   // The material jacobian
     124   256109760 :   Real J = grad_test.doubleContraction(_material_jacobian[_qp] * (_f_inv[_qp] * grad_trial));
     125             : 
     126             :   // The geometric jacobian
     127   256109760 :   if (_large_kinematics)
     128             :   {
     129             :     // No stablization in the geometric jacobian
     130   172146624 :     const auto grad_trial_ust = gradTrialUnstabilized(beta);
     131   344293248 :     J += _stress[_qp].doubleContraction(grad_test) * grad_trial_ust.trace() -
     132   172146624 :          _stress[_qp].doubleContraction(grad_test * grad_trial_ust);
     133             :   }
     134             : 
     135   256109760 :   return J;
     136             : }
     137             : 
     138             : template <class G>
     139             : Real
     140      196608 : UpdatedLagrangianStressDivergenceBase<G>::computeQpJacobianTemperature(unsigned int cvar)
     141             : {
     142             :   // Multiple eigenstrains may depend on the same coupled var
     143      196608 :   RankTwoTensor total_deigen;
     144      393216 :   for (const auto deigen_darg : _deigenstrain_dargs[cvar])
     145      196608 :     total_deigen += (*deigen_darg)[_qp];
     146             : 
     147      196608 :   RankFourTensor C = _material_jacobian[_qp];
     148      196608 :   RankFourTensor Csym = 0.5 * (C + C.transposeMajor().transposeIj().transposeMajor());
     149             : 
     150      196608 :   return -(Csym * total_deigen).doubleContraction(gradTest(_alpha)) * _temperature->phi()[_j][_qp];
     151             : }
     152             : 
     153             : template class UpdatedLagrangianStressDivergenceBase<GradientOperatorCartesian>;