Merge pull request #24097 from jessecarterMOOSE/rotate_umat

Rotate UMAT

modules/solid_mechanics/doc/content/source/materials/abaqus/AbaqusUMATStress.md

@@ -96,6 +96,14 @@ An example of how to pass the user object in an input file is given below:
 Note that the step capability is three-pronged: 1) It allows to pass the step number to the UMAT
 routine via the present object, 2) It allows to pass the step number to the [AbaqusUExternalDB](/AbaqusUExternalDB.md) plugin, and 3) It allows to directly drive controls via step number in [StepPeriod](/StepPeriod.md).
 
+## UMAT Orientation
+
+Anisotropic material models typically have a specified local coordinate system that may need to be
+oriented differently than the global system. The `orientation` parameter takes a vector of 3 Euler
+angles that defines the rotation between the local and global coordinate systems. More information
+on the use of Euler angles for rotating material models in MOOSE can be found on the
+[ComputeElasticityTensor](/ComputeElasticityTensor.md) page.
+
 ## Example input file
 
 !listing modules/solid_mechanics/test/tests/umat/elastic_hardening/linear_strain_hardening.i block=Materials/constant

modules/solid_mechanics/include/materials/abaqus/AbaqusUMATStress.h

@@ -12,6 +12,7 @@
 #include "ComputeGeneralStressBase.h"
 #include "DynamicLibraryLoader.h"
 #include "ComputeFiniteStrain.h"
+#include "RotationTensor.h"
 #include "StepUOInterface.h"
 
 class StepUserObject;
@@ -223,6 +224,7 @@ class AbaqusUMATStress : public ComputeGeneralStressBase, public StepUOInterface
 
   // Time step rotation increment
   const OptionalMaterialProperty<RankTwoTensor> & _rotation_increment;
+  const OptionalMaterialProperty<RankTwoTensor> & _rotation_increment_old;
 
   // Coupled temperature field
   const VariableValue & _temperature;
@@ -246,6 +248,12 @@ class AbaqusUMATStress : public ComputeGeneralStressBase, public StepUOInterface
   /// parameter to assist with the transition to 1-based indexing
   const bool _use_one_based_indexing;
 
+  /// Rotation information
+  const bool _use_orientation;
+  const RotationTensor _R;
+  MaterialProperty<RankTwoTensor> & _total_rotation;
+  const MaterialProperty<RankTwoTensor> & _total_rotation_old;
+
 private:
   /// Method being used to compute strain and rotation increments
   const ComputeFiniteStrain::DecompMethod _decomposition_method;

modules/solid_mechanics/src/materials/ComputeIncrementalStrainBase.C

@@ -51,6 +51,7 @@ ComputeIncrementalStrainBase::initQpStatefulProperties()
   _mechanical_strain[_qp].zero();
   _total_strain[_qp].zero();
   _deformation_gradient[_qp].setToIdentity();
+  _rotation_increment[_qp].setToIdentity();
 }
 
 void

modules/solid_mechanics/src/materials/abaqus/AbaqusUMATStress.C

@@ -51,6 +51,9 @@ AbaqusUMATStress::validParams()
       "displaced mesh for computing displacements and quantities based on the deformed state.");
   params.addParam<UserObjectName>(
       "step_user_object", "The StepUserObject that provides times from simulation loading steps.");
+  params.addParam<RealVectorValue>(
+      "orientation",
+      "Euler angles that describe the orientation of the local material coordinate system.");
   return params;
 }
 
@@ -81,6 +84,8 @@ AbaqusUMATStress::AbaqusUMATStress(const InputParameters & parameters)
     _material_timestep(declareProperty<Real>(_base_name + "material_timestep_limit")),
     _rotation_increment(
         getOptionalMaterialProperty<RankTwoTensor>(_base_name + "rotation_increment")),
