Merge pull request #28742 from Kavan-Khaledi/shell_local_coordinates

Calculating local quantities for Q4 shell elements

modules/solid_mechanics/doc/content/modules/solid_mechanics/ShellElements.md

@@ -39,7 +39,7 @@ If $^t \mathbf{g_i} = \partial ^t \mathbf{x}/\partial r_i$ are the covariant bas
 \end{equation}
 
 The expressions for the small and large strains in the 11, 22, 12, 13 and 23 directions are provided in equations 21 to 24 of [!cite](dvorkin1984continuum). The expressions for the shear strains in the 13 and 23 directions also include a correction for shear locking. In MOOSE the small strain increments are computed in
-[ADComputeIncrementalShellStrain](/ADComputeIncrementalShellStrain.md) and the small + large strain increments are computed in [ADComputeFiniteShellStrain](/ADComputeFiniteShellStrain.md)). Note that strain in the 33 direction is not computed.
+[ADComputeIncrementalShellStrain](/ADComputeIncrementalShellStrain.md) and the small + large strain increments are computed in [ADComputeFiniteShellStrain](/ADComputeFiniteShellStrain.md). Note that strain in the 33 direction is not computed.
 
 Apart from the strain increment, two other matrices (B and BNL) are computed by the strain class. These two matrices connect the nodal displacements/rotations to the small and large strains, respectively, i.e., $e=B u$ and $\eta = B_{NL} u$. These two matrices are then used to compute the residuals in [ADStressDivergenceShell](/ADStressDivergenceShell.md).
 
@@ -63,13 +63,13 @@ g^2=\dfrac{g_3 \times g_1}{g_1 \cdot (g_2 \times g_3)}
 \begin{equation}
 g^3=\dfrac{g_1 \times g_2}{g_1 \cdot (g_2 \times g_3)}
 \end{equation}
-These contravariant base vectors satisfy the orthogonality condition: $g^i \cdot g_j =\delta_{ij}$, where $\delta_{ij}$ is Kronecker delta. The $\hat{e}_i$ vectors are the local cartesian coordinate system computed using the covariant base vectors, i.e.,
+These contravariant base vectors satisfy the orthogonality condition: $g^i \cdot g_j =\delta_{ij}$, where $\delta_{ij}$ is Kronecker delta. The $\hat{e}_i$ vectors are the local cartesian coordinate system defined based on a user-defined reference first local direction vector $\hat{e}_{1r}$, which is defined by the `reference_first_local_direction` parameter, which defaults to the $X$ direction (i.e. $\hat{e}_{1r}=(1, 0, 0)$). The projection of $\hat{e}_{1r}$ onto the local element's plane provides the first local direction, so that the $\hat{e}$ vectors are computed as:
 
 \begin{equation}
-\hat{e}_3 = \frac{g_3}{|g_3|}; \hat{e}_1 = \frac{g_2 \times \hat{e}^3}{|g_2 \times \hat{e}^3|}; \hat{e}_2 = \hat{e}_3 \times \hat{e}_1
+\hat{e}_3 = V_n; \hat{e}_2 = \frac{\hat{e}_3 \times \hat{e}_{1r}}{|\hat{e}_3 \times \hat{e}_{1r}|}; \hat{e}_1 = \hat{e}_2 \times \hat{e}_3
 \end{equation}
 
-The elasticity tensor at each quadrature point (qp) is computed in [ADComputeIsotropicElasticityTensorShell](/ADComputeIsotropicElasticityTensorShell.md).
+In this equation, $V_n$ represents the element's normal vector. The second local vector $e_2$ is obtained via the cross product of $e_3$ and $\hat{e}_{1r}$. The projection of $\hat{e}_{1r}$ onto each element ensures that $\hat{e}_1$ lies in-plane and is perpendicular to the element's normal. The projection of $\hat{e}_{1r}$ on the shell element is then determined via the cross product of $e_2$ and $e_3$. The elasticity tensor at each quadrature point (qp) is computed in [ADComputeIsotropicElasticityTensorShell](/ADComputeIsotropicElasticityTensorShell.md).
 
 ## Stress calculation
 

modules/solid_mechanics/doc/content/source/auxkernels/ShellLocalCoordinatesAux.md

@@ -0,0 +1,18 @@
+# ShellLocalCoordinatesAux
+
+!syntax description /AuxKernels/ShellLocalCoordinatesAux
+
+!alert note
+The three Cartesian local vectors for each shell element are indexed as follows: the first vector is indexed by 0, the second vector by 1, and the normal vector by 2. The convention used to define the direction of these vectors is explained in [ShellElements](/ShellElements.md)
+
+
+
+## Example Input Syntax
+
+!listing modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_local_stress.i block=AuxKernels/first_axis_x
+
+!syntax parameters /AuxKernels/ShellLocalCoordinatesAux
+
+!syntax inputs /AuxKernels/ShellLocalCoordinatesAux
+
+!syntax children /AuxKernels/ShellLocalCoordinatesAux

modules/solid_mechanics/doc/content/source/auxkernels/ShellResultantsAux.md

