diff --git a/modules/solid_mechanics/doc/content/modules/solid_mechanics/ShellElements.md b/modules/solid_mechanics/doc/content/modules/solid_mechanics/ShellElements.md index ca9e740d78be..0c28fffe2327 100644 --- a/modules/solid_mechanics/doc/content/modules/solid_mechanics/ShellElements.md +++ b/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 diff --git a/modules/solid_mechanics/doc/content/source/auxkernels/ShellLocalCoordinatesAux.md b/modules/solid_mechanics/doc/content/source/auxkernels/ShellLocalCoordinatesAux.md new file mode 100644 index 000000000000..572b84f48fb4 --- /dev/null +++ b/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 diff --git a/modules/solid_mechanics/doc/content/source/auxkernels/ShellResultantsAux.md b/modules/solid_mechanics/doc/content/source/auxkernels/ShellResultantsAux.md new file mode 100644 index 000000000000..032df332d9c1 --- /dev/null +++ b/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 diff --git a/modules/solid_mechanics/include/auxkernels/ShellLocalCoordinatesAux.h b/modules/solid_mechanics/include/auxkernels/ShellLocalCoordinatesAux.h new file mode 100644 index 000000000000..7363bc0559f4 --- /dev/null +++ b/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 * _local_coordinates; +}; diff --git a/modules/solid_mechanics/include/auxkernels/ShellResultantsAux.h b/modules/solid_mechanics/include/auxkernels/ShellResultantsAux.h new file mode 100644 index 000000000000..26f4c1c32505 --- /dev/null +++ b/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 _t_qrule; + + /// Quadrature points and weights in the out of plane direction in isoparametric coordinate system + std::vector _t_points; + std::vector _t_weights; + + /// The local stress tensor + std::vector *> _local_stress_t_points; +}; diff --git a/modules/solid_mechanics/include/materials/ADComputeIncrementalShellStrain.h b/modules/solid_mechanics/include/materials/ADComputeIncrementalShellStrain.h index 4cd4df117bc6..6f707b0776cb 100644 --- a/modules/solid_mechanics/include/materials/ADComputeIncrementalShellStrain.h +++ b/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 *> _J_map_old; + /// Transformation matrix to map stresses from global coordinate to the element's local coordinate + std::vector *> _local_transformation_matrix; + std::vector *> _local_transformation_matrix_old; + /// Covariant base vector matrix material property to transform stress std::vector *> _covariant_transformation_matrix; std::vector *> _covariant_transformation_matrix_old; @@ -177,6 +184,7 @@ class ADComputeIncrementalShellStrain : public Material ADRealVectorValue _x3; const NumericVector * const & _sol; const NumericVector & _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; }; diff --git a/modules/solid_mechanics/include/materials/ADComputeShellStress.h b/modules/solid_mechanics/include/materials/ADComputeShellStress.h index b6a80639aac5..6afb590d09b4 100644 --- a/modules/solid_mechanics/include/materials/ADComputeShellStress.h +++ b/modules/solid_mechanics/include/materials/ADComputeShellStress.h @@ -47,12 +47,18 @@ class ADComputeShellStress : public Material /// Quadrature points along thickness std::vector _t_points; + /// Transformation matrix to map the global stress to the element's local coordinate + std::vector *> _local_transformation_matrix; + /// Covariant base vector matrix material property to transform stress std::vector *> _covariant_transformation_matrix; /// Global stress tensor material property std::vector *> _global_stress; + /// local stress tensor material property + std::vector *> _local_stress; + /// Real value of stress in the local coordinate system RankTwoTensor _unrotated_stress; }; diff --git a/modules/solid_mechanics/src/auxkernels/ShellLocalCoordinatesAux.C b/modules/solid_mechanics/src/auxkernels/ShellLocalCoordinatesAux.C new file mode 100644 index 000000000000..c1bdacb5e32b --- /dev/null +++ b/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("base_name", "Mechanical property base name"); + + MooseEnum property("first_local_vector second_local_vector normal_local_vector"); + params.addRequiredParam( + "property", + property, + "The local axis to output: first_local_vector, second_local_vector or normal_local_vector"); + params.addRequiredParam( + "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("base_name") + "_" : ""), + _property(getParam("property").getEnum()), + _component(getParam("component")) + +{ + _local_coordinates = + &getMaterialProperty(_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; +} diff --git a/modules/solid_mechanics/src/auxkernels/ShellResultantsAux.C b/modules/solid_mechanics/src/auxkernels/ShellResultantsAux.C new file mode 100644 index 000000000000..00e0f42f4e17 --- /dev/null +++ b/modules/solid_mechanics/src/auxkernels/ShellResultantsAux.C @@ -0,0 +1,142 @@ +//* 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 "ShellResultantsAux.h" +#include "ArbitraryQuadrature.h" +#include "libmesh/quadrature.h" +#include "libmesh/utility.h" +#include "libmesh/enum_quadrature_type.h" +#include "libmesh/fe_type.h" +#include "libmesh/string_to_enum.h" +#include "libmesh/quadrature_gauss.h" + +registerMooseObject("SolidMechanicsApp", ShellResultantsAux); + +InputParameters +ShellResultantsAux::validParams() +{ + InputParameters params = AuxKernel::validParams(); + params.addClassDescription( + "Computes the local forces, bending moments and shear forces acting on shell elements"); + params.addParam("base_name", "Mechanical property base name"); + params.addRequiredCoupledVar( + "thickness", + "Thickness of the shell. Can be supplied as either a number or a variable name."); + params.addRequiredParam("through_thickness_order", + "Quadrature order in out of plane direction"); + MooseEnum stress_resultant( + "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"); + params.addRequiredParam( + "stress_resultant", + stress_resultant, + "The stress resultant type to output, calculated on the shell element."); + + return params; +} + +ShellResultantsAux::ShellResultantsAux(const InputParameters & parameters) + : AuxKernel(parameters), + _base_name(isParamValid("base_name") ? getParam("base_name") + "_" : ""), + _thickness(coupledValue("thickness")), + _resultant(getParam("stress_resultant").getEnum()) +{ + _t_qrule = std::make_unique( + 1, Utility::string_to_enum(getParam("through_thickness_order"))); + _t_points = _t_qrule->get_points(); + _t_weights = _t_qrule->get_weights(); + _local_stress_t_points.resize(_t_points.size()); + + for (const auto t : index_range(_t_points)) + _local_stress_t_points[t] = &getMaterialProperty( + _base_name + "local_stress_t_points_" + std::to_string(t)); +} + +Real +ShellResultantsAux::computeValue() +{ + Real _shell_resultant = 0.0; + + switch (_resultant) + { + case ResultantType::axial_force_0: + for (unsigned int i = 0; i < _t_points.size(); ++i) + { + _shell_resultant += + (*_local_stress_t_points[i])[_qp](0, 0) * _t_weights[i] * (_thickness[_qp] / 2); + } + break; + + case ResultantType::axial_force_1: + for (unsigned int i = 0; i < _t_points.size(); ++i) + { + _shell_resultant += + (*_local_stress_t_points[i])[_qp](1, 1) * _t_weights[i] * (_thickness[_qp] / 2); + } + break; + + case ResultantType::normal_force: + for (unsigned int i = 0; i < _t_points.size(); ++i) + { + _shell_resultant += + (*_local_stress_t_points[i])[_qp](2, 2) * _t_weights[i] * (_thickness[_qp] / 2); + } + break; + + case ResultantType::bending_moment_0: + for (unsigned int i = 0; i < _t_points.size(); ++i) + { + _shell_resultant -= (*_local_stress_t_points[i])[_qp](1, 1) * _t_points[i](0) * + _t_weights[i] * (_thickness[_qp] / 2) * (_thickness[_qp] / 2); + } + break; + + case ResultantType::bending_moment_1: + for (unsigned int i = 0; i < _t_points.size(); ++i) + { + _shell_resultant -= (*_local_stress_t_points[i])[_qp](0, 0) * _t_points[i](0) * + _t_weights[i] * (_thickness[_qp] / 2) * (_thickness[_qp] / 2); + } + break; + + case ResultantType::bending_moment_01: + for (unsigned int i = 0; i < _t_points.size(); ++i) + { + _shell_resultant -= (*_local_stress_t_points[i])[_qp](0, 1) * _t_points[i](0) * + _t_weights[i] * (_thickness[_qp] / 2) * (_thickness[_qp] / 2); + } + break; + + case ResultantType::shear_force_01: + for (unsigned int i = 0; i < _t_points.size(); ++i) + { + _shell_resultant += + (*_local_stress_t_points[i])[_qp](0, 1) * _t_weights[i] * (_thickness[_qp] / 2); + } + break; + + case ResultantType::shear_force_02: + for (unsigned int i = 0; i < _t_points.size(); ++i) + { + _shell_resultant += + (*_local_stress_t_points[i])[_qp](0, 2) * _t_weights[i] * (_thickness[_qp] / 2); + } + break; + + case ResultantType::shear_force_12: + for (unsigned int i = 0; i < _t_points.size(); ++i) + { + _shell_resultant += + (*_local_stress_t_points[i])[_qp](1, 2) * _t_weights[i] * (_thickness[_qp] / 2); + } + break; + } + + return _shell_resultant; +} diff --git a/modules/solid_mechanics/src/materials/ADComputeIncrementalShellStrain.C b/modules/solid_mechanics/src/materials/ADComputeIncrementalShellStrain.C index 4996ad2245c3..19e9bed7ac17 100644 --- a/modules/solid_mechanics/src/materials/ADComputeIncrementalShellStrain.C +++ b/modules/solid_mechanics/src/materials/ADComputeIncrementalShellStrain.C @@ -41,6 +41,10 @@ ADComputeIncrementalShellStrain::validParams() "Quadrature order in out of plane direction"); params.addParam( "large_strain", false, "Set to true to turn on finite strain calculations."); + params.addParam("reference_first_local_direction", + RealVectorValue(1, 0, 0), + "Vector that is projected onto an element to define the " + "direction for the first local coordinate direction e1."); return params; } @@ -77,13 +81,17 @@ ADComputeIncrementalShellStrain::ADComputeIncrementalShellStrain(const InputPara _ge_old(), _J_map(), _J_map_old(), + _local_transformation_matrix(), + _local_transformation_matrix_old(), _covariant_transformation_matrix(), _covariant_transformation_matrix_old(), _contravariant_transformation_matrix(), _contravariant_transformation_matrix_old(), _total_global_strain(), _sol(_nonlinear_sys.currentSolution()), - _sol_old(_nonlinear_sys.solutionOld()) + _sol_old(_nonlinear_sys.solutionOld()), + _ref_first_local_dir(getParam("reference_first_local_direction")) + { // Checking for consistency between length of the provided displacements and rotations vector if (_ndisp != 3 || _nrot != 2) @@ -95,6 +103,11 @@ ADComputeIncrementalShellStrain::ADComputeIncrementalShellStrain(const InputPara mooseError( "ADComputeIncrementalShellStrain: Shell element is implemented only for linear elements."); + // Checking if the size of the first local vector is within an acceptable range + if (MooseUtils::absoluteFuzzyEqual(_ref_first_local_dir.norm(), 0.0, 1e-3)) + mooseError( + "The norm of the defined first local axis is less than the acceptable telerance (1e-3)"); + // fetch coupled variables and gradients (as stateful properties if necessary) for (unsigned int i = 0; i < _ndisp; ++i) { @@ -126,6 +139,8 @@ ADComputeIncrementalShellStrain::ADComputeIncrementalShellStrain(const InputPara _dxyz_dxi_old.resize(_t_points.size()); _dxyz_deta_old.resize(_t_points.size()); _dxyz_dzeta_old.resize(_t_points.size()); + _local_transformation_matrix.resize(_t_points.size()); + _local_transformation_matrix_old.resize(_t_points.size()); _covariant_transformation_matrix.resize(_t_points.size()); _covariant_transformation_matrix_old.resize(_t_points.size()); _contravariant_transformation_matrix.resize(_t_points.size()); @@ -159,6 +174,10 @@ ADComputeIncrementalShellStrain::ADComputeIncrementalShellStrain(const InputPara _dxyz_dzeta_old[i] = &getMaterialPropertyOldByName("dxyz_dzeta_t_points_" + std::to_string(i)); // Create rotation matrix and total strain global for output purposes only + _local_transformation_matrix[i] = + &declareProperty("local_transformation_t_points_" + std::to_string(i)); + _local_transformation_matrix_old[i] = &getMaterialPropertyOldByName( + "local_transformation_t_points_" + std::to_string(i)); _covariant_transformation_matrix[i] = &declareProperty("covariant_transformation_t_points_" + std::to_string(i)); _covariant_transformation_matrix_old[i] = &getMaterialPropertyOldByName( @@ -174,6 +193,9 @@ ADComputeIncrementalShellStrain::ADComputeIncrementalShellStrain(const InputPara // used later for computing local coordinate system _x2(1) = 1; _x3(2) = 1; + _e1(0) = 1; + _e2(1) = 1; + _e3(2) = 1; } void @@ -188,6 +210,7 @@ ADComputeIncrementalShellStrain::initQpStatefulProperties() if (_qrule->get_points().size() != 4) mooseError("ADComputeIncrementalShellStrain: Shell element needs to have exactly four " "quadrature points."); + computeGMatrix(); computeBMatrix(); } @@ -221,6 +244,7 @@ ADComputeIncrementalShellStrain::computeProperties() (*_dxyz_deta[j])[i] = (*_dxyz_deta_old[j])[i]; (*_dxyz_dzeta[j])[i] = (*_dxyz_dzeta_old[j])[i]; (*_B[j])[i] = (*_B_old[j])[i]; + (*_local_transformation_matrix[j])[i] = (*_local_transformation_matrix_old[j])[i]; (*_covariant_transformation_matrix[j])[i] = (*_covariant_transformation_matrix_old[j])[i]; (*_contravariant_transformation_matrix[j])[i] = (*_contravariant_transformation_matrix_old[j])[i]; @@ -280,7 +304,7 @@ ADComputeIncrementalShellStrain::computeGMatrix() _dphideta_map = fe->get_fe_map().get_dphideta_map(); _phi_map = fe->get_fe_map().get_phi_map(); - for (unsigned int i = 0; i < 4; ++i) + for (unsigned int i = 0; i < _nodes.size(); ++i) _nodes[i] = _current_elem->node_ptr(i); ADRealVectorValue x = (*_nodes[1] - *_nodes[0]); @@ -295,6 +319,7 @@ ADComputeIncrementalShellStrain::computeGMatrix() ADDenseMatrix b(5, 20); ADRealVectorValue c; RankTwoTensor d; + RealVectorValue f; for (unsigned int t = 0; t < _t_points.size(); ++t) { (*_strain_increment[t])[_qp] = a; @@ -305,6 +330,7 @@ ADComputeIncrementalShellStrain::computeGMatrix() (*_dxyz_dxi[t])[_qp] = c; (*_dxyz_deta[t])[_qp] = c; (*_dxyz_dzeta[t])[_qp] = c; + (*_local_transformation_matrix[t])[_qp] = d; (*_covariant_transformation_matrix[t])[_qp] = d; (*_contravariant_transformation_matrix[t])[_qp] = d; } @@ -388,12 +414,17 @@ ADComputeIncrementalShellStrain::computeGMatrix() // _transformation_matrix is an AD property, but we're not setting the derivatives! (*_transformation_matrix)[i] = J; - // calculate ge - ADRealVectorValue e3 = (*_dxyz_dzeta[j])[i] / (*_dxyz_dzeta[j])[i].norm(); - ADRealVectorValue e1 = (*_dxyz_deta[j])[i].cross(e3); - e1 /= e1.norm(); - ADRealVectorValue e2 = e3.cross(e1); - e2 /= e2.norm(); + // compute element's local coordinate + computeLocalCoordinates(); + + // calculate the local transformation matrix to be used to map the global stresses to the + // local element coordinate + const auto local_rotation_mat = ADRankTwoTensor::initializeFromRows(_e1, _e2, _e3); + + for (const auto ii : make_range(Moose::dim)) + for (const auto jj : make_range(Moose::dim)) + (*_local_transformation_matrix[j])[i](ii, jj) = + MetaPhysicL::raw_value(local_rotation_mat(ii, jj)); // Calculate the contravariant base vectors g^0, g^1, g^2 // The base vectors for the strain tensor in the convected coordinate @@ -410,19 +441,51 @@ ADComputeIncrementalShellStrain::computeGMatrix() // the strain tensor (with g_1, g2_, g_3 base vectors) to // the stress tensor (with base vectors e1, e2, e3) - (*_ge[j])[i](0, 0) = gi0 * e1; - (*_ge[j])[i](0, 1) = gi0 * e2; - (*_ge[j])[i](0, 2) = gi0 * e3; - (*_ge[j])[i](1, 0) = gi1 * e1; - (*_ge[j])[i](1, 1) = gi1 * e2; - (*_ge[j])[i](1, 2) = gi1 * e3; - (*_ge[j])[i](2, 0) = gi2 * e1; - (*_ge[j])[i](2, 1) = gi2 * e2; - (*_ge[j])[i](2, 2) = gi2 * e3; + (*_ge[j])[i](0, 0) = gi0 * _e1; + (*_ge[j])[i](0, 1) = gi0 * _e2; + (*_ge[j])[i](0, 2) = gi0 * _e3; + (*_ge[j])[i](1, 0) = gi1 * _e1; + (*_ge[j])[i](1, 1) = gi1 * _e2; + (*_ge[j])[i](1, 2) = gi1 * _e3; + (*_ge[j])[i](2, 0) = gi2 * _e1; + (*_ge[j])[i](2, 1) = gi2 * _e2; + (*_ge[j])[i](2, 2) = gi2 * _e3; } } } +void +ADComputeIncrementalShellStrain::computeLocalCoordinates() +{ + // default option for the first local vector:the global X-axis is projected on the shell plane + // alternatively the reference first local vector provided by the user can be used to define the + // 1st local axis + + // All nodes in an element have the same normal vector. Therefore, here, we use only the normal + // vecor of the first node as "the element's normal vector" + _e3 = _node_normal[0]; + + _e1 = _ref_first_local_dir; + + _e1 /= _e1.norm(); + + // The reference first local axis and the normal are considered parallel if the angle between them + // is less than 0.05 degrees + if (MooseUtils::absoluteFuzzyEqual(std::abs(_e1 * _e3), 1.0, 0.05 * libMesh::pi / 180.0)) + { + mooseError("The reference first local axis must not be perpendicular to any of the shell " + "elements "); + } + + // we project the reference first local vector on the shell element and calculate the in-plane e1 + // and e2 vectors + _e2 = _e3.cross(_e1); + _e2 /= _e2.norm(); + + _e1 = _e2.cross(_e3); + _e1 /= _e1.norm(); +} + void ADComputeIncrementalShellStrain::computeNodeNormal() { @@ -434,17 +497,11 @@ void ADComputeIncrementalShellStrain::computeBMatrix() { // compute nodal local axis + computeLocalCoordinates(); for (unsigned int k = 0; k < _nodes.size(); ++k) { - _v1[k] = _x2.cross(_node_normal[k]); - - // If x2 is parallel to node normal, set V1 to x3 - if (MooseUtils::absoluteFuzzyEqual(_v1[k].norm(), 0.0, 1e-6)) - _v1[k] = _x3; - - _v1[k] /= _v1[k].norm(); - - _v2[k] = _node_normal[k].cross(_v1[k]); + _v1[k] = _e1; + _v2[k] = _e2; } // compute B matrix rows correspond to [ux1, ux2, ux3, ux4, uy1, uy2, uy3, uy4, uz1, uz2, uz3, diff --git a/modules/solid_mechanics/src/materials/ADComputeShellStress.C b/modules/solid_mechanics/src/materials/ADComputeShellStress.C index 1c945e383cfd..a9f918aa694a 100644 --- a/modules/solid_mechanics/src/materials/ADComputeShellStress.C +++ b/modules/solid_mechanics/src/materials/ADComputeShellStress.C @@ -42,8 +42,10 @@ ADComputeShellStress::ADComputeShellStress(const InputParameters & parameters) _stress.resize(_t_points.size()); _stress_old.resize(_t_points.size()); _strain_increment.resize(_t_points.size()); + _local_transformation_matrix.resize(_t_points.size()); _covariant_transformation_matrix.resize(_t_points.size()); _global_stress.resize(_t_points.size()); + _local_stress.resize(_t_points.size()); for (unsigned int t = 0; t < _t_points.size(); ++t) { _elasticity_tensor[t] = @@ -54,10 +56,14 @@ ADComputeShellStress::ADComputeShellStress(const InputParameters & parameters) _strain_increment[t] = &getADMaterialProperty("strain_increment_t_points_" + std::to_string(t)); // rotation matrix and stress for output purposes only + _local_transformation_matrix[t] = + &getMaterialProperty("local_transformation_t_points_" + std::to_string(t)); _covariant_transformation_matrix[t] = &getMaterialProperty( "covariant_transformation_t_points_" + std::to_string(t)); _global_stress[t] = &declareProperty("global_stress_t_points_" + std::to_string(t)); + _local_stress[t] = + &declareProperty("local_stress_t_points_" + std::to_string(t)); } } @@ -82,5 +88,7 @@ ADComputeShellStress::computeQpProperties() _unrotated_stress(ii, jj) = MetaPhysicL::raw_value((*_stress[i])[_qp](ii, jj)); (*_global_stress[i])[_qp] = (*_covariant_transformation_matrix[i])[_qp].transpose() * _unrotated_stress * (*_covariant_transformation_matrix[i])[_qp]; + (*_local_stress[i])[_qp] = (*_local_transformation_matrix[i])[_qp] * (*_global_stress[i])[_qp] * + (*_local_transformation_matrix[i])[_qp].transpose(); } } diff --git a/modules/solid_mechanics/test/tests/shell/static/Plate_Concentrated_Loads.msh b/modules/solid_mechanics/test/tests/shell/static/Plate_Concentrated_Loads.msh new file mode 100644 index 000000000000..ab796af573a4 --- /dev/null +++ b/modules/solid_mechanics/test/tests/shell/static/Plate_Concentrated_Loads.msh @@ -0,0 +1,213 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$PhysicalNames +5 +1 5 "right" +1 6 "left" +1 7 "top" +1 8 "bottom" +2 9 "shell" +$EndPhysicalNames +$Entities +4 4 1 0 +1 0 0 0 0 +2 9 0 0 0 +3 9 0 1 0 +4 0 0 1 0 +1 -1.000000002804313e-07 -1e-07 -1e-07 9.000000099999999 1e-07 1e-07 1 8 2 1 -2 +2 8.999999900000001 -1e-07 -9.999999994736442e-08 9.000000099999999 1e-07 1.0000001 1 5 2 2 -3 +3 -1.000000002804313e-07 -1e-07 0.9999999000000001 9.000000099999999 1e-07 1.0000001 1 7 2 3 -4 +4 -1e-07 -1e-07 -9.999999994736442e-08 1e-07 1e-07 1.0000001 1 6 2 4 -1 +1 -0.4500001000000005 -1e-07 -0.05000010000000088 9.450000100000004 1e-07 1.050000100000001 1 9 4 4 1 2 3 +$EndEntities +$Nodes +9 52 1 52 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +9 0 0 +0 3 0 1 +3 +9 0 1 +0 4 0 1 +4 +0 0 1 +1 1 0 11 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +0.75 0 0 +1.5 0 0 +2.25 0 0 +3 0 0 +3.75 0 0 +4.5 0 0 +5.25 0 0 +6 0 0 +6.75 0 0 +7.5 0 0 +8.25 0 0 +1 2 0 2 +16 +17 +9 0 0.3333333333333333 +9 0 0.6666666666666666 +1 3 0 11 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +8.25 0 1 +7.5 0 1 +6.75 0 1 +6 0 1 +5.25 0 1 +4.5 0 1 +3.75 0 1 +3 0 1 +2.25 0 1 +1.5 0 1 +0.75 0 1 +1 4 0 2 +29 +30 +0 0 0.6666666666666667 +0 0 0.3333333333333334 +2 1 0 22 +31 +32 +33 +34 +35 +36 +37 +38 +39 +40 +41 +42 +43 +44 +45 +46 +47 +48 +49 +50 +51 +52 +0.7499999999999982 0 0.6666666666666664 +1.499999999999999 0 0.6666666666666664 +2.249999999999999 0 0.666666666666667 +2.999999999999998 0 0.6666666666666669 +3.750000000000001 0 0.6666666666666672 +4.5 0 0.6666666666666672 +5.250000000000001 0 0.6666666666666674 +6 0 0.6666666666666674 +6.750000000000004 0 0.6666666666666675 +7.499999999999999 0 0.6666666666666672 +8.250000000000002 0 0.666666666666667 +0.7499999999999984 0 0.3333333333333331 +1.5 0 0.3333333333333333 +2.249999999999998 0 0.3333333333333333 +2.999999999999999 0 0.3333333333333336 +3.749999999999998 0 0.3333333333333339 +4.499999999999999 0 0.3333333333333341 +5.249999999999999 0 0.3333333333333341 +6.000000000000002 0 0.3333333333333343 +6.75 0 0.3333333333333341 +7.500000000000001 0 0.3333333333333341 +8.25 0 0.333333333333334 +$EndNodes +$Elements +5 66 1 66 +1 1 1 12 +1 1 5 +2 5 6 +3 6 7 +4 7 8 +5 8 9 +6 9 10 +7 10 11 +8 11 12 +9 12 13 +10 13 14 +11 14 15 +12 15 2 +1 2 1 3 +13 2 16 +14 16 17 +15 17 3 +1 3 1 12 +16 3 18 +17 18 19 +18 19 20 +19 20 21 +20 21 22 +21 22 23 +22 23 24 +23 24 25 +24 25 26 +25 26 27 +26 27 28 +27 28 4 +1 4 1 3 +28 4 29 +29 29 30 +30 30 1 +2 1 3 36 +31 4 29 31 28 +32 28 31 32 27 +33 27 32 33 26 +34 26 33 34 25 +35 25 34 35 24 +36 24 35 36 23 +37 23 36 37 22 +38 22 37 38 21 +39 21 38 39 20 +40 20 39 40 19 +41 19 40 41 18 +42 18 41 17 3 +43 29 30 42 31 +44 31 42 43 32 +45 32 43 44 33 +46 33 44 45 34 +47 34 45 46 35 +48 35 46 47 36 +49 36 47 48 37 +50 37 48 49 38 +51 38 49 50 39 +52 39 50 51 40 +53 40 51 52 41 +54 41 52 16 17 +55 30 1 5 42 +56 42 5 6 43 +57 43 6 7 44 +58 44 7 8 45 +59 45 8 9 46 +60 46 9 10 47 +61 47 10 11 48 +62 48 11 12 49 +63 49 12 13 50 +64 50 13 14 51 +65 51 14 15 52 +66 52 15 2 16 +$EndElements diff --git a/modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD.i b/modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD.i index dff6df49e398..ceca73273044 100644 --- a/modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD.i +++ b/modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD.i @@ -43,103 +43,103 @@ [] [Variables] - [./disp_x] + [disp_x] order = FIRST family = LAGRANGE - [../] - [./disp_y] + [] + [disp_y] order = FIRST family = LAGRANGE - [../] - [./disp_z] + [] + [disp_z] order = FIRST family = LAGRANGE - [../] - [./rot_x] + [] + [rot_x] order = FIRST family = LAGRANGE - [../] - [./rot_y] + [] + [rot_y] order = FIRST family = LAGRANGE - [../] + [] [] [AuxVariables] - [./stress_yy] + [stress_yy] order = CONSTANT family = MONOMIAL - [../] - [./stress_yz] + [] + [stress_yz] order = CONSTANT family = MONOMIAL - [../] + [] [] [AuxKernels] - [./stress_yy] + [stress_yy] type = RankTwoAux variable = stress_yy rank_two_tensor = global_stress_t_points_0 index_i = 1 index_j = 1 - [../] - [./stress_yz] + [] + [stress_yz] type = RankTwoAux variable = stress_yz rank_two_tensor = global_stress_t_points_0 index_i = 1 index_j = 2 - [../] + [] [] [BCs] - [./fixy1] + [fixy1] type = DirichletBC variable = disp_y boundary = 'bottom' value = 0.0 - [../] - [./fixz1] + [] + [fixz1] type = DirichletBC variable = disp_z boundary = 'bottom' value = 0.0 - [../] - [./fixr1] + [] + [fixr1] type = DirichletBC variable = rot_x boundary = 'bottom' value = 0.0 - [../] - [./fixr2] + [] + [fixr2] type = DirichletBC variable = rot_y boundary = 'bottom' value = 0.0 - [../] - [./fixx1] + [] + [fixx1] type = DirichletBC variable = disp_x boundary = 'bottom' value = 0.0 - [../] + [] [] [NodalKernels] - [./force_y2] + [force_y2] type = ConstantRate variable = disp_z boundary = 'top' rate = 0.5 - [../] + [] [] [Preconditioning] - [./smp] + [smp] type = SMP full = true - [../] + [] [] [Executioner] @@ -155,117 +155,117 @@ [] [Kernels] - [./solid_disp_x] + [solid_disp_x] type = ADStressDivergenceShell block = '0' component = 0 variable = disp_x through_thickness_order = SECOND - [../] - [./solid_disp_y] + [] + [solid_disp_y] type = ADStressDivergenceShell block = '0' component = 1 variable = disp_y through_thickness_order = SECOND - [../] - [./solid_disp_z] + [] + [solid_disp_z] type = ADStressDivergenceShell block = '0' component = 2 variable = disp_z through_thickness_order = SECOND - [../] - [./solid_rot_x] + [] + [solid_rot_x] type = ADStressDivergenceShell block = '0' component = 3 variable = rot_x through_thickness_order = SECOND - [../] - [./solid_rot_y] + [] + [solid_rot_y] type = ADStressDivergenceShell block = '0' component = 4 variable = rot_y through_thickness_order = SECOND - [../] + [] [] [Materials] - [./elasticity] + [elasticity] type = ADComputeIsotropicElasticityTensorShell youngs_modulus = 2100000 poissons_ratio = 0.0 block = 0 through_thickness_order = SECOND - [../] - [./strain] + [] + [strain] type = ADComputeIncrementalShellStrain block = '0' displacements = 'disp_x disp_y disp_z' rotations = 'rot_x rot_y' thickness = 0.1 through_thickness_order = SECOND - [../] - [./stress] + [] + [stress] type = ADComputeShellStress block = 0 through_thickness_order = SECOND - [../] + [] [] [Postprocessors] - [./disp_z_tip] + [disp_z_tip] type = PointValue point = '1.0 10.0 0.0' variable = disp_z - [../] - [./rot_x_tip] + [] + [rot_x_tip] type = PointValue point = '0.0 10.0 0.0' variable = rot_x - [../] - [./stress_yy_el_0] + [] + [stress_yy_el_0] type = ElementalVariableValue elementid = 0 variable = stress_yy - [../] - [./stress_yy_el_1] + [] + [stress_yy_el_1] type = ElementalVariableValue elementid = 1 variable = stress_yy - [../] - [./stress_yy_el_2] + [] + [stress_yy_el_2] type = ElementalVariableValue elementid = 2 variable = stress_yy - [../] - [./stress_yy_el_3] + [] + [stress_yy_el_3] type = ElementalVariableValue elementid = 3 variable = stress_yy - [../] - [./stress_yz_el_0] + [] + [stress_yz_el_0] type = ElementalVariableValue elementid = 0 variable = stress_yz - [../] - [./stress_yz_el_1] + [] + [stress_yz_el_1] type = ElementalVariableValue elementid = 1 variable = stress_yz - [../] - [./stress_yz_el_2] + [] + [stress_yz_el_2] type = ElementalVariableValue elementid = 2 variable = stress_yz - [../] - [./stress_yz_el_3] + [] + [stress_yz_el_3] type = ElementalVariableValue elementid = 3 variable = stress_yz - [../] + [] [] [Outputs] diff --git a/modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD_2.i b/modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD_2.i index c96b1681c70e..6a7a6acaf983 100644 --- a/modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD_2.i +++ b/modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD_2.i @@ -32,7 +32,7 @@ # y = -0.57735 * (t/2) is 10 Pa at any location along the length of the beam. [Mesh] - [./gen] + [gen] type = GeneratedMeshGenerator dim = 2 nx = 1 @@ -42,112 +42,112 @@ ymin = 0.0 ymax = 10.0 [] - [./rotate] + [rotate] type = TransformGenerator input = gen transform = ROTATE vector_value = '0 90 0' - [../] + [] [] [Variables] - [./disp_x] + [disp_x] order = FIRST family = LAGRANGE - [../] - [./disp_y] + [] + [disp_y] order = FIRST family = LAGRANGE - [../] - [./disp_z] + [] + [disp_z] order = FIRST family = LAGRANGE - [../] - [./rot_x] + [] + [rot_x] order = FIRST family = LAGRANGE - [../] - [./rot_y] + [] + [rot_y] order = FIRST family = LAGRANGE - [../] + [] [] [AuxVariables] - [./stress_zz] + [stress_zz] order = CONSTANT family = MONOMIAL - [../] - [./stress_yz] + [] + [stress_yz] order = CONSTANT family = MONOMIAL - [../] + [] [] [AuxKernels] - [./stress_zz] + [stress_zz] type = RankTwoAux variable = stress_zz rank_two_tensor = global_stress_t_points_0 index_i = 2 index_j = 2 - [../] - [./stress_yz] + [] + [stress_yz] type = RankTwoAux variable = stress_yz rank_two_tensor = global_stress_t_points_0 index_i = 1 index_j = 2 - [../] + [] [] [BCs] - [./fixy1] + [fixy1] type = DirichletBC variable = disp_y boundary = 'bottom' value = 0.0 - [../] - [./fixz1] + [] + [fixz1] type = DirichletBC variable = disp_z boundary = 'bottom' value = 0.0 - [../] - [./fixr1] + [] + [fixr1] type = DirichletBC variable = rot_x boundary = 'bottom' value = 0.0 - [../] - [./fixr2] + [] + [fixr2] type = DirichletBC variable = rot_y boundary = 'bottom' value = 0.0 - [../] - [./fixx1] + [] + [fixx1] type = DirichletBC variable = disp_x boundary = 'bottom' value = 0.0 - [../] + [] [] [NodalKernels] - [./force_y2] + [force_y2] type = ConstantRate variable = disp_y boundary = 'top' rate = 0.5 - [../] + [] [] [Preconditioning] - [./smp] + [smp] type = SMP full = true - [../] + [] [] [Executioner] @@ -163,117 +163,117 @@ [] [Kernels] - [./solid_disp_x] + [solid_disp_x] type = ADStressDivergenceShell block = '0' component = 0 variable = disp_x through_thickness_order = SECOND - [../] - [./solid_disp_y] + [] + [solid_disp_y] type = ADStressDivergenceShell block = '0' component = 1 variable = disp_y through_thickness_order = SECOND - [../] - [./solid_disp_z] + [] + [solid_disp_z] type = ADStressDivergenceShell block = '0' component = 2 variable = disp_z through_thickness_order = SECOND - [../] - [./solid_rot_x] + [] + [solid_rot_x] type = ADStressDivergenceShell block = '0' component = 3 variable = rot_x through_thickness_order = SECOND - [../] - [./solid_rot_y] + [] + [solid_rot_y] type = ADStressDivergenceShell block = '0' component = 4 variable = rot_y through_thickness_order = SECOND - [../] + [] [] [Materials] - [./elasticity] + [elasticity] type = ADComputeIsotropicElasticityTensorShell youngs_modulus = 2100000 poissons_ratio = 0.0 block = 0 through_thickness_order = SECOND - [../] - [./strain] + [] + [strain] type = ADComputeIncrementalShellStrain block = '0' displacements = 'disp_x disp_y disp_z' rotations = 'rot_x rot_y' thickness = 0.1 through_thickness_order = SECOND - [../] - [./stress] + [] + [stress] type = ADComputeShellStress block = 0 through_thickness_order = SECOND - [../] + [] [] [Postprocessors] - [./disp_z_tip] + [disp_z_tip] type = PointValue point = '1.0 0.0 10.0' variable = disp_y - [../] - [./rot_y_tip] + [] + [rot_y_tip] type = PointValue - point = '0.0 0.0 10.0' - variable = rot_y - [../] - [./stress_zz_el_0] + point = '1.0 0.0 10.0' + variable = rot_x + [] + [stress_zz_el_0] type = ElementalVariableValue elementid = 0 variable = stress_zz - [../] - [./stress_zz_el_1] + [] + [stress_zz_el_1] type = ElementalVariableValue elementid = 1 variable = stress_zz - [../] - [./stress_zz_el_2] + [] + [stress_zz_el_2] type = ElementalVariableValue elementid = 2 variable = stress_zz - [../] - [./stress_zz_el_3] + [] + [stress_zz_el_3] type = ElementalVariableValue elementid = 3 variable = stress_zz - [../] - [./stress_yz_el_0] + [] + [stress_yz_el_0] type = ElementalVariableValue elementid = 0 variable = stress_yz - [../] - [./stress_yz_el_1] + [] + [stress_yz_el_1] type = ElementalVariableValue elementid = 1 variable = stress_yz - [../] - [./stress_yz_el_2] + [] + [stress_yz_el_2] type = ElementalVariableValue elementid = 2 variable = stress_yz - [../] - [./stress_yz_el_3] + [] + [stress_yz_el_3] type = ElementalVariableValue elementid = 3 variable = stress_yz - [../] + [] [] [Outputs] diff --git a/modules/solid_mechanics/test/tests/shell/static/gold/beam_bending_moment_AD_2_out.e b/modules/solid_mechanics/test/tests/shell/static/gold/beam_bending_moment_AD_2_out.e index d79a4f49faf6..df1b4a1a9140 100644 Binary files a/modules/solid_mechanics/test/tests/shell/static/gold/beam_bending_moment_AD_2_out.e and b/modules/solid_mechanics/test/tests/shell/static/gold/beam_bending_moment_AD_2_out.e differ diff --git a/modules/solid_mechanics/test/tests/shell/static/gold/inclined_straintest_local_stress_out.e b/modules/solid_mechanics/test/tests/shell/static/gold/inclined_straintest_local_stress_out.e new file mode 100644 index 000000000000..31ee8e922e54 Binary files /dev/null and b/modules/solid_mechanics/test/tests/shell/static/gold/inclined_straintest_local_stress_out.e differ diff --git a/modules/solid_mechanics/test/tests/shell/static/gold/large_strain_m_40_AD_out.e b/modules/solid_mechanics/test/tests/shell/static/gold/large_strain_m_40_AD_out.e index 5fc23e189b09..a3b812c6757d 100644 Binary files a/modules/solid_mechanics/test/tests/shell/static/gold/large_strain_m_40_AD_out.e and b/modules/solid_mechanics/test/tests/shell/static/gold/large_strain_m_40_AD_out.e differ diff --git a/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_local_stress_out.csv b/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_local_stress_out.csv new file mode 100644 index 000000000000..1ce32f81831f --- /dev/null +++ b/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_local_stress_out.csv @@ -0,0 +1,3 @@ +time,disp_x1,disp_x2,disp_y1,disp_y2 +0,0,0,0,0 +1,-0.13469146310401,0,-4.4206019516245e-19,-0.00092841071532347 diff --git a/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_local_stress_out.e b/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_local_stress_out.e new file mode 100644 index 000000000000..667e2b135d6a Binary files /dev/null and b/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_local_stress_out.e differ diff --git a/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_unstructured_out.e b/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_unstructured_out.e index 5c2e019a1d13..70c871a66dfd 100644 Binary files a/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_unstructured_out.e and b/modules/solid_mechanics/test/tests/shell/static/gold/pinched_cylinder_symm_unstructured_out.e differ diff --git a/modules/solid_mechanics/test/tests/shell/static/gold/plate_cantilever_out.e b/modules/solid_mechanics/test/tests/shell/static/gold/plate_cantilever_out.e new file mode 100644 index 000000000000..7afeb0b2bcfd Binary files /dev/null and b/modules/solid_mechanics/test/tests/shell/static/gold/plate_cantilever_out.e differ diff --git a/modules/solid_mechanics/test/tests/shell/static/gold/plate_concentrated_loads_out.e b/modules/solid_mechanics/test/tests/shell/static/gold/plate_concentrated_loads_out.e new file mode 100644 index 000000000000..673c27115e8a Binary files /dev/null and b/modules/solid_mechanics/test/tests/shell/static/gold/plate_concentrated_loads_out.e differ diff --git a/modules/solid_mechanics/test/tests/shell/static/inclined_straintest_local_stress.i b/modules/solid_mechanics/test/tests/shell/static/inclined_straintest_local_stress.i new file mode 100644 index 000000000000..e43251a3232b --- /dev/null +++ b/modules/solid_mechanics/test/tests/shell/static/inclined_straintest_local_stress.i @@ -0,0 +1,280 @@ +# Static test for the inclined shell element. +# A single shell element is oriented at a 45 deg. angle with respect to the Y axis. +# One end of the shell is fixed and an axial deformation to the shell element is +# applied at the other end by resolving the deformation into Y and Z direction. +# The local stresses are computed and stored in aux variables. +# The local stress_22 should be zero (because of plane stress condition). + +[Mesh] + type = FileMesh + file = shell_inclined.e + displacements = 'disp_x disp_y disp_z' +[] + +[Variables] + [disp_x] + [] + [disp_y] + [] + [disp_z] + [] + [rot_x] + [] + [rot_y] + [] +[] + +[AuxVariables] + [stress_00] + order = CONSTANT + family = MONOMIAL + [] + [stress_11] + order = CONSTANT + family = MONOMIAL + [] + [stress_22] + order = CONSTANT + family = MONOMIAL + [] + [stress_01] + order = CONSTANT + family = MONOMIAL + [] + [stress_02] + order = CONSTANT + family = MONOMIAL + [] + [stress_12] + order = CONSTANT + family = MONOMIAL + [] +[] + +[AuxKernels] + [stress_00] + type = RankTwoAux + variable = stress_00 + rank_two_tensor = local_stress_t_points_0 + index_i = 0 + index_j = 0 + execute_on = TIMESTEP_END + [] + [stress_11] + type = RankTwoAux + variable = stress_11 + rank_two_tensor = local_stress_t_points_0 + index_i = 1 + index_j = 1 + execute_on = TIMESTEP_END + [] + [stress_22] + type = RankTwoAux + variable = stress_22 + rank_two_tensor = local_stress_t_points_0 + index_i = 2 + index_j = 2 + execute_on = TIMESTEP_END + [] + [stress_01] + type = RankTwoAux + variable = stress_01 + rank_two_tensor = local_stress_t_points_0 + index_i = 0 + index_j = 1 + execute_on = TIMESTEP_END + [] + [stress_02] + type = RankTwoAux + variable = stress_02 + rank_two_tensor = local_stress_t_points_0 + index_i = 0 + index_j = 2 + execute_on = TIMESTEP_END + [] + [stress_12] + type = RankTwoAux + variable = stress_12 + rank_two_tensor = local_stress_t_points_0 + index_i = 1 + index_j = 2 + execute_on = TIMESTEP_END + [] +[] + +[BCs] + [fixy1] + type = DirichletBC + variable = disp_y + boundary = '0' + value = 0.0 + [] + [fixz1] + type = DirichletBC + variable = disp_z + boundary = '0' + value = 0.0 + [] + [fixr1] + type = DirichletBC + variable = rot_x + boundary = '0' + value = 0.0 + [] + [fixr2] + type = DirichletBC + variable = rot_y + boundary = '0' + value = 0.0 + [] + [fixx1] + type = DirichletBC + variable = disp_x + boundary = '0' + value = 0.0 + [] + [dispz] + type = FunctionDirichletBC + variable = disp_z + boundary = '2' + function = force_function + [] + [dispy] + type = FunctionDirichletBC + variable = disp_y + boundary = '2' + function = force_function + [] +[] + +[Functions] + [force_function] + type = PiecewiseLinear + x = '0.0 1' + y = '0.0 0.33535534' + [] +[] + +[Kernels] + [solid_disp_x] + type = ADStressDivergenceShell + block = '0' + component = 0 + variable = disp_x + through_thickness_order = SECOND + [] + [solid_disp_y] + type = ADStressDivergenceShell + block = '0' + component = 1 + variable = disp_y + through_thickness_order = SECOND + [] + [solid_disp_z] + type = ADStressDivergenceShell + block = '0' + component = 2 + variable = disp_z + through_thickness_order = SECOND + [] + [solid_rot_x] + type = ADStressDivergenceShell + block = '0' + component = 3 + variable = rot_x + through_thickness_order = SECOND + [] + [solid_rot_y] + type = ADStressDivergenceShell + block = '0' + component = 4 + variable = rot_y + through_thickness_order = SECOND + [] +[] + +[Materials] + [elasticity] + type = ADComputeIsotropicElasticityTensorShell + youngs_modulus = 5 + poissons_ratio = 0.0 + block = 0 + through_thickness_order = SECOND + [] + [strain] + type = ADComputeIncrementalShellStrain + block = '0' + displacements = 'disp_x disp_y disp_z' + rotations = 'rot_x rot_y' + thickness = 0.1 + through_thickness_order = SECOND + [] + [stress] + type = ADComputeShellStress + block = 0 + through_thickness_order = SECOND + [] +[] + +[Postprocessors] + [stress_11_el_0] + type = ElementalVariableValue + elementid = 0 + variable = stress_11 + [] + + [stress_12_el_0] + type = ElementalVariableValue + elementid = 0 + variable = stress_12 + [] + + [stress_00_el_0] + type = ElementalVariableValue + elementid = 0 + variable = stress_00 + [] + + [stress_01_el_0] + type = ElementalVariableValue + elementid = 0 + variable = stress_01 + [] + + [stress_02_el_0] + type = ElementalVariableValue + elementid = 0 + variable = stress_02 + [] + + [stress_22_el_0] + type = ElementalVariableValue + elementid = 0 + variable = stress_22 + [] +[] + +[Preconditioning] + [smp] + type = SMP + full = true + [] +[] + +[Executioner] + type = Transient + solve_type = PJFNK + l_tol = 1e-11 + nl_max_its = 15 + nl_rel_tol = 1e-11 + nl_abs_tol = 1e-10 + l_max_its = 20 + dt = 1 + dtmin = 0.01 + timestep_tolerance = 2e-13 + end_time = 1 +[] + +[Outputs] + exodus = true +[] diff --git a/modules/solid_mechanics/test/tests/shell/static/large_strain_m_40_AD.i b/modules/solid_mechanics/test/tests/shell/static/large_strain_m_40_AD.i index e7c69f41fe93..78bb2059534b 100644 --- a/modules/solid_mechanics/test/tests/shell/static/large_strain_m_40_AD.i +++ b/modules/solid_mechanics/test/tests/shell/static/large_strain_m_40_AD.i @@ -14,8 +14,8 @@ # For PL^2/EI = 3, disp_z at tip = 0.6L = 24 m & disp_y at tip = 0.76*L-L = -9.6 m # The FEM solution at tip of cantilever is: -# disp_z = 24.85069 m; relative error = 3.54 % -# disp_y = -9.125937 m; relative error = 5.19 % +# disp_z = 25.2 m; relative error = 5 % +# disp_y = -9.42 m; relative error = 1.87 % [Mesh] type = GeneratedMesh @@ -29,84 +29,84 @@ [] [Variables] - [./disp_x] + [disp_x] order = FIRST family = LAGRANGE - [../] - [./disp_y] + [] + [disp_y] order = FIRST family = LAGRANGE - [../] - [./disp_z] + [] + [disp_z] order = FIRST family = LAGRANGE - [../] - [./rot_x] + [] + [rot_x] order = FIRST family = LAGRANGE - [../] - [./rot_y] + [] + [rot_y] order = FIRST family = LAGRANGE - [../] + [] [] [BCs] - [./fixy1] + [fixy1] type = DirichletBC variable = disp_y boundary = bottom value = 0.0 - [../] - [./fixz1] + [] + [fixz1] type = DirichletBC variable = disp_z boundary = bottom value = 0.0 - [../] - [./fixr1] + [] + [fixr1] type = DirichletBC variable = rot_x boundary = bottom value = 0.0 - [../] - [./fixr2] + [] + [fixr2] type = DirichletBC variable = rot_y boundary = bottom value = 0.0 - [../] - [./fixx1] + [] + [fixx1] type = DirichletBC variable = disp_x boundary = bottom value = 0.0 - [../] + [] [] [NodalKernels] - [./force_y2] + [force_y2] type = UserForcingFunctionNodalKernel variable = disp_z boundary = top function = force_y - [../] + [] [] [Functions] - [./force_y] + [force_y] type = PiecewiseLinear x = '0.0 1.0 3.0' y = '0.0 1.0 1.0' scale_factor = 0.140625 - [../] + [] [] [Preconditioning] - [./smp] + [smp] type = SMP full = true - [../] + [] [] [Executioner] @@ -126,82 +126,82 @@ [] [Kernels] - [./solid_disp_x] + [solid_disp_x] type = ADStressDivergenceShell block = '0' component = 0 variable = disp_x through_thickness_order = SECOND large_strain = true - [../] - [./solid_disp_y] + [] + [solid_disp_y] type = ADStressDivergenceShell block = '0' component = 1 variable = disp_y through_thickness_order = SECOND large_strain = true - [../] - [./solid_disp_z] + [] + [solid_disp_z] type = ADStressDivergenceShell block = '0' component = 2 variable = disp_z through_thickness_order = SECOND large_strain = true - [../] - [./solid_rot_x] + [] + [solid_rot_x] type = ADStressDivergenceShell block = '0' component = 3 variable = rot_x through_thickness_order = SECOND large_strain = true - [../] - [./solid_rot_y] + [] + [solid_rot_y] type = ADStressDivergenceShell block = '0' component = 4 variable = rot_y through_thickness_order = SECOND large_strain = true - [../] + [] [] [Materials] - [./elasticity] + [elasticity] type = ADComputeIsotropicElasticityTensorShell youngs_modulus = 1800 poissons_ratio = 0.0 block = 0 through_thickness_order = SECOND - [../] - [./strain] + [] + [strain] type = ADComputeFiniteShellStrain block = '0' displacements = 'disp_x disp_y disp_z' rotations = 'rot_x rot_y' thickness = 1.0 through_thickness_order = SECOND - [../] - [./stress] + [] + [stress] type = ADComputeShellStress block = 0 through_thickness_order = SECOND - [../] + [] [] [Postprocessors] - [./disp_z2] + [disp_z2] type = PointValue point = '1.0 40.0 0.0' variable = disp_z - [../] - [./disp_y2] + [] + [disp_y2] type = PointValue point = '1.0 40.0 0.0' variable = disp_y - [../] + [] [] [Outputs] diff --git a/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm.i b/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm.i index a3da7c2a1187..ea4b314a7fa1 100644 --- a/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm.i +++ b/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm.i @@ -27,7 +27,6 @@ # 40 x 40 | 0.99 | - # 80 x 80 | 1.01 | - - [Mesh] [mesh] type = FileMeshGenerator @@ -84,13 +83,13 @@ [simply_support_rot_x] type = DirichletBC variable = rot_x - boundary = 'AD BC' + boundary = 'AB' value = 0.0 [] [simply_support_rot_y] type = DirichletBC variable = rot_y - boundary = 'AB' + boundary = 'AD BC' value = 0.0 [] [] diff --git a/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_local_stress.i b/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_local_stress.i new file mode 100644 index 000000000000..97ce7e756028 --- /dev/null +++ b/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_local_stress.i @@ -0,0 +1,461 @@ +# test for displacement of pinched cylinder with user-defined local vectors +# everything is similar to the pinch_cylinder_symm.i, except the local coordinates. +# in the original test the first local axis is '0 0 1' +# in this test, the first local vector is defined by the user : first_local_vector_ref='1 -1 0' +# the given vector by the user is projected on the shell elements +# The rotational BCs are switched in order to get same results. + +# Moreover, axiliary variables are added in this test to visualize local coordinates +# The local stresses, forces and bending moments are also calcualted +# The local stress_22 should be zero for all elements + +[Mesh] + [mesh] + type = FileMeshGenerator + file = cyl_sym_10x10.e + [] +[] + +[Variables] + [disp_x] + order = FIRST + family = LAGRANGE + [] + [disp_y] + order = FIRST + family = LAGRANGE + [] + [disp_z] + order = FIRST + family = LAGRANGE + [] + [rot_x] + order = FIRST + family = LAGRANGE + [] + [rot_y] + order = FIRST + family = LAGRANGE + [] +[] + +[BCs] + [simply_support_x] + type = DirichletBC + variable = disp_x + boundary = 'CD AD' + value = 0.0 + [] + [simply_support_y] + type = DirichletBC + variable = disp_y + boundary = 'CD BC' + value = 0.0 + [] + [simply_support_z] + type = DirichletBC + variable = disp_z + boundary = 'AB' + value = 0.0 + [] + [simply_support_rot_x] + type = DirichletBC + variable = rot_x + boundary = 'AB' + value = 0.0 + [] + [simply_support_rot_y] + type = DirichletBC + variable = rot_y + boundary = 'AD BC' + value = 0.0 + [] +[] + +[DiracKernels] + [point] + type = ConstantPointSource + variable = disp_x + point = '1 0 1' + value = -2.5 # P = 10 + [] +[] + +[AuxVariables] + [stress_00] + order = CONSTANT + family = MONOMIAL + [] + [stress_11] + order = CONSTANT + family = MONOMIAL + [] + [stress_22] + order = CONSTANT + family = MONOMIAL + [] + [stress_01] + order = CONSTANT + family = MONOMIAL + [] + [stress_02] + order = CONSTANT + family = MONOMIAL + [] + [stress_12] + order = CONSTANT + family = MONOMIAL + [] + [force_1] + order = CONSTANT + family = MONOMIAL + [] + [force_2] + order = CONSTANT + family = MONOMIAL + [] + [moment_11] + order = CONSTANT + family = MONOMIAL + [] + [moment_22] + order = CONSTANT + family = MONOMIAL + [] + [moment_12] + order = CONSTANT + family = MONOMIAL + [] + [shear_12] + order = CONSTANT + family = MONOMIAL + [] + [shear_13] + order = CONSTANT + family = MONOMIAL + [] + [shear_23] + order = CONSTANT + family = MONOMIAL + [] + [first_axis_x] + order = CONSTANT + family = MONOMIAL + [] + [first_axis_y] + order = CONSTANT + family = MONOMIAL + [] + [first_axis_z] + order = CONSTANT + family = MONOMIAL + [] + [second_axis_x] + order = CONSTANT + family = MONOMIAL + [] + [second_axis_y] + order = CONSTANT + family = MONOMIAL + [] + [second_axis_z] + order = CONSTANT + family = MONOMIAL + [] + [normal_axis_x] + order = CONSTANT + family = MONOMIAL + [] + [normal_axis_y] + order = CONSTANT + family = MONOMIAL + [] + [normal_axis_z] + order = CONSTANT + family = MONOMIAL + [] +[] + +[AuxKernels] + [stress_00] + type = RankTwoAux + variable = stress_00 + rank_two_tensor = local_stress_t_points_0 + index_i = 0 + index_j = 0 + execute_on = TIMESTEP_END + [] + [stress_11] + type = RankTwoAux + variable = stress_11 + rank_two_tensor = local_stress_t_points_0 + index_i = 1 + index_j = 1 + execute_on = TIMESTEP_END + [] + [stress_22] + type = RankTwoAux + variable = stress_22 + rank_two_tensor = local_stress_t_points_0 + index_i = 2 + index_j = 2 + execute_on = TIMESTEP_END + [] + [stress_01] + type = RankTwoAux + variable = stress_01 + rank_two_tensor = local_stress_t_points_0 + index_i = 0 + index_j = 1 + execute_on = TIMESTEP_END + [] + [stress_02] + type = RankTwoAux + variable = stress_02 + rank_two_tensor = local_stress_t_points_0 + index_i = 0 + index_j = 2 + execute_on = TIMESTEP_END + [] + [stress_12] + type = RankTwoAux + variable = stress_12 + rank_two_tensor = local_stress_t_points_0 + index_i = 1 + index_j = 2 + execute_on = TIMESTEP_END + [] + + [force_1] + type = ShellResultantsAux + variable = force_1 + stress_resultant = axial_force_0 + thickness = 0.01 + through_thickness_order = SECOND + execute_on = TIMESTEP_END + [] + [force_2] + type = ShellResultantsAux + variable = force_2 + stress_resultant = axial_force_1 + thickness = 0.01 + through_thickness_order = SECOND + execute_on = TIMESTEP_END + [] + [moment_11] + type = ShellResultantsAux + variable = moment_11 + stress_resultant = bending_moment_0 + thickness = 0.01 + through_thickness_order = SECOND + execute_on = TIMESTEP_END + [] + [moment_22] + type = ShellResultantsAux + variable = moment_22 + stress_resultant = bending_moment_1 + thickness = 0.01 + through_thickness_order = SECOND + execute_on = TIMESTEP_END + [] + + [moment_12] + type = ShellResultantsAux + variable = moment_12 + stress_resultant = bending_moment_01 + thickness = 0.01 + through_thickness_order = SECOND + execute_on = TIMESTEP_END + [] + [shear_12] + type = ShellResultantsAux + variable = shear_12 + stress_resultant = shear_force_01 + thickness = 0.01 + through_thickness_order = SECOND + execute_on = TIMESTEP_END + [] + [shear_13] + type = ShellResultantsAux + variable = shear_13 + stress_resultant = shear_force_02 + thickness = 0.01 + through_thickness_order = SECOND + execute_on = TIMESTEP_END + [] + [shear_23] + type = ShellResultantsAux + variable = shear_23 + stress_resultant = shear_force_12 + thickness = 0.01 + through_thickness_order = SECOND + execute_on = TIMESTEP_END + [] + [first_axis_x] + type = ShellLocalCoordinatesAux + variable = first_axis_x + property = first_local_vector + component = 0 + [] + [first_axis_y] + type = ShellLocalCoordinatesAux + variable = first_axis_y + property = first_local_vector + component = 1 + [] + [first_axis_z] + type = ShellLocalCoordinatesAux + variable = first_axis_z + property = first_local_vector + component = 2 + [] + + [second_axis_x] + type = ShellLocalCoordinatesAux + variable = second_axis_x + property = second_local_vector + component = 0 + [] + [second_axis_y] + type = ShellLocalCoordinatesAux + variable = second_axis_y + property = second_local_vector + component = 1 + [] + [second_axis_z] + type = ShellLocalCoordinatesAux + variable = second_axis_z + property = second_local_vector + component = 2 + [] + + [normal_axis_x] + type = ShellLocalCoordinatesAux + variable = normal_axis_x + property = normal_local_vector + component = 0 + [] + [normal_axis_y] + type = ShellLocalCoordinatesAux + variable = normal_axis_y + property = normal_local_vector + component = 1 + [] + [normal_axis_z] + type = ShellLocalCoordinatesAux + variable = normal_axis_z + property = normal_local_vector + component = 2 + [] +[] + +[Preconditioning] + [smp] + type = SMP + full = true + [] +[] + +[Executioner] + type = Transient + solve_type = NEWTON + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-pc_type' + petsc_options_value = ' lu' + line_search = 'none' + nl_rel_tol = 1e-10 + nl_abs_tol = 1e-8 + dt = 1.0 + dtmin = 1.0 + end_time = 1.0 +[] + +[Kernels] + [solid_disp_x] + type = ADStressDivergenceShell + block = '100' + component = 0 + variable = disp_x + through_thickness_order = SECOND + [] + [solid_disp_y] + type = ADStressDivergenceShell + block = '100' + component = 1 + variable = disp_y + through_thickness_order = SECOND + [] + [solid_disp_z] + type = ADStressDivergenceShell + block = '100' + component = 2 + variable = disp_z + through_thickness_order = SECOND + [] + [solid_rot_x] + type = ADStressDivergenceShell + block = '100' + component = 3 + variable = rot_x + through_thickness_order = SECOND + [] + [solid_rot_y] + type = ADStressDivergenceShell + block = '100' + component = 4 + variable = rot_y + through_thickness_order = SECOND + [] +[] + +[Materials] + [elasticity] + type = ADComputeIsotropicElasticityTensorShell + youngs_modulus = 1e6 + poissons_ratio = 0.3 + block = '100' + through_thickness_order = SECOND + [] + [strain] + type = ADComputeIncrementalShellStrain + block = '100' + displacements = 'disp_x disp_y disp_z' + rotations = 'rot_x rot_y' + thickness = 0.01 + through_thickness_order = SECOND + reference_first_local_direction = '1 -1 0' + [] + [stress] + type = ADComputeShellStress + block = '100' + through_thickness_order = SECOND + [] +[] + +[Postprocessors] + [disp_x1] + type = PointValue + point = '1 0 1' + variable = disp_x + [] + [disp_y1] + type = PointValue + point = '1 0 1' + variable = disp_y + [] + [disp_x2] + type = PointValue + point = '0 1 1' + variable = disp_x + [] + [disp_y2] + type = PointValue + point = '0 1 1' + variable = disp_y + [] +[] + +[Outputs] + exodus = true + csv = true +[] diff --git a/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_unstructured.i b/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_unstructured.i index 174fdce6ec54..b97a526d80cf 100644 --- a/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_unstructured.i +++ b/modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_unstructured.i @@ -111,19 +111,19 @@ [simply_support_z] type = DirichletBC variable = disp_z - boundary = 'CD AB' + boundary = 'AB' value = 0.0 [] [simply_support_rot_x] type = DirichletBC variable = rot_x - boundary = 'CD BC' + boundary = 'AD BC' value = 0.0 [] [simply_support_rot_y] type = DirichletBC variable = rot_y - boundary = 'CD AD' + boundary = 'AB' value = 0.0 [] [] @@ -211,6 +211,7 @@ rotations = 'rot_x rot_y' thickness = 0.01 through_thickness_order = SECOND + reference_first_local_direction = '0 0 1' [] [stress] type = ADComputeShellStress diff --git a/modules/solid_mechanics/test/tests/shell/static/plate_cantilever.i b/modules/solid_mechanics/test/tests/shell/static/plate_cantilever.i new file mode 100644 index 000000000000..b4c97ffbdebc --- /dev/null +++ b/modules/solid_mechanics/test/tests/shell/static/plate_cantilever.i @@ -0,0 +1,267 @@ +#constant bending of 0.05 applied to the tip of a Plate_Cantilever +#Analytical bending=ML/EI, deflection=ML^2/2EI +#E=200e9, I=bh3/12=2e-4 +#Therefore, analytical solution M22=2e5, uz=0.25 +#Numerical results M22=2e5, uz=0.25 + +[Mesh] + type = GeneratedMesh + dim = 2 + nx = 10 + ny = 1 + xmin = 0.0 + xmax = 10 + zmin = 0.0 + zmax = 1 +[] + +[Variables] + [disp_x] + order = FIRST + family = LAGRANGE + [] + [disp_y] + order = FIRST + family = LAGRANGE + [] + [disp_z] + order = FIRST + family = LAGRANGE + [] + [rot_x] + order = FIRST + family = LAGRANGE + [] + [rot_y] + order = FIRST + family = LAGRANGE + [] +[] + +[BCs] + [symm_left_rot] + type = DirichletBC + variable = rot_y + boundary = 'left' + value = 0.0 + [] + + [simply_support_x] + type = DirichletBC + variable = disp_x + boundary = 'left' + value = 0.0 + [] + [simply_support_y] + type = DirichletBC + variable = disp_z + boundary = 'left' + value = 0.0 + [] + [simply_moment_x] + type = DirichletBC + variable = rot_y + boundary = 'right' + value = 0.05 + [] +[] + +[Preconditioning] + [smp] + type = SMP + full = true + [] +[] + +[Executioner] + type = Transient + solve_type = NEWTON + petsc_options = '-snes_ksp_ew' + petsc_options_iname = '-pc_type -pc_factor_mat_solver_package' + petsc_options_value = ' lu mumps' + line_search = 'none' + nl_rel_tol = 1e-8 + nl_abs_tol = 1e-5 + dt = 1 + dtmin = 1 + end_time = 1. +[] + +[Kernels] + [solid_disp_x] + type = ADStressDivergenceShell + component = 0 + variable = disp_x + through_thickness_order = SECOND + [] + [solid_disp_y] + type = ADStressDivergenceShell + component = 1 + variable = disp_y + through_thickness_order = SECOND + [] + [solid_disp_z] + type = ADStressDivergenceShell + component = 2 + variable = disp_z + through_thickness_order = SECOND + [] + [solid_rot_x] + type = ADStressDivergenceShell + component = 3 + variable = rot_x + through_thickness_order = SECOND + [] + [solid_rot_y] + type = ADStressDivergenceShell + component = 4 + variable = rot_y + through_thickness_order = SECOND + [] +[] + +[Materials] + [elasticity] + type = ADComputeIsotropicElasticityTensorShell + youngs_modulus = 2e11 + poissons_ratio = 0.0 + through_thickness_order = SECOND + [] + [strain] + type = ADComputeIncrementalShellStrain + displacements = 'disp_x disp_y disp_z' + rotations = 'rot_x rot_y' + thickness = 0.133887 + through_thickness_order = SECOND + [] + [stress] + type = ADComputeShellStress + through_thickness_order = SECOND + [] +[] + +[Postprocessors] + [disp_z2] + type = PointValue + point = '10.0 0.0 0.0' + variable = disp_z + [] +[] + +[AuxVariables] + + [moment_22] + order = CONSTANT + family = MONOMIAL + [] + + [first_axis_x] + order = CONSTANT + family = MONOMIAL + [] + [first_axis_y] + order = CONSTANT + family = MONOMIAL + [] + [first_axis_z] + order = CONSTANT + family = MONOMIAL + [] + [second_axis_x] + order = CONSTANT + family = MONOMIAL + [] + [second_axis_y] + order = CONSTANT + family = MONOMIAL + [] + [second_axis_z] + order = CONSTANT + family = MONOMIAL + [] + [normal_axis_x] + order = CONSTANT + family = MONOMIAL + [] + [normal_axis_y] + order = CONSTANT + family = MONOMIAL + [] + [normal_axis_z] + order = CONSTANT + family = MONOMIAL + [] +[] + +[AuxKernels] + + [moment_22] + type = ShellResultantsAux + variable = moment_22 + stress_resultant = bending_moment_1 + thickness = 0.133887 + through_thickness_order = SECOND + execute_on = TIMESTEP_END + [] + + [first_axis_x] + type = ShellLocalCoordinatesAux + variable = first_axis_x + property = first_local_vector + component = 0 + [] + [first_axis_y] + type = ShellLocalCoordinatesAux + variable = first_axis_y + property = first_local_vector + component = 1 + [] + [first_axis_z] + type = ShellLocalCoordinatesAux + variable = first_axis_z + property = first_local_vector + component = 2 + [] + + [second_axis_x] + type = ShellLocalCoordinatesAux + variable = second_axis_x + property = second_local_vector + component = 0 + [] + [second_axis_y] + type = ShellLocalCoordinatesAux + variable = second_axis_y + property = second_local_vector + component = 1 + [] + [second_axis_z] + type = ShellLocalCoordinatesAux + variable = second_axis_z + property = second_local_vector + component = 2 + [] + + [normal_axis_x] + type = ShellLocalCoordinatesAux + variable = normal_axis_x + property = normal_local_vector + component = 0 + [] + [normal_axis_y] + type = ShellLocalCoordinatesAux + variable = normal_axis_y + property = normal_local_vector + component = 1 + [] + [normal_axis_z] + type = ShellLocalCoordinatesAux + variable = normal_axis_z + property = normal_local_vector + component = 2 + [] +[] + +[Outputs] + exodus = true +[] diff --git a/modules/solid_mechanics/test/tests/shell/static/plate_concentrated_loads.i b/modules/solid_mechanics/test/tests/shell/static/plate_concentrated_loads.i new file mode 100644 index 000000000000..d406c7d39c5c --- /dev/null +++ b/modules/solid_mechanics/test/tests/shell/static/plate_concentrated_loads.i @@ -0,0 +1,298 @@ +# A simply supported plate with a length of 9m and width of 1m is loaded by two equal concentrated loads (F=10000 N/m) +# The concentrated loads are symmetrically applied at x=3 and x=6 + +# Analytical solution: maximum diplacement at the center= 6.469e-3 +# Numerical model: maximum diplacement at the center=6.436e-3 + +# Analytical solution: maximum bending moment (m22) at the center =30000 +# Numerical model: maximum bending moment (m22) at the center =30000 + +# Analytical solution: out of plane shear force (q13) for 0