diff --git a/modules/solid_mechanics/doc/content/source/materials/EigenDecompositionMaterial.md b/modules/solid_mechanics/doc/content/source/materials/EigenDecompositionMaterial.md new file mode 100644 index 000000000000..abc0b210ab69 --- /dev/null +++ b/modules/solid_mechanics/doc/content/source/materials/EigenDecompositionMaterial.md @@ -0,0 +1,17 @@ +# EigenDecompositionMaterial + +!syntax description /Materials/EigenDecompositionMaterial + +## Description + +This class reads in a symmetric `RankTwoTensor` material given by [!param](/Materials/EigenDecompositionMaterial/rank_two_tensor) and performs an eigendecomposition on it. The results of the decomposition are scalar materials named `max_eigen_value`, `mid_eigen_value`, and `min_eigen_value` and the corresponding vector material properties named `max_eigen_vector`, `mid_eigen_vector`, and `min_eigen_vector`. These names can be preceeded by the [!param](/Materials/EigenDecompositionMaterial/base_name) to allow for more than one `EigenDecompositionMaterial` per block. These material properties can be output using the Material Outputs system with the [!param](/Materials/EigenDecompositionMaterial/output_properties) as shown in the below example. This material is useful for visualizing the elemental maximum principal stress and direction as a vector. An error will be produced if an unsymmetric `RankTwoTensor` is decomposed which can happen for the deformation gradient. + +## Example Input File Syntax + +!listing modules/solid_mechanics/test/tests/eigen_decomp_material/prescribed_strain_3D.i block=Materials/nonADeig_decomp + +!syntax parameters /Materials/EigenDecompositionMaterial + +!syntax inputs /Materials/EigenDecompositionMaterial + +!syntax children /Materials/EigenDecompositionMaterial diff --git a/modules/solid_mechanics/include/materials/EigenDecompositionMaterial.h b/modules/solid_mechanics/include/materials/EigenDecompositionMaterial.h new file mode 100644 index 000000000000..abe33fc2948e --- /dev/null +++ b/modules/solid_mechanics/include/materials/EigenDecompositionMaterial.h @@ -0,0 +1,44 @@ +//* 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 "Material.h" +#include "RankTwoTensorForward.h" +/** + * Perform eigendecomposition on a RankTwoTensor material property + */ +template +class EigenDecompositionMaterialTempl : public Material +{ +public: + static InputParameters validParams(); + + EigenDecompositionMaterialTempl(const InputParameters & parameters); + +protected: + virtual void computeQpProperties() override; + /// Base name to allow multiple tensors to be decomposed + const std::string _base_name; + /// Rank two tensor for eigen decomposition + const GenericMaterialProperty & _tensor; + + /// eigen vectors + GenericMaterialProperty & _max_eigen_vector; + GenericMaterialProperty & _mid_eigen_vector; + GenericMaterialProperty & _min_eigen_vector; + + /// eigen values + GenericMaterialProperty & _max_eigen_value; + GenericMaterialProperty & _mid_eigen_value; + GenericMaterialProperty & _min_eigen_value; +}; + +typedef EigenDecompositionMaterialTempl EigenDecompositionMaterial; +typedef EigenDecompositionMaterialTempl ADEigenDecompositionMaterial; diff --git a/modules/solid_mechanics/src/materials/EigenDecompositionMaterial.C b/modules/solid_mechanics/src/materials/EigenDecompositionMaterial.C new file mode 100644 index 000000000000..ded508bbfec3 --- /dev/null +++ b/modules/solid_mechanics/src/materials/EigenDecompositionMaterial.C @@ -0,0 +1,72 @@ +//* 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 "EigenDecompositionMaterial.h" +#include "RankTwoTensor.h" + +registerMooseObject("SolidMechanicsApp", EigenDecompositionMaterial); +registerMooseObject("SolidMechanicsApp", ADEigenDecompositionMaterial); + +template +InputParameters +EigenDecompositionMaterialTempl::validParams() +{ + InputParameters params = Material::validParams(); + params.addClassDescription("Emits material properties for the eigenvalues and eigenvectors of a " + "symmetric rank two tensor."); + params.addRequiredParam( + "rank_two_tensor", + "The name of the symmetric rank two tensor to used in eigen decomposition."); + params.addParam( + "base_name", + "Optional parameter to allow multiple tensors to be decomposed on the same block."); + return params; +} + +template +EigenDecompositionMaterialTempl::EigenDecompositionMaterialTempl( + const InputParameters & parameters) + : Material(parameters), + _base_name(isParamValid("base_name") ? getParam("base_name") + "_" : ""), + _tensor(getGenericMaterialProperty("rank_two_tensor")), + _max_eigen_vector( + declareGenericProperty(_base_name + "max_eigen_vector")), + _mid_eigen_vector( + declareGenericProperty(_base_name + "mid_eigen_vector")), + _min_eigen_vector( + declareGenericProperty(_base_name + "min_eigen_vector")), + _max_eigen_value(declareGenericProperty(_base_name + "max_eigen_value")), + _mid_eigen_value(declareGenericProperty(_base_name + "mid_eigen_value")), + _min_eigen_value(declareGenericProperty(_base_name + "min_eigen_value")) +{ + if (LIBMESH_DIM != 3) + mooseError("EigenDecompositionMaterial is only defined for LIBMESH_DIM=3"); +} + +template +void +EigenDecompositionMaterialTempl::computeQpProperties() +{ + + if (!_tensor[_qp].isSymmetric()) + mooseError("EigenDecompositionMaterial will only operate on symmetric rank two tensors."); + + std::vector> eigval(3, 0.0); + GenericRankTwoTensor eigvec; + + _tensor[_qp].symmetricEigenvaluesEigenvectors(eigval, eigvec); + + _max_eigen_vector[_qp] = eigvec.column(2); + _mid_eigen_vector[_qp] = eigvec.column(1); + _min_eigen_vector[_qp] = eigvec.column(0); + + _max_eigen_value[_qp] = eigval[2]; + _mid_eigen_value[_qp] = eigval[1]; + _min_eigen_value[_qp] = eigval[0]; +} diff --git a/modules/solid_mechanics/test/tests/eigen_decomp_material/gold/pipe_2d_out.e b/modules/solid_mechanics/test/tests/eigen_decomp_material/gold/pipe_2d_out.e new file mode 100644 index 000000000000..59ddffe1781e Binary files /dev/null and b/modules/solid_mechanics/test/tests/eigen_decomp_material/gold/pipe_2d_out.e differ diff --git a/modules/solid_mechanics/test/tests/eigen_decomp_material/gold/prescribed_strain_3D_out.e b/modules/solid_mechanics/test/tests/eigen_decomp_material/gold/prescribed_strain_3D_out.e new file mode 100644 index 000000000000..b646464ffc5b Binary files /dev/null and b/modules/solid_mechanics/test/tests/eigen_decomp_material/gold/prescribed_strain_3D_out.e differ diff --git a/modules/solid_mechanics/test/tests/eigen_decomp_material/pipe_2d.i b/modules/solid_mechanics/test/tests/eigen_decomp_material/pipe_2d.i new file mode 100644 index 000000000000..49af09f8fd22 --- /dev/null +++ b/modules/solid_mechanics/test/tests/eigen_decomp_material/pipe_2d.i @@ -0,0 +1,125 @@ + +# pipe with outer diameter = 24", wall thickness 0.979" +# units in MPa, meters +[GlobalParams] + order = SECOND + family = LAGRANGE + displacements = 'disp_x disp_y' +[] + +[Mesh] + coord_type = RZ + [pipe] + type = GeneratedMeshGenerator + dim = 2 + xmin = 0.2799334 #12"-0.979" + xmax = .3048 #12" + ymin = 0 + ymax = 0.00497 #0.979"/5 + nx = 5 + ny = 1 + elem_type = quad9 + [] +[] + +[Variables] + [disp_x] + [] + [disp_y] + [] +[] + +[Physics] + [SolidMechanics] + [QuasiStatic] + [all] + strain = SMALL + add_variables = true + new_system = true + formulation = TOTAL + volumetric_locking_correction = false + [] + [] + [] +[] + +[Functions] + [inner_pressure] + type = ConstantFunction + value = 3 + [] + [big_pipe] + type = ParsedFunction + expression = -(p*11.021*11.021)/(12*12-11.021*11.021) + symbol_names = 'p' + symbol_values = 3 + [] +[] + +[BCs] + [fixBottom] + type = DirichletBC + variable = disp_y + boundary = bottom + value = 0.0 + preset = true + [] + [Pressure] + [inside] + boundary = left + function = inner_pressure + [] + [axial] + boundary = top + function = big_pipe + [] + [] +[] + +[Constraints] + [top] + type = EqualValueBoundaryConstraint + variable = disp_y + secondary = top + penalty = 1e+14 + formulation = penalty + [] +[] + +[Materials] + [elastic_tensor] + type = ComputeIsotropicElasticityTensor + youngs_modulus = 168000 + poissons_ratio = 0.31 + [] + [compute_stress] + type = ComputeLagrangianLinearElasticStress + [] + [nonADeig_decomp] + type = EigenDecompositionMaterial + rank_two_tensor = cauchy_stress + outputs = exodus + output_properties = "max_eigen_vector mid_eigen_vector min_eigen_vector " + "max_eigen_valuemid_eigen_value min_eigen_value" + [] +[] + +[Executioner] + type = Steady + solve_type = 'NEWTON' + petsc_options_iname = '-pc_type -pc_factor_mat_solver_package' + petsc_options_value = 'lu superlu_dist' + nl_rel_tol = 5e-8 +[] + +[Postprocessors] + [eigval_max] + type = ElementAverageMaterialProperty + mat_prop = max_eigen_value + [] +[] + +[Outputs] + exodus = on + console = true +[] diff --git a/modules/solid_mechanics/test/tests/eigen_decomp_material/prescribed_strain_3D.i b/modules/solid_mechanics/test/tests/eigen_decomp_material/prescribed_strain_3D.i new file mode 100644 index 000000000000..cc3d406baa51 --- /dev/null +++ b/modules/solid_mechanics/test/tests/eigen_decomp_material/prescribed_strain_3D.i @@ -0,0 +1,156 @@ +[Debug] + show_material_props = true +[] + +[Problem] + kernel_coverage_check = false + solve = false +[] + +[Mesh] + type = GeneratedMesh + dim = 3 + nx = 1 + ny = 1 + nz = 1 + displacements = 'disp_x disp_y disp_z' +[] + +[AuxVariables] + [disp_y] + initial_condition = 0 + [] + [disp_x] + initial_condition = 0 + [] + [disp_z] + initial_condition = 0 + [] +[] + +[AuxKernels] + # The applied displacements will cause the max eigenvector to change directions. + # At t=5, the body undergoes simple shear, producing a nonsymmetric deformation gradient. + [disp_x] + execute_on = 'TIMESTEP_BEGIN' + type = ParsedAux + variable = disp_x + use_xyzt = true + expression = "if(t<4.1,4e-1*x*t,x)" + [] + [disp_y] + execute_on = 'TIMESTEP_BEGIN' + type = ParsedAux + variable = disp_y + use_xyzt = true + expression = "if(t<4.1,3e-1*y*t^2,1e-1*y*t+1e-1*x*t)" + [] + [disp_z] + execute_on = 'TIMESTEP_BEGIN' + type = ParsedAux + variable = disp_z + use_xyzt = true + expression = "if(t<4.1,1e-1*z*t^3,z)" + [] +[] + +[Materials] + [compute_strain] + type = ComputeLagrangianStrain + displacements = 'disp_x disp_y disp_z' + large_kinematics = true + [] + [nonAD_strain] + type = RankTwoTensorMaterialADConverter + reg_props_in = mechanical_strain + ad_props_out = AD_mechanical_strain + [] + [eig_decomp] + type = ADEigenDecompositionMaterial + rank_two_tensor = AD_mechanical_strain + outputs = exodus + output_properties = "max_eigen_vector mid_eigen_vector min_eigen_vector max_eigen_value " + "mid_eigen_value min_eigen_value" + [] + [nonADeig_decomp] + type = EigenDecompositionMaterial + rank_two_tensor = mechanical_strain + base_name = nonAD + outputs = exodus + output_properties = "nonAD_max_eigen_vector nonAD_mid_eigen_vector nonAD_min_eigen_vector " + "nonAD_max_eigen_value nonAD_mid_eigen_value nonAD_min_eigen_value" + [] + + [non_symmetric_eig_decomp_error] + type = EigenDecompositionMaterial + rank_two_tensor = deformation_gradient + base_name = nonSym + [] +[] + +[BCs] +[] + +[Executioner] + type = Transient + solve_type = LINEAR + dt = 1 + end_time = 5 +[] + +[Postprocessors] + [sxx] + type = MaterialTensorAverage + rank_two_tensor = mechanical_strain + index_i = 0 + index_j = 0 + execute_on = 'TIMESTEP_END' + [] + [sxy] + type = MaterialTensorAverage + rank_two_tensor = mechanical_strain + index_i = 0 + index_j = 1 + execute_on = 'TIMESTEP_END' + [] + [syy] + type = MaterialTensorAverage + rank_two_tensor = mechanical_strain + index_i = 1 + index_j = 1 + [] + [szz] + type = MaterialTensorAverage + rank_two_tensor = mechanical_strain + index_i = 1 + index_j = 0 + [] + [AD_eigval_max] + type = ADElementAverageMaterialProperty + mat_prop = max_eigen_value + [] + [AD_eigval_mid] + type = ADElementAverageMaterialProperty + mat_prop = mid_eigen_value + [] + [AD_eigval_min] + type = ADElementAverageMaterialProperty + mat_prop = min_eigen_value + [] + [nonAD_eigval_max] + type = ElementAverageMaterialProperty + mat_prop = nonAD_max_eigen_value + [] + [nonAD_eigval_mid] + type = ElementAverageMaterialProperty + mat_prop = nonAD_mid_eigen_value + [] + [nonAD_eigval_min] + type = ElementAverageMaterialProperty + mat_prop = nonAD_min_eigen_value + [] +[] + +[Outputs] + exodus = true +[] diff --git a/modules/solid_mechanics/test/tests/eigen_decomp_material/tests b/modules/solid_mechanics/test/tests/eigen_decomp_material/tests new file mode 100644 index 000000000000..57c61578bfa7 --- /dev/null +++ b/modules/solid_mechanics/test/tests/eigen_decomp_material/tests @@ -0,0 +1,23 @@ +[Tests] + design = 'source/materials/EigenDecompositionMaterial.md' + issues = '#29179' + [3D] + type = Exodiff + input = prescribed_strain_3D.i + exodiff = prescribed_strain_3D_out.e + cli_args = 'Materials/inactive=non_symmetric_eig_decomp_error' + requirement = 'The system shall output an eigen decomposition of a symmetric rank two tensor in a 3D simulation with an evolving strain field.' + [] + [nonSymmetricError] + type = RunException + input = 'prescribed_strain_3D.i' + expect_err = 'EigenDecompositionMaterial will only operate on symmetric rank two tensors.' + requirement = 'The system shall produce an error if the rank two tensor is unsymmetric.' + [] + [2Drz] + type = Exodiff + input = pipe_2d.i + exodiff =pipe_2d_out.e + requirement = 'The system shall output an eigen decomposition of a symmetric rank two tensor in a 2D radial symmetry simulation of a pressurized pipe.' + [] +[]