@@ -0,0 +1,65 @@
+# ShellResultantsAux
+
+!syntax description /AuxKernels/ShellResultantsAux
+
+!alert note
+The three Cartesian local vectors for each shell element are indexed as follows: the first vector is indexed by 0, the second vector by 1, and the normal vector by 2. The convention used to define the direction of these vectors is explained in [ShellElements](/ShellElements.md)
+
+The following stress resultants are computed using this auxiliary kernel:
+
+axial_force_0: The in-plane axial force in the direction of the first local coordinate axis:
+
+  $$F_{0} = \int_{-t/2}^{t/2} \sigma_{00} dz$$
+
+
+axial_force_1: The in-plane axial force in the direction of the second local coordinate axis:
+
+ $$ F_{1} = \int_{-t/2}^{t/2} \sigma_{11} dz$$
+
+
+normal_force: The normal force applied to the shell element in the thickness direction. This force is expected to be always zero due to the plane stress assumption used in the shell element formulation, which disregards out-of-plane stresses:
+
+  $$F_{N} = \int_{-t/2}^{t/2} \sigma_{22} dz$$
+
+
+bending_moments_0: The bending moment about the first local coordinate axis:
+
+  $$M_{0} = \int_{-t/2}^{t/2} \sigma_{11} z dz$$
+
+
+bending_moment_1: The bending moment about the second local coordinate axis:
+
+ $$ M_{1} = \int_{-t/2}^{t/2} \sigma_{00} z dz$$
+
+
+bending_moment_01: The in-plane bending moment:
+
+$$  M_{01} =M_{10}= \int_{-t/2}^{t/2} \sigma_{01} z dz$$
+
+
+shear_force_01: The in-plane shear force:
+
+ $$ Q_{01} =Q_{10}= \int_{-t/2}^{t/2} \sigma_{01} dz$$
+
+
+shear_force_02: The transverse shear force:
+
+$$Q_{02} =Q_{20}= \int_{-t/2}^{t/2} \sigma_{02} dz$$
+
+
+shear_force_12: The transverse shear force:
+
+$$  Q_{12} =Q_{21}= \int_{-t/2}^{t/2} \sigma_{12} dz$$
+
+
+
+
+## Example Input Syntax
+
+!listing modules/solid_mechanics/test/tests/shell/static/plate_cantilever.i block=AuxKernels/moment_22
+
+!syntax parameters /AuxKernels/ShellResultantsAux
+
+!syntax inputs /AuxKernels/ShellResultantsAux
+
+!syntax children /AuxKernels/ShellResultantsAux

modules/solid_mechanics/include/auxkernels/ShellLocalCoordinatesAux.h

@@ -0,0 +1,39 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "AuxKernel.h"
+#include "RankTwoTensor.h"
+#include "MooseEnum.h"
+
+class ShellLocalCoordinatesAux : public AuxKernel
+{
+public:
+  static InputParameters validParams();
+
+  ShellLocalCoordinatesAux(const InputParameters & parameters);
+
+protected:
+  virtual Real computeValue() override;
+
+  /// Base name of the material system used to calculate the elastic energy
+  const std::string _base_name;
+
+  enum class PropertyType
+  {
+    first_local_vector,
+    second_local_vector,
+    normal_local_vector
+  } _property;
+  unsigned int _component;
+
+  /// The local stress tensor
+  const MaterialProperty<RankTwoTensor> * _local_coordinates;
+};

modules/solid_mechanics/include/auxkernels/ShellResultantsAux.h

@@ -0,0 +1,59 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "AuxKernel.h"
+#include "RankTwoTensor.h"
+#include "MooseEnum.h"
+
+// Forward declarations
+namespace libMesh
+{
+class QGauss;
+}
+
+class ShellResultantsAux : public AuxKernel
+{
+public:
+  static InputParameters validParams();
+
+  ShellResultantsAux(const InputParameters & parameters);
+
+protected:
+  virtual Real computeValue();
+
+  /// Base name of the material system used to calculate the elastic energy
+  const std::string _base_name;
+
+  /// Coupled variable for the shell thickness
+  const VariableValue & _thickness;
+
+  enum class ResultantType
+  {
+    axial_force_0,
+    axial_force_1,
+    normal_force,
+    bending_moment_0,
+    bending_moment_1,
+    bending_moment_01,
+    shear_force_01,
+    shear_force_02,
+    shear_force_12
+  } _resultant;
+  /// Quadrature rule in the out of plane direction
+  std::unique_ptr<QGauss> _t_qrule;
+
+  /// Quadrature points and weights in the out of plane direction in isoparametric coordinate system
+  std::vector<Point> _t_points;
+  std::vector<Real> _t_weights;
+
+  /// The local stress tensor
+  std::vector<const MaterialProperty<RankTwoTensor> *> _local_stress_t_points;
+};

modules/solid_mechanics/include/materials/ADComputeIncrementalShellStrain.h