+    _rotation_increment_old(
+        getOptionalMaterialPropertyOld<RankTwoTensor>(_base_name + "rotation_increment")),
     _temperature(coupledValue("temperature")),
     _temperature_old(coupledValueOld("temperature")),
     _external_fields(isCoupled("external_fields") ? coupledValues("external_fields")
@@ -93,6 +98,10 @@ AbaqusUMATStress::AbaqusUMATStress(const InputParameters & parameters)
     _external_properties(_number_external_properties),
     _external_properties_old(_number_external_properties),
     _use_one_based_indexing(getParam<bool>("use_one_based_indexing")),
+    _use_orientation(isParamValid("orientation")),
+    _R(_use_orientation ? getParam<RealVectorValue>("orientation") : RealVectorValue(0.0)),
+    _total_rotation(declareProperty<RankTwoTensor>("total_rotation")),
+    _total_rotation_old(getMaterialPropertyOld<RankTwoTensor>("total_rotation")),
     _decomposition_method(
         getParam<MooseEnum>("decomposition_method").getEnum<ComputeFiniteStrain::DecompMethod>())
 {
@@ -158,6 +167,9 @@ AbaqusUMATStress::initQpStatefulProperties()
   _state_var[_qp].resize(_aqNSTATV);
   for (const auto i : make_range(_aqNSTATV))
     _state_var[_qp][i] = 0.0;
+
+  // Initialize total rotation tensor
+  _total_rotation[_qp] = _R;
 }
 
 void
@@ -201,11 +213,49 @@ AbaqusUMATStress::computeQpStress()
   // therefore, all unsymmetric matrices must be transposed before passing them to Fortran
   RankTwoTensor FBar_old_fortran = _Fbar_old[_qp].transpose();
   RankTwoTensor FBar_fortran = _Fbar[_qp].transpose();
-  RankTwoTensor DROT_fortran = _rotation_increment[_qp].transpose();
+
+  // DROT needed by UMAT will depend on kinematics and whether or not an intermediate configuration
+  // is used
+  RankTwoTensor DROT_fortran;
+  if (_use_orientation)
+  {
+    DROT_fortran = RankTwoTensor::Identity();
+  }
+  else
+  {
+    if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
+      DROT_fortran = _rotation_increment[_qp].transpose();
+    else
+      DROT_fortran = _rotation_increment_old[_qp].transpose();
+  }
+
   const Real * myDFGRD0 = &(FBar_old_fortran(0, 0));
   const Real * myDFGRD1 = &(FBar_fortran(0, 0));
   const Real * myDROT = &(DROT_fortran(0, 0));
 
+  // More local copies of materials so we can (optionally) rotate
+  RankTwoTensor stress_old = _stress_old[_qp];
+  RankTwoTensor total_strain_old = _total_strain_old[_qp];
+  RankTwoTensor strain_increment = _strain_increment[_qp];
+
+  // check if we need to rotate to intermediate configuration
+  if (_use_orientation)
+  {
+    // keep track of total rotation
+    _total_rotation[_qp] = _rotation_increment[_qp] * _total_rotation_old[_qp];
+    // rotate stress/strain/increment from reference configuration to intermediate configuration
+    stress_old.rotate(_total_rotation_old[_qp].transpose());
+    total_strain_old.rotate(_total_rotation_old[_qp].transpose());
+
+    if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
+      strain_increment.rotate(_total_rotation[_qp].transpose());
+    else
+      strain_increment.rotate(_total_rotation_old[_qp].transpose());
+  }
+  else if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
+    // rotate old stress to reference configuration
+    stress_old.rotate(_rotation_increment[_qp]);
+
   // copy because UMAT does not guarantee constness
   for (const auto i : make_range(9))
   {
@@ -225,18 +275,13 @@ AbaqusUMATStress::computeQpStress()
   static const std::array<std::pair<unsigned int, unsigned int>, 6> component{
       {{0, 0}, {1, 1}, {2, 2}, {0, 1}, {0, 2}, {1, 2}}};
 
-  // rotate old stress if HughesWinget
-  RankTwoTensor stress_old = _stress_old[_qp];
-  if (_decomposition_method == ComputeFiniteStrain::DecompMethod::HughesWinget)
-    stress_old.rotate(_rotation_increment[_qp]);
-
   for (const auto i : make_range(_aqNTENS))
   {
     const auto a = component[i].first;
     const auto b = component[i].second;
     _aqSTRESS[i] = stress_old(a, b);
-    _aqSTRAN[i] = _total_strain_old[_qp](a, b) * strain_factor[i];
-    _aqDSTRAN[i] = _strain_increment[_qp](a, b) * strain_factor[i];
+    _aqSTRAN[i] = total_strain_old(a, b) * strain_factor[i];
+    _aqDSTRAN[i] = strain_increment(a, b) * strain_factor[i];
   }
 
   // current coordinates
@@ -342,10 +387,6 @@ AbaqusUMATStress::computeQpStress()
   _stress[_qp] = RankTwoTensor(
       _aqSTRESS[0], _aqSTRESS[1], _aqSTRESS[2], _aqSTRESS[5], _aqSTRESS[4], _aqSTRESS[3]);
 
-  // Rotate the stress state to the current configuration, unless HughesWinget
-  if (_decomposition_method != ComputeFiniteStrain::DecompMethod::HughesWinget)
-    _stress[_qp].rotate(_rotation_increment[_qp]);
-
   // Build Jacobian matrix from UMAT's Voigt non-standard order to fourth order tensor.
   const unsigned int N = Moose::dim;
   const unsigned int ntens = N * (N + 1) / 2;
@@ -365,4 +406,14 @@ AbaqusUMATStress::computeQpStress()
                 k == l ? _aqDDSDDE[(nskip + i + j) * ntens + k]
                        : _aqDDSDDE[(nskip + i + j) * ntens + k + nskip + l];
         }
+
+  // check if we need to rotate from intermediate reference frame
+  if (_use_orientation)
+  {
+    // rotate to current configuration
+    _stress[_qp].rotate(_total_rotation[_qp]);
+    _jacobian_mult[_qp].rotate(_total_rotation[_qp]);
+  }
+  else if (_decomposition_method != ComputeFiniteStrain::DecompMethod::HughesWinget)
+    _stress[_qp].rotate(_rotation_increment[_qp]);
 }

