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         594 : UpdatedLagrangianStressDivergenceBase<G>::UpdatedLagrangianStressDivergenceBase(
      16             :     const InputParameters & parameters)
      17             :   : LagrangianStressDivergenceBase(parameters),
      18         594 :     _stress(getMaterialPropertyByName<RankTwoTensor>(_base_name + "cauchy_stress")),
      19        1188 :     _material_jacobian(getMaterialPropertyByName<RankFourTensor>(_base_name + "cauchy_jacobian")),
      20             : 
      21             :     // Assembly quantities in the reference frame for stabilization
      22         594 :     _assembly_undisplaced(_fe_problem.assembly(_tid, this->_sys.number())),
      23         594 :     _grad_phi_undisplaced(_assembly_undisplaced.gradPhi()),
      24         594 :     _JxW_undisplaced(_assembly_undisplaced.JxW()),
      25         594 :     _coord_undisplaced(_assembly_undisplaced.coordTransformation()),
      26         594 :     _q_point_undisplaced(_assembly_undisplaced.qPoints())
      27             : {
      28             :   // The use_displaced_mesh <-> large_kinematics consistency check lives in initialSetup(), since
      29             :   // large_kinematics is now derived from the strain calculator's guarantee (not available in the
      30             :   // constructor).
      31             : 
      32             :   // TODO: add weak plane stress support
      33         594 :   if (_out_of_plane_strain)
      34           0 :     mooseError("The UpdatedLagrangianStressDivergence kernels do not yet support the weak plane "
      35             :                "stress formulation. Please use the TotalLagrangianStressDivergecen kernels.");
      36             : 
      37             :   // The UL push-forward in `precalculateJacobianDisplacement` uses `_F_avg` (= avg(F_ust))
      38             :   // for the spatial-to-reference frame mapping; combining it with the incremental F-bar
      39             :   // chain (which expects `_avg_grad_trial` to represent deltaf_avg, not deltaF_avg) hasn't been
      40             :   // worked out yet. Use Total Lagrangian for incremental F-bar.
      41         594 :   if (_stabilize_strain && _F_bar_mode == FBarMode::Incremental)
      42           0 :     mooseError("`F_bar_mode = incremental` is not yet supported with the UpdatedLagrangian "
      43             :                "kernels; use TotalLagrangianStressDivergence (or `F_bar_mode = total`).");
      44         594 : }
      45             : 
      46             : template <class G>
      47             : RankTwoTensor
      48   464359424 : UpdatedLagrangianStressDivergenceBase<G>::gradTest(unsigned int component)
      49             : {
      50             :   // F-bar doesn't modify the test function
      51   464359424 :   return G::gradOp(component, _grad_test[_i][_qp], _test[_i][_qp], _q_point[_qp]);
      52             : }
      53             : 
      54             : template <class G>
      55             : RankTwoTensor
      56           0 : UpdatedLagrangianStressDivergenceBase<G>::gradTrial(unsigned int component)
      57             : {
      58             :   // F-bar stabilization is handled explicitly in computeQpJacobianDisplacement via the stored
      59             :   // partial derivatives _d_F_stab_d_F_ust and _d_F_stab_d_F_avg, so gradTrial always returns
      60             :   // the unstabilized spatial gradient.
      61           0 :   return gradTrialUnstabilized(component);
      62             : }
      63             : 
      64             : template <class G>
      65             : RankTwoTensor
      66   188517888 : UpdatedLagrangianStressDivergenceBase<G>::gradTrialUnstabilized(unsigned int component)
      67             : {
      68             :   // Without F-bar stabilization, simply return the gradient of the trial functions
      69   188517888 :   return G::gradOp(component, _grad_phi[_j][_qp], _phi[_j][_qp], _q_point[_qp]);
      70             : }
      71             : 
      72             : template <class G>
      73             : void
      74       61208 : UpdatedLagrangianStressDivergenceBase<G>::precalculateJacobianDisplacement(unsigned int component)
      75             : {
      76             :   // For updated Lagrangian, the averaging is taken on the reference frame. If large kinematics is
      77             :   // used, the averaged gradients should be pushed forward to the current frame.
      78      536760 :   for (auto j : make_range(_phi.size()))
      79             :   {
      80      475552 :     _avg_grad_trial[component][j] = StabilizationUtils::elementAverage(
      81      475552 :         [this, component, j](unsigned int qp)
      82             :         {
      83     3748928 :           return G::gradOp(
      84     3748928 :               component, _grad_phi_undisplaced[j][qp], _phi[j][qp], _q_point_undisplaced[qp]);
      85             :         },
      86             :         _JxW_undisplaced,
      87             :         _coord_undisplaced);
      88      475552 :     if (_large_kinematics)
      89             :       // Push forward to the current frame.
      90             :       // The average deformation gradient is the same at all qps.
      91      130000 :       _avg_grad_trial[component][j] *= _F_avg[0].inverse();
      92             :   }
      93       61208 : }
      94             : 
      95             : template <class G>
      96             : Real
      97   275595776 : UpdatedLagrangianStressDivergenceBase<G>::computeQpResidual()
      98             : {
      99   275595776 :   return gradTest(_alpha).doubleContraction(_stress[_qp]);
     100             : }
     101             : 
     102             : template <class G>
     103             : Real
     104   188517888 : UpdatedLagrangianStressDivergenceBase<G>::computeQpJacobianDisplacement(unsigned int alpha,
     105             :                                                                         unsigned int beta)
     106             : {
     107   188517888 :   const auto grad_test = gradTest(alpha);
     108   188517888 :   const auto grad_trial = gradTrialUnstabilized(beta);
     109             : 
     110             :   //           J^{alpha beta} = J^{alpha beta}_material + J^{alpha beta}_geometric
     111             :   //  J^{alpha beta}_material = phi^alpha : T : d(dL)/d(grad u) * grad_trial
     112             :   // J^{alpha beta}_geometric = sigma_{ij} (phi^alpha_{k, k} psi^beta_{i, j} -
     113             :   //                                        phi^alpha_{k, j} psi^beta_{i, k})
     114             : 
     115             :   // Local contribution to delta(F_ust): pull the spatial trial gradient back through the
     116             :   // literal n+1 deformation gradient (use_displaced_mesh = true uses F_actual, regardless
     117             :   // of alpha or F-bar), then chain to F_ust via _d_F_d_grad_u (= alpha * I^(4)).
     118   188517888 :   const RankTwoTensor delta_grad_u_local =
     119   188517888 :       _large_kinematics ? grad_trial * _F_actual[_qp] : grad_trial;
     120   188517888 :   const RankTwoTensor delta_F_ust_local = _d_F_d_grad_u[_qp] * delta_grad_u_local;
     121             : 
     122             :   // Non-local F-bar contribution to delta(F_avg): _avg_grad_trial is stored as
     123             :   // avg(gradTrial_reference) * F_avg^{-1} (see precalculateJacobianDisplacement). Multiply
     124             :   // by F_avg to recover the reference-frame averaged gradient, then chain to F_ust the
     125             :   // same way.
     126   188517888 :   RankTwoTensor delta_F_avg;
     127   188517888 :   if (_stabilize_strain)
     128             :   {
     129    29771392 :     const RankTwoTensor delta_grad_u_avg =
     130    29771392 :         _large_kinematics ? _avg_grad_trial[beta][_j] * _F_avg[_qp] : _avg_grad_trial[beta][_j];
     131    29771392 :     delta_F_avg = _d_F_d_grad_u[_qp] * delta_grad_u_avg;
     132             :   }
     133             : 
     134             :   // Stabilized delta(F_stab) through the F-bar tangent. With F-bar off, the partials are
     135             :   // _d_F_stab_d_F_ust = I^(4) and _d_F_stab_d_F_avg = 0, so this reduces to delta_F_ust_local.
     136   188517888 :   const RankTwoTensor delta_F_stab =
     137   188517888 :       _d_F_stab_d_F_ust[_qp] * delta_F_ust_local + _d_F_stab_d_F_avg[_qp] * delta_F_avg;
     138             : 
     139   188517888 :   const RankTwoTensor delta_dL = _d_deformation_gradient_increment_d_F[_qp] * delta_F_stab;
     140             : 
     141             :   // The material jacobian
     142   188517888 :   Real J = grad_test.doubleContraction(_material_jacobian[_qp] * delta_dL);
     143             : 
     144             :   // The geometric jacobian (no F-bar stabilization in this term)
     145   188517888 :   if (_large_kinematics)
     146             :   {
     147   239116032 :     J += _stress[_qp].doubleContraction(grad_test) * grad_trial.trace() -
     148   119558016 :          _stress[_qp].doubleContraction(grad_test * grad_trial);
     149             :   }
     150             : 
     151   188517888 :   return J;
     152             : }
     153             : 
     154             : template <class G>
     155             : Real
     156      245760 : UpdatedLagrangianStressDivergenceBase<G>::computeQpJacobianTemperature(unsigned int cvar)
     157             : {
     158             :   // Multiple eigenstrains may depend on the same coupled var
     159      245760 :   RankTwoTensor total_deigen;
     160      491520 :   for (const auto deigen_darg : _deigenstrain_dargs[cvar])
     161      245760 :     total_deigen += (*deigen_darg)[_qp];
     162             : 
     163      245760 :   RankFourTensor C = _material_jacobian[_qp];
     164      245760 :   RankFourTensor Csym = 0.5 * (C + C.transposeMajor().transposeIj().transposeMajor());
     165             : 
     166      245760 :   return -(Csym * total_deigen).doubleContraction(gradTest(_alpha)) * _temperature->phi()[_j][_qp];
     167             : }
     168             : 
     169             : template class UpdatedLagrangianStressDivergenceBase<GradientOperatorCartesian>;