@@ -47,6 +47,9 @@ class ADComputeIncrementalShellStrain : public Material
   /// Computes the transformation matrix from natural coordinates to local cartesian coordinates for elasticity tensor transformation
   virtual void computeGMatrix();
 
+  /// Computes the element's local coordinates and store in _e1, _e2 and _e3
+  virtual void computeLocalCoordinates();
+
   /// Number of coupled rotational variables
   unsigned int _nrot;
 
@@ -158,6 +161,10 @@ class ADComputeIncrementalShellStrain : public Material
   /// Old material property containing jacobian of transformation
   std::vector<const MaterialProperty<Real> *> _J_map_old;
 
+  /// Transformation matrix to map stresses from global coordinate to the element's local coordinate
+  std::vector<MaterialProperty<RankTwoTensor> *> _local_transformation_matrix;
+  std::vector<const MaterialProperty<RankTwoTensor> *> _local_transformation_matrix_old;
+
   /// Covariant base vector matrix material property to transform stress
   std::vector<MaterialProperty<RankTwoTensor> *> _covariant_transformation_matrix;
   std::vector<const MaterialProperty<RankTwoTensor> *> _covariant_transformation_matrix_old;
@@ -177,6 +184,7 @@ class ADComputeIncrementalShellStrain : public Material
   ADRealVectorValue _x3;
   const NumericVector<Number> * const & _sol;
   const NumericVector<Number> & _sol_old;
+  const RealVectorValue _ref_first_local_dir;
   ADRealVectorValue _g3_a;
   ADRealVectorValue _g3_c;
   ADRealVectorValue _g3_b;
@@ -186,4 +194,7 @@ class ADComputeIncrementalShellStrain : public Material
   ADRealVectorValue _g2_b;
   ADRealVectorValue _g2_d;
   RankTwoTensor _unrotated_total_strain;
+  ADRealVectorValue _e1;
+  ADRealVectorValue _e2;
+  ADRealVectorValue _e3;
 };

modules/solid_mechanics/include/materials/ADComputeShellStress.h

@@ -47,12 +47,18 @@ class ADComputeShellStress : public Material
   /// Quadrature points along thickness
   std::vector<Point> _t_points;
 
+  /// Transformation matrix to map the global stress to the element's local coordinate
+  std::vector<const MaterialProperty<RankTwoTensor> *> _local_transformation_matrix;
+
   /// Covariant base vector matrix material property to transform stress
   std::vector<const MaterialProperty<RankTwoTensor> *> _covariant_transformation_matrix;
 
   /// Global stress tensor material property
   std::vector<MaterialProperty<RankTwoTensor> *> _global_stress;
 
+  /// local stress tensor material property
+  std::vector<MaterialProperty<RankTwoTensor> *> _local_stress;
+
   /// Real value of stress in the local coordinate system
   RankTwoTensor _unrotated_stress;
 };

modules/solid_mechanics/src/auxkernels/ShellLocalCoordinatesAux.C

@@ -0,0 +1,71 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "ShellLocalCoordinatesAux.h"
+#include "libmesh/utility.h"
+#include "libmesh/string_to_enum.h"
+
+registerMooseObject("SolidMechanicsApp", ShellLocalCoordinatesAux);
+
+InputParameters
+ShellLocalCoordinatesAux::validParams()
+{
+  InputParameters params = AuxKernel::validParams();
+  params.addClassDescription(
+      "This AuxKernel stores a specific component of a shell element's local coordinate "
+      "vector in an auxiliary variable.");
+  params.addParam<std::string>("base_name", "Mechanical property base name");
+
+  MooseEnum property("first_local_vector second_local_vector normal_local_vector");
+  params.addRequiredParam<MooseEnum>(
+      "property",
+      property,
+      "The local axis to output: first_local_vector, second_local_vector or normal_local_vector");
+  params.addRequiredParam<unsigned int>(
+      "component", "The vector component of the local coordinate vector: 0, 1 or 2");
+
+  return params;
+}
+
+ShellLocalCoordinatesAux::ShellLocalCoordinatesAux(const InputParameters & parameters)
+  : AuxKernel(parameters),
+    _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
+    _property(getParam<MooseEnum>("property").getEnum<PropertyType>()),
+    _component(getParam<unsigned int>("component"))
+
+{
+  _local_coordinates =
+      &getMaterialProperty<RankTwoTensor>(_base_name + "local_transformation_t_points_0");
+
+  if (_component > 2)
+    mooseError("Invalid component: ",
+               _component,
+               ". The component index of a shell local vector must be 0, 1, or 2.");
+}
+
+Real
+ShellLocalCoordinatesAux::computeValue()
+{
+  Real output_value = 0.0;
+
+  switch (_property)
+  {
+    case PropertyType::first_local_vector:
+      output_value = (*_local_coordinates)[_qp](0, _component);
+      break;
+    case PropertyType::second_local_vector:
+      output_value = (*_local_coordinates)[_qp](1, _component);
+      break;
+    case PropertyType::normal_local_vector:
+      output_value = (*_local_coordinates)[_qp](2, _component);
+      break;
+  }
+
+  return output_value;
+}