modules/solid_mechanics/test/plugins/elastic_print.f

@@ -76,7 +76,12 @@ SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
 C     PRINTING FOR VERIFICATION PURPOSES
 C
       IF (TIME(1).GT.1.D0 .AND. COORDS(1).GT.0.25D0 .AND.
-     1   COORDS(2).GE.0.5D0 .AND. COORDS(3).GT.0.25D0) THEN
+     1   COORDS(2).GE.0.25D0 .AND. COORDS(3).GT.0.25D0) THEN
+
+        DO K1=1, NTENS
+          WRITE(*,110)  K1, STRESS(K1)
+        END DO
+
         DO K1=1, NTENS
           WRITE(*,120)  K1, STRAN(K1)
           WRITE(*,125)  K1, DSTRAN(K1)
@@ -121,6 +126,7 @@ SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
       ENDIF
 
 
+110   FORMAT ( 1X 'STRESS_', I2, :, 3X, 4F10.7 )
 120   FORMAT ( 1X 'STRAIN_', I2, :, 3X, 4F10.7 )
 125   FORMAT ( 1X 'DSTRAIN_', I2, :, 3X, 4F10.7 )
 130   FORMAT ( 1X 'COORDS_', I2, :, 3X, 4F10.7 )

modules/solid_mechanics/test/tests/umat/orient_umat/reference/shear_top_umat.inp

@@ -0,0 +1,71 @@
+** Reference Abaqus input for "hw_reference_umat_output" MOOSE test
+*NODE, NSET=NALL
+   1, -0.5, -0.5, -0.5
+   2,  0.5, -0.5, -0.5
+   3,  0.5,  0.5, -0.5
+   4, -0.5,  0.5, -0.5
+   5, -0.5, -0.5,  0.5
+   6,  0.5, -0.5,  0.5
+   7,  0.5,  0.5,  0.5
+   8, -0.5,  0.5,  0.5
+   9,  0.0, -0.5, -0.5
+  10,  0.5,  0.0, -0.5
+  11,  0.0,  0.5, -0.5
+  12, -0.5,  0.0, -0.5
+  13,  0.0, -0.5,  0.5
+  14,  0.5,  0.0,  0.5
+  15,  0.0,  0.5,  0.5
+  16, -0.5,  0.0,  0.5
+  17, -0.5, -0.5,  0.0
+  18,  0.5, -0.5,  0.0
+  19,  0.5,  0.5,  0.0
+  20, -0.5,  0.5,  0.0
+*ELEMENT, TYPE=C3D20, ELSET=EALL
+1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
+11, 12, 13, 14, 15, 16, 17, 18, 19, 20
+**
+** Use this for no material model rotation
+** *SOLID SECTION, ELSET=EALL, MATERIAL=UMAT
+** ,
+**
+** Use this section for a 30 degree material model rotation
+** *ORIENTATION, NAME=30_DEGREE_Z, SYSTEM=RECTANGULAR
+** 0.8660254, 0.5, 0.,   -0.5, 0.8660254, 0.,   0., 0., 0.
+** 3, 0.0
+** *SOLID SECTION, ELSET=EALL, MATERIAL=UMAT, ORIENTATION=30_DEGREE_Z
+** ,
+**
+** Use this section for a 0 degree material model rotation but stress calculated in rotated configuration
+*ORIENTATION, NAME=0_DEGREE_Z, SYSTEM=RECTANGULAR
+1., 0., 0.,   0., 1., 0.,   0., 0., 0.
+3, 0.0
+*SOLID SECTION, ELSET=EALL, MATERIAL=UMAT, ORIENTATION=0_DEGREE_Z
+,
+**
+** Node Set for bottom face
+*NSET, NSET=BOTTOM
+   1, 2, 5, 6, 9, 13, 17, 18
+** Node Set for top face
+*NSET, NSET=TOP
+   3, 4, 7, 8, 11, 15, 19, 20
+** Hold bottom fixed
+*BOUNDARY
+BOTTOM, 1, 3, 0.0
+** Call UMAT
+*MATERIAL, NAME=UMAT
+*USER MATERIAL, CONSTANTS=2
+1000.0, 0.3
+** Shear top face in +x direction by 0.005 over 5 seconds
+*STEP, NAME=STEP-1, NLGEOM=YES
+*STATIC
+1.0, 5.0, 1.0, 1.0
+*BOUNDARY
+TOP, 1, 1, 0.005
+*OUTPUT, FIELD
+*NODE OUTPUT, NSET=NALL
+U
+COORD
+*ELEMENT OUTPUT, ELSET=EALL, POSITION=INTEGRATION POINTS
+S
+COORD
+*END STEP

modules/solid_mechanics/test/tests/umat/orient_umat/shear_top_umat.i

@@ -0,0 +1,103 @@
+[GlobalParams]
+  displacements = 'disp_x disp_y disp_z'
+  decomposition_method = HughesWinget
+[]
+
+[Mesh]
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 3
+    xmin = -0.5
+    xmax = 0.5
+    ymin = -0.5
+    ymax = 0.5
+    zmin = -0.5
+    zmax = 0.5
+    elem_type = HEX20
+  []
+[]
+
+[Functions]
+  [top_pull]
+    type = ParsedFunction
+    expression = t/1000
+  []
+[]
+
+[Physics/SolidMechanics/QuasiStatic]
+  [all]
+    add_variables = true
+    strain = FINITE
+    generate_output = 'stress_xx stress_yy stress_xy'
+  []
+[]
+
+[BCs]
+  [x_pull]
+    type = FunctionDirichletBC
+    variable = disp_x
+    boundary = top
+    function = top_pull
+  []
+  [x_bot]
+    type = DirichletBC
+    variable = disp_x
+    boundary = bottom
+    value = 0.0
+  []
+  [y_bot]
+    type = DirichletBC
+    variable = disp_y
+    boundary = bottom
+    value = 0.0
+  []
+  [z_bot]
+    type = DirichletBC
+    variable = disp_z
+    boundary = bottom
+    value = 0.0
+  []
+[]
+
+[Materials]
+  [umat]
+    type = AbaqusUMATStress
+    constant_properties = '1000 0.3'
+    plugin = '../../../plugins/elastic_print'
+    num_state_vars = 0
+    use_one_based_indexing = true
+    use_displaced_mesh = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'PJFNK'
+
+  petsc_options = '-snes_ksp_ew'
+  petsc_options_iname = '-ksp_gmres_restart'
+  petsc_options_value = '101'
+
+  line_search = 'none'
+
+  l_max_its = 100
+  nl_max_its = 100
+  nl_rel_tol = 1e-12
+  nl_abs_tol = 1e-10
+  l_tol = 1e-9
+  start_time = 0.0
+  end_time = 5
+
+  dt = 1.0
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+
+[Outputs]
+  exodus = true
+[]