- componentAn integer corresponding to the degree of freedom this kernel acts on. (0 for disp_x, 1 for disp_y, 2 for disp_z, 3 for rot_x, 4 for rot_y)
C++ Type:unsigned int
Controllable:No
Description:An integer corresponding to the degree of freedom this kernel acts on. (0 for disp_x, 1 for disp_y, 2 for disp_z, 3 for rot_x, 4 for rot_y)
- through_thickness_orderQuadrature order in out of plane direction
C++ Type:std::string
Controllable:No
Description:Quadrature order in out of plane direction
- variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on
ADStressDivergenceShell
Quasi-static stress divergence kernel for Shell element
Description
The ADStressDivergenceShell kernel calculates the contribution to the residual from the stress divergence for shell elements. Forward mode automatic differentiation is used to compute an exact Jacobian. Please refer to ShellElements for details.
Example Input File syntax
Input Parameters
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- displacementsThe displacements
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacements
- large_strainFalseSet to true to turn on finite strain calculations.
Default:False
C++ Type:bool
Controllable:No
Description:Set to true to turn on finite strain calculations.
- matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)
Optional Parameters
- absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
- extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
- extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
- matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
- vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill
Contribution To Tagged Field Data Parameters
- control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
- diag_save_inThe name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- save_inThe name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
- use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Advanced Parameters
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
- use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Material Property Retrieval Parameters
Input Files
- (modules/solid_mechanics/test/tests/shell/static/finite_straintest.i)
- (modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment_free_orientation_inclined_hht.i)
- (modules/solid_mechanics/test/tests/shell/static/inclined_straintest.i)
- (modules/solid_mechanics/test/tests/shell/static/tank_shell.i)
- (modules/solid_mechanics/test/tests/shell/static/plate_bending.i)
- (modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment_free.i)
- (modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD_2.i)
- (modules/solid_mechanics/test/tests/shell/static/tapered.i)
- (modules/solid_mechanics/test/tests/shell/static/plate_bending2.i)
- (modules/solid_mechanics/test/tests/shell/static/tank_shell_rotated.i)
- (modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD.i)
- (modules/solid_mechanics/test/tests/shell/static/large_strain_m_40_AD.i)
- (modules/solid_mechanics/test/tests/shell/static/scordelis_lo_roof_shell.i)
- (modules/solid_mechanics/test/tests/shell/static/clamped_plate_flat.i)
- (modules/solid_mechanics/test/tests/shell/static/plate_cantilever.i)
- (modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment.i)
- (modules/solid_mechanics/test/tests/shell/static/pressure_error.i)
- (modules/solid_mechanics/test/tests/shell/static/inclined_straintest_local_stress.i)
- (modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm.i)
- (modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment_free_orientation_inclined.i)
- (modules/solid_mechanics/test/tests/shell/static/straintest_shear.i)
- (modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_unstructured.i)
- (modules/solid_mechanics/test/tests/shell/static/straintest.i)
- (modules/solid_mechanics/test/tests/shell/static/plate_concentrated_loads.i)
- (modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_local_stress.i)
(modules/solid_mechanics/test/tests/shell/static/finite_straintest.i)
# Test for the axial stress and strain output for single shell element
# for 2D planar shell with uniform mesh.
# A single shell 1 mm x 1 mm element having Young's Modulus of 5 N/mm^2
# and poissons ratio of 0 is fixed at the left end and
# an axial displacement of 0.2 mm is applied at the right.
# Theoretical value of axial stress and strain are 1 N/mm^2 and 0.2.
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 1.0
[]
[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
[]
[]
[AuxVariables]
[stress_xx]
order = CONSTANT
family = MONOMIAL
[]
[strain_xx]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_xx]
type = RankTwoAux
variable = stress_xx
rank_two_tensor = global_stress_t_points_1
index_i = 0
index_j = 0
[]
[strain_xx]
type = RankTwoAux
variable = strain_xx
rank_two_tensor = total_global_strain_t_points_1
index_i = 0
index_j = 0
[]
[]
[BCs]
[fixx]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[]
[fixy]
type = DirichletBC
variable = disp_y
boundary = left
value = 0.0
[]
[fixz]
type = DirichletBC
variable = disp_z
boundary = left
value = 0.0
[]
[fixr1]
type = DirichletBC
variable = rot_x
boundary = left
value = 0.0
[]
[fixr2]
type = DirichletBC
variable = rot_y
boundary = left
value = 0.0
[]
[disp]
type = FunctionDirichletBC
variable = disp_x
boundary = 'right'
function = displacement
[]
[]
[Functions]
[displacement]
type = PiecewiseLinear
x = '0.0 1.0'
y = '0.0 0.2'
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
line_search = 'none'
nl_rel_tol = 1e-12
nl_abs_tol = 1e-14
dt = 1
dtmin = 1
end_time = 1
[]
[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.0
poissons_ratio = 0.0
block = 0
through_thickness_order = SECOND
[]
[strain]
type = ADComputeFiniteShellStrain
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]
[disp_x]
type = PointValue
point = '0.5 0.0 0.0'
variable = disp_z
[]
[stress_xx_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_xx
[]
[strain_xx_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_xx
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment_free_orientation_inclined_hht.i)
# Test to verify the fundamental natural frequency of a one element ADComputeShellStress
# BCs: Clamped on one end, free on others.
# Initial perturbation applied to edge of the beam. After that, the shell vibrates freely.
#
# Results have been compared for various thicknesses with the following approximate Results
# (Moose results were obtained with 8 elements along the length)
# Thickness = 0.1. Reference freq: 10.785 Hz, Moose freq: 10.612 Hz
# Thickness = 0.05. Reference freq: 5.393 Hz, Moose freq: 5.335 Hz
# Thickness = 0.025. Reference freq: 2.696 Hz, Moose freq: 2.660 Hz
#
# Reference values have been obtained from Robert Blevins, "Formulas for Dynamics, Acoustics and Vibration",
# Table 5.3 case 11. Formula looks like: f = lambda^2/(2*pi*a^2) * sqrt(E*h^2/(12*(1-nu*nu))), where lambda
# changes as a function of shell dimensions.
# This test uses one single element for speed reasons.
# Here, the shell, instead of being on the XY plane, is oriented at a 45 deg. angle
# with respect to the Y axis.
[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_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz]
order = CONSTANT
family = MONOMIAL
[../]
# aux variables for dynamics
[./vel_x]
[../]
[./vel_y]
[../]
[./vel_z]
[../]
[./accel_x]
[../]
[./accel_y]
[../]
[./accel_z]
[../]
[./rot_vel_x]
[../]
[./rot_vel_y]
[../]
[./rot_accel_x]
[../]
[./rot_accel_y]
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
variable = stress_yy
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 1
[../]
[./stress_yz]
type = RankTwoAux
variable = stress_yz
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 2
[../]
# Kernels for dynamics
[./accel_x]
type = NewmarkAccelAux
variable = accel_x
displacement = disp_x
velocity = vel_x
beta = 0.25
execute_on = timestep_end
[../]
[./vel_x]
type = NewmarkVelAux
variable = vel_x
acceleration = accel_x
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_y]
type = NewmarkAccelAux
variable = accel_y
displacement = disp_y
velocity = vel_y
beta = 0.25
execute_on = timestep_end
[../]
[./vel_y]
type = NewmarkVelAux
variable = vel_y
acceleration = accel_y
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_z]
type = NewmarkAccelAux
variable = accel_z
displacement = disp_z
velocity = vel_z
beta = 0.25
execute_on = timestep_end
[../]
[./vel_z]
type = NewmarkVelAux
variable = vel_z
acceleration = accel_z
gamma = 0.5
execute_on = timestep_end
[../]
[./rot_accel_x]
type = NewmarkAccelAux
variable = rot_accel_x
displacement = rot_x
velocity = rot_vel_x
beta = 0.25
execute_on = timestep_end
[../]
[./rot_vel_x]
type = NewmarkVelAux
variable = rot_vel_x
acceleration = rot_accel_x
gamma = 0.5
execute_on = timestep_end
[../]
[./rot_accel_y]
type = NewmarkAccelAux
variable = rot_accel_y
displacement = rot_y
velocity = rot_vel_y
beta = 0.25
execute_on = timestep_end
[../]
[./rot_vel_y]
type = NewmarkVelAux
variable = rot_vel_y
acceleration = rot_accel_y
gamma = 0.5
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
[../]
[]
[Functions]
[./force_function]
type = PiecewiseLinear
x = '0.0 0.01 0.15 10.0'
y = '0.0 0.01 0.0 0.0'
[../]
[]
[NodalKernels]
[./force_y2]
type = UserForcingFunctorNodalKernel
variable = disp_z
boundary = '2'
functor = force_function
[../]
[]
[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
[../]
[./inertial_force_x]
type = ADInertialForceShell
use_displaced_mesh = true
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 0
variable = disp_x
thickness = 0.1
eta = 0.0
alpha = 0.0
[../]
[./inertial_force_y]
type = ADInertialForceShell
use_displaced_mesh = true
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 1
variable = disp_y
thickness = 0.1
eta = 0.0
alpha = 0.0
[../]
[./inertial_force_z]
type = ADInertialForceShell
use_displaced_mesh = true
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 2
variable = disp_z
thickness = 0.1
eta = 0.0
alpha = 0.0
[../]
[./inertial_force_rot_x]
type = ADInertialForceShell
use_displaced_mesh = true
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 3
variable = rot_x
thickness = 0.1
eta = 0.0
alpha = 0.0
[../]
[./inertial_force_rot_y]
type = ADInertialForceShell
use_displaced_mesh = true
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 4
variable = rot_y
thickness = 0.1
eta = 0.0
alpha = 0.0
[../]
[]
[Materials]
[./elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 2100000
poissons_ratio = 0.3
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
[../]
[./density]
type = GenericConstantMaterial
block = 0
prop_names = 'density'
prop_values = '1.0'
[../]
[]
[Postprocessors]
[./disp_z_tip]
type = PointValue
point = '0.0 1.06 1.06'
variable = disp_z
[../]
[./rot_x_tip]
type = PointValue
point = '0.0 1.06 1.06'
variable = rot_x
[../]
[./stress_yy_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yy
[../]
[./stress_yy_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yy
[../]
[./stress_yy_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yy
[../]
[./stress_yy_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yy
[../]
[./stress_yz_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yz
[../]
[./stress_yz_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yz
[../]
[./stress_yz_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yz
[../]
[./stress_yz_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yz
[../]
[]
[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 = 0.005
dtmin = 0.005
timestep_tolerance = 2e-13
end_time = 0.5
[./TimeIntegrator]
type = NewmarkBeta
beta = 0.25
gamma = 0.5
[../]
[]
[Outputs]
perf_graph = true
csv = true
[]
(modules/solid_mechanics/test/tests/shell/static/inclined_straintest.i)
# 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 stress and strain result in the global orientation when transformed to
# the shell oriention gives the correct value of the axial stress and strain.
[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_xx]
order = CONSTANT
family = MONOMIAL
[]
[strain_xx]
order = CONSTANT
family = MONOMIAL
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[stress_xy]
order = CONSTANT
family = MONOMIAL
[]
[strain_xy]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[strain_yz]
order = CONSTANT
family = MONOMIAL
[]
[stress_xz]
order = CONSTANT
family = MONOMIAL
[]
[strain_xz]
order = CONSTANT
family = MONOMIAL
[]
[stress_zz]
order = CONSTANT
family = MONOMIAL
[]
[strain_zz]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_xx]
type = RankTwoAux
variable = stress_xx
selected_qp = 0
rank_two_tensor = global_stress_t_points_0
index_i = 0
index_j = 0
[]
[strain_xx]
type = RankTwoAux
variable = strain_xx
rank_two_tensor = total_global_strain_t_points_0
selected_qp = 0
index_i = 0
index_j = 0
[]
[stress_yy]
type = RankTwoAux
variable = stress_yy
rank_two_tensor = global_stress_t_points_0
selected_qp = 0
index_i = 1
index_j = 1
[]
[strain_yy]
type = RankTwoAux
variable = strain_yy
rank_two_tensor = total_global_strain_t_points_0
selected_qp = 0
index_i = 1
index_j = 1
[]
[stress_xy]
type = RankTwoAux
variable = stress_xy
rank_two_tensor = global_stress_t_points_0
selected_qp = 0
index_i = 0
index_j = 1
[]
[strain_xy]
type = RankTwoAux
variable = strain_xy
rank_two_tensor = total_global_strain_t_points_0
selected_qp = 0
index_i = 0
index_j = 1
[]
[stress_yz]
type = RankTwoAux
variable = stress_yz
rank_two_tensor = global_stress_t_points_0
selected_qp = 0
index_i = 1
index_j = 2
[]
[strain_yz]
type = RankTwoAux
variable = strain_yz
rank_two_tensor = total_global_strain_t_points_0
selected_qp = 0
index_i = 1
index_j = 2
[]
[stress_xz]
type = RankTwoAux
variable = stress_xz
rank_two_tensor = global_stress_t_points_0
selected_qp = 0
index_i = 0
index_j = 2
[]
[strain_xz]
type = RankTwoAux
variable = strain_xz
rank_two_tensor = total_global_strain_t_points_0
selected_qp = 0
index_i = 0
index_j = 2
[]
[stress_zz]
type = RankTwoAux
variable = stress_zz
rank_two_tensor = global_stress_t_points_0
selected_qp = 0
index_i = 2
index_j = 2
[]
[strain_zz]
type = RankTwoAux
variable = strain_zz
rank_two_tensor = total_global_strain_t_points_0
selected_qp = 0
index_i = 2
index_j = 2
[]
[]
[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_yy_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yy
[]
[strain_yy_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_yy
[]
[stress_yz_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yz
[]
[strain_yz_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_yz
[]
[stress_xx_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_xx
[]
[strain_xx_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_xx
[]
[stress_xy_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_xy
[]
[strain_xy_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_xy
[]
[stress_xz_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_xz
[]
[strain_xz_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_xz
[]
[stress_zz_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_zz
[]
[strain_zz_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_zz
[]
[]
[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
[]
(modules/solid_mechanics/test/tests/shell/static/tank_shell.i)
# Test for distributed load on shell (pressure)
# A long cylindrical tank (3m) with a wall thickness of t=0.03 m and a radius of 0.5m is subjected to an internal presure of p=40MPa.
# The lower part of the cylinder is constrained in the x-y-z directions
# Theorically, assuming a thin_walled cylinder t/r <0.1, the hoop stress is sigma_t=p*r/t
# Therefore, in-plane force in the circumference of the cylinder is F=sigma_t*t= p*r=0.5*40=20 MN (independent of material properties of the shell)
# Analytical solution for the radial displacement : u=p*r^2/(E*t)=0.00167 m
# We check the axial_force_1 at the upper part of the cylinder (far from the lower boundary to avoid boundary effects)
# The numerical modeling results in axial_force_1 =19.882 MPa (0.6% error) and radial displacement u=0.00165 (1.1% error)
[Mesh]
[gmg]
type = FileMeshGenerator
file = tank_shell.msh
[]
[]
[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
[]
[]
[AuxVariables]
[axial_force_1]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[axial_force_1]
type = ShellResultantsAux
variable = axial_force_1
stress_resultant = axial_force_1
thickness = 0.03
through_thickness_order = SECOND
execute_on = TIMESTEP_END
[]
[]
[BCs]
[simply_support_x]
type = DirichletBC
variable = disp_x
boundary = 'lower_circle'
value = 0.0
[]
[simply_support_y]
type = DirichletBC
variable = disp_y
boundary = 'lower_circle'
value = 0.0
[]
[simply_support_z]
type = DirichletBC
variable = disp_z
boundary = 'lower_circle'
value = 0.0
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-snes_ksp_ew'
# best overall
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
line_search = 'none'
nl_rel_tol = 1e-12
nl_abs_tol = 1e-13
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
[]
[load_x]
type = ADDistributedLoadShell
function = '40'
variable = disp_x
project_load_to_normal = true
displacements = 'disp_x disp_y disp_z'
[]
[load_y]
type = ADDistributedLoadShell
function = '40'
variable = disp_y
project_load_to_normal = true
displacements = 'disp_x disp_y disp_z'
[]
[load_z]
type = ADDistributedLoadShell
function = '40'
variable = disp_z
project_load_to_normal = true
displacements = 'disp_x disp_y disp_z'
[]
[]
[Materials]
[elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 2e5
poissons_ratio = 0.3
through_thickness_order = SECOND
[]
[strain]
type = ADComputeIncrementalShellStrain
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
thickness = 0.03
reference_first_local_direction = ' 0 0 1'
through_thickness_order = SECOND
[]
[stress]
type = ADComputeShellStress
through_thickness_order = SECOND
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/plate_bending.i)
# Test for simply supported plate under uniform pressure
# One quarter of a 50 m x 50 m x 1m plate is modeled in this test.
# Pressure loading is applied on the top surface using nodal forces
# of magnitude -10 N on all nodes. This corresponds to a pressure (q) of
# -10.816 N/m^2.
# The FEM solution at (0,0), which is at the center of the full plate
# is -2.997084e-03 m.
# The analytical solution for displacement at center of plate obtained
# using a thin plate assumption for a square plate is
# w = 16 q a^4/(D*pi^6) \sum_{m = 1,3,5, ..}^\inf \sum_{n = 1,3,5, ..}^\inf (-1)^{(m+n-2)/2}/(mn*(m^2+n^2)^2)
# The above solution is the Navier's series solution from the "Theory of plates
# and shells" by Timoshenko and Woinowsky-Krieger (1959).
# where a = 50 m, q = -10.816 N/m^2 and D = E/(12(1-v^2))
# The analytical series solution converges to 2.998535904e-03 m
# when the first 16 terms of the series are considered (i.e., until
# m & n = 7).
# The resulting relative error between FEM and analytical solution is
# 0.048%.
[Mesh]
[./gmg]
type = GeneratedMeshGenerator
dim = 2
nx = 25
ny = 25
xmin = 0.0
xmax = 25.0
ymin = 0.0
ymax = 25.0
[../]
[./allnodes]
type = BoundingBoxNodeSetGenerator
input = gmg
bottom_left = '0.0 0.0 0.0'
top_right = '25.0 25.0 0.0'
new_boundary = 101
[../]
[]
[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
[../]
[./symm_bottom_rot]
type = DirichletBC
variable = rot_x
boundary = bottom
value = 0.0
[../]
[./simply_support_x]
type = DirichletBC
variable = disp_x
boundary = 'right top bottom left'
value = 0.0
[../]
[./simply_support_y]
type = DirichletBC
variable = disp_y
boundary = 'right top bottom left'
value = 0.0
[../]
[./simply_support_z]
type = DirichletBC
variable = disp_z
boundary = 'right top'
value = 0.0
[../]
[]
[NodalKernels]
[./force_y2]
type = ConstantRate
variable = disp_z
boundary = 101
rate = -10.0
[../]
[]
[Preconditioning]
[./smp]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = NEWTON
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 = '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 = 1e9
poissons_ratio = 0.3
block = 0
through_thickness_order = SECOND
[../]
[./strain]
type = ADComputeIncrementalShellStrain
block = '0'
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
thickness = 1.0
through_thickness_order = SECOND
[../]
[./stress]
type = ADComputeShellStress
block = 0
through_thickness_order = SECOND
[../]
[]
[Postprocessors]
[./disp_z2]
type = PointValue
point = '0.0 0.0 0.0'
variable = disp_z
[../]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment_free.i)
# Test to verify the fundamental natural frequency of a one element ADComputeShellStress
# BCs: Clamped on one end, free on others.
# Initial perturbation applied to edge of the beam. After that, the shell vibrates freely.
#
# Results have been compared for various thicknesses with the following approximate Results
# (Moose results were obtained with 8 elements along the length)
# Thickness = 0.1. Reference freq: 10.785 Hz, Moose freq: 10.612 Hz
# Thickness = 0.05. Reference freq: 5.393 Hz, Moose freq: 5.335 Hz
# Thickness = 0.025. Reference freq: 2.696 Hz, Moose freq: 2.660 Hz
#
# Reference values have been obtained from Robert Blevins, "Formulas for Dynamics, Acoustics and Vibration",
# Table 5.3 case 11. Formula looks like: f = lambda^2/(2*pi*a^2) * sqrt(E*h^2/(12*(1-nu*nu))), where lambda
# changes as a function of shell dimensions.
# This test uses one single element for speed reasons.
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1 # 1
ny = 1# 4
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 1.5
displacements = 'disp_x disp_y disp_z'
[]
[Variables]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./rot_x]
[../]
[./rot_y]
[../]
[]
[AuxVariables]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz]
order = CONSTANT
family = MONOMIAL
[../]
# aux variables for dynamics
[./vel_x]
[../]
[./vel_y]
[../]
[./vel_z]
[../]
[./accel_x]
[../]
[./accel_y]
[../]
[./accel_z]
[../]
[./rot_vel_x]
[../]
[./rot_vel_y]
[../]
[./rot_accel_x]
[../]
[./rot_accel_y]
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
variable = stress_yy
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 1
[../]
[./stress_yz]
type = RankTwoAux
variable = stress_yz
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 2
[../]
# Kernels for dynamics
[./accel_x]
type = NewmarkAccelAux
variable = accel_x
displacement = disp_x
velocity = vel_x
beta = 0.25
execute_on = timestep_end
[../]
[./vel_x]
type = NewmarkVelAux
variable = vel_x
acceleration = accel_x
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_y]
type = NewmarkAccelAux
variable = accel_y
displacement = disp_y
velocity = vel_y
beta = 0.25
execute_on = timestep_end
[../]
[./vel_y]
type = NewmarkVelAux
variable = vel_y
acceleration = accel_y
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_z]
type = NewmarkAccelAux
variable = accel_z
displacement = disp_z
velocity = vel_z
beta = 0.25
execute_on = timestep_end
[../]
[./vel_z]
type = NewmarkVelAux
variable = vel_z
acceleration = accel_z
gamma = 0.5
execute_on = timestep_end
[../]
[./rot_accel_x]
type = NewmarkAccelAux
variable = rot_accel_x
displacement = rot_x
velocity = rot_vel_x
beta = 0.25
execute_on = timestep_end
[../]
[./rot_vel_x]
type = NewmarkVelAux
variable = rot_vel_x
acceleration = rot_accel_x
gamma = 0.5
execute_on = timestep_end
[../]
[./rot_accel_y]
type = NewmarkAccelAux
variable = rot_accel_y
displacement = rot_y
velocity = rot_vel_y
beta = 0.25
execute_on = timestep_end
[../]
[./rot_vel_y]
type = NewmarkVelAux
variable = rot_vel_y
acceleration = rot_accel_y
gamma = 0.5
execute_on = timestep_end
[../]
[]
[BCs]
[./fixy1]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[../]
[./fixz1]
type = DirichletBC
variable = disp_z
boundary = 'bottom'
value = 0.0
[../]
[./fixr1]
type = DirichletBC
variable = rot_x
boundary = 'bottom'
value = 0.0
[../]
[./fixr2]
type = DirichletBC
variable = rot_y
boundary = 'bottom'
value = 0.0
[../]
[./fixx1]
type = DirichletBC
variable = disp_x
boundary = 'bottom'
value = 0.0
[../]
[]
[Functions]
[./force_function]
type = PiecewiseLinear
x = '0.0 0.01 0.15 10.0'
y = '0.0 0.01 0.0 0.0'
[../]
[]
[NodalKernels]
[./force_z2]
type = UserForcingFunctorNodalKernel
variable = disp_z
boundary = 'top'
functor = force_function
[../]
[]
[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
[../]
[./inertial_force_x]
type = ADInertialForceShell
# use_displaced_mesh = true
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 0
variable = disp_x
thickness = 0.1
[../]
[./inertial_force_y]
type = ADInertialForceShell
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 1
variable = disp_y
thickness = 0.1
[../]
[./inertial_force_z]
type = ADInertialForceShell
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 2
variable = disp_z
thickness = 0.1
[../]
[./inertial_force_rot_x]
type = ADInertialForceShell
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 3
variable = rot_x
thickness = 0.1
[../]
[./inertial_force_rot_y]
type = ADInertialForceShell
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 4
variable = rot_y
thickness = 0.1
[../]
[]
[Materials]
[./elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 2100000
poissons_ratio = 0.3
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
[../]
[./density]
type = GenericConstantMaterial
block = 0
prop_names = 'density'
prop_values = '1.0'
[../]
[]
[Postprocessors]
[./disp_z_tip]
type = PointValue
point = '1.0 1.0 0.0'
variable = disp_z
[../]
[./rot_x_tip]
type = PointValue
point = '0.0 1.0 0.0'
variable = rot_x
[../]
[./stress_yy_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yy
[../]
[./stress_yy_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yy
[../]
[./stress_yy_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yy
[../]
[./stress_yy_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yy
[../]
[./stress_yz_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yz
[../]
[./stress_yz_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yz
[../]
[./stress_yz_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yz
[../]
[./stress_yz_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yz
[../]
[]
[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 = 0.005
dtmin = 0.005
timestep_tolerance = 2e-13
end_time = 0.5
[./TimeIntegrator]
type = NewmarkBeta
beta = 0.25
gamma = 0.5
[../]
[]
[Outputs]
perf_graph = true
csv = true
[]
(modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD_2.i)
# Test that models bending of a rotated cantilever beam using shell elements
# A cantilever beam of length 10 m (in Z direction) and cross-section
# 1 m x 0.1 m is modeled using 4 shell elements placed along the length
# (Figure 6a from Dvorkin and Bathe, 1984). All displacements and
# X rotations are fixed on the bottom boundary. E = 2100000 and v = 0.0.
# A load of 0.5 N (in the Y direction) is applied at each node on the top
# boundary resulting in a total load of 1 N.
# The analytical solution for displacement at tip using small strain/rotations # is PL^3/3EI + PL/AG = 1.90485714 m
# The FEM solution using 4 shell elements is 1.875095 m with a relative error
# of 1.5%.
# Similarly, the analytical solution for slope at tip is PL^2/2EI = 0.285714286
# The FEM solution is 0.2857143 and the relative error is 5e-6%.
# The stress_zz for the four elements at y = -0.57735 * (t/2) (first qp below mid-surface of shell) are:
# 3031.089 Pa, 2165.064 Pa, 1299.038 Pa and 433.0127 Pa.
# Note the above values are the average stresses in each element.
# Analytically, stress_zz decreases linearly from z = 0 to z = 10 m.
# The maximum value of stress_zz at z = 0 is My/I = PL * 0.57735*(t/2)/I = 3464.1 Pa
# Therefore, the analytical value of stress at y = -0.57735 * (t/2) at the mid-point
# of the four elements are:
# 3031.0875 Pa, 2165.0625 Pa, 1299.0375 Pa ,433.0125 Pa
# The relative error in stress_zz is in the order of 5e-5%.
# The stress_yz at y = -0.57735 * (t/2) at all four elements from the simulation is 10 Pa.
# The analytical solution for the shear stress is: V/2/I *((t^2)/4 - y^2), where the shear force (V)
# is 1 N at any z along the length of the beam. Therefore, the analytical shear stress at
# y = -0.57735 * (t/2) is 10 Pa at any location along the length of the beam.
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 2
nx = 1
ny = 4
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 10.0
[]
[rotate]
type = TransformGenerator
input = gen
transform = ROTATE
vector_value = '0 90 0'
[]
[]
[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
[]
[]
[AuxVariables]
[stress_zz]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_zz]
type = RankTwoAux
variable = stress_zz
rank_two_tensor = global_stress_t_points_0
index_i = 2
index_j = 2
[]
[stress_yz]
type = RankTwoAux
variable = stress_yz
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 2
[]
[]
[BCs]
[fixy1]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[fixz1]
type = DirichletBC
variable = disp_z
boundary = 'bottom'
value = 0.0
[]
[fixr1]
type = DirichletBC
variable = rot_x
boundary = 'bottom'
value = 0.0
[]
[fixr2]
type = DirichletBC
variable = rot_y
boundary = 'bottom'
value = 0.0
[]
[fixx1]
type = DirichletBC
variable = disp_x
boundary = 'bottom'
value = 0.0
[]
[]
[NodalKernels]
[force_y2]
type = ConstantRate
variable = disp_y
boundary = 'top'
rate = 0.5
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
nl_max_its = 2
nl_rel_tol = 1e-10
nl_abs_tol = 5e-4
dt = 1
dtmin = 1
end_time = 1
[]
[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 = 2100000
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]
[disp_z_tip]
type = PointValue
point = '1.0 0.0 10.0'
variable = disp_y
[]
[rot_y_tip]
type = PointValue
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]
type = ElementalVariableValue
elementid = 1
variable = stress_zz
[]
[stress_zz_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_zz
[]
[stress_zz_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_zz
[]
[stress_yz_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yz
[]
[stress_yz_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yz
[]
[stress_yz_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yz
[]
[stress_yz_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yz
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/tapered.i)
# Test for the stress and strain output for tapered shell elements.
# A tapered beam is represented with shell elements in XY plane
# having Young's Modulus of 210000 and poissons ratio of 0.3.
# The displacement in X direction is constrained in the left end and the
# displacement of center node of the left end is constrained in Y direction.
# A uniform displacement is applied at the right end.
# The problem is symmetric about Y-axis and the results are symmetric.
[Mesh]
[input]
type = FileMeshGenerator
file = taperedmesh.e
[]
[]
[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_10]
order = CONSTANT
family = MONOMIAL
[]
[stress_02]
order = CONSTANT
family = MONOMIAL
[]
[stress_20]
order = CONSTANT
family = MONOMIAL
[]
[stress_12]
order = CONSTANT
family = MONOMIAL
[]
[stress_21]
order = CONSTANT
family = MONOMIAL
[]
[strain_00]
order = CONSTANT
family = MONOMIAL
[]
[strain_11]
order = CONSTANT
family = MONOMIAL
[]
[strain_22]
order = CONSTANT
family = MONOMIAL
[]
[strain_01]
order = CONSTANT
family = MONOMIAL
[]
[strain_10]
order = CONSTANT
family = MONOMIAL
[]
[strain_02]
order = CONSTANT
family = MONOMIAL
[]
[strain_20]
order = CONSTANT
family = MONOMIAL
[]
[strain_12]
order = CONSTANT
family = MONOMIAL
[]
[strain_21]
order = CONSTANT
family = MONOMIAL
[]
[]
[Kernels]
[solid_disp_x]
type = ADStressDivergenceShell
block = 1
component = 0
variable = disp_x
through_thickness_order = SECOND
[]
[solid_disp_y]
type = ADStressDivergenceShell
block = 1
component = 1
variable = disp_y
through_thickness_order = SECOND
[]
[solid_disp_z]
type = ADStressDivergenceShell
block = 1
component = 2
variable = disp_z
through_thickness_order = SECOND
[]
[solid_rot_x]
type = ADStressDivergenceShell
block = 1
component = 3
variable = rot_x
through_thickness_order = SECOND
[]
[solid_rot_y]
type = ADStressDivergenceShell
block = 1
component = 4
variable = rot_y
through_thickness_order = SECOND
[]
[]
[AuxKernels]
[stress_00]
type = RankTwoAux
variable = stress_00
rank_two_tensor = global_stress_t_points_0
index_i = 0
index_j = 0
execute_on = TIMESTEP_END
[]
[strain_00]
type = RankTwoAux
variable = strain_00
rank_two_tensor = total_global_strain_t_points_0
index_i = 0
index_j = 0
execute_on = TIMESTEP_END
[]
[stress_11]
type = RankTwoAux
variable = stress_11
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 1
execute_on = TIMESTEP_END
[]
[strain_11]
type = RankTwoAux
variable = strain_11
rank_two_tensor = total_global_strain_t_points_0
index_i = 1
index_j = 1
execute_on = TIMESTEP_END
[]
[stress_22]
type = RankTwoAux
variable = stress_22
rank_two_tensor = global_stress_t_points_0
index_i = 2
index_j = 2
execute_on = TIMESTEP_END
[]
[strain_22]
type = RankTwoAux
variable = strain_22
rank_two_tensor = total_global_strain_t_points_0
index_i = 2
index_j = 2
execute_on = TIMESTEP_END
[]
[stress_01]
type = RankTwoAux
variable = stress_01
rank_two_tensor = global_stress_t_points_0
index_i = 0
index_j = 1
execute_on = TIMESTEP_END
[]
[strain_01]
type = RankTwoAux
variable = strain_01
rank_two_tensor = total_global_strain_t_points_0
index_i = 0
index_j = 1
execute_on = TIMESTEP_END
[]
[stress_10]
type = RankTwoAux
variable = stress_10
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 0
execute_on = TIMESTEP_END
[]
[strain_10]
type = RankTwoAux
variable = strain_10
rank_two_tensor = total_global_strain_t_points_0
index_i = 1
index_j = 0
execute_on = TIMESTEP_END
[]
[stress_02]
type = RankTwoAux
variable = stress_02
rank_two_tensor = global_stress_t_points_0
index_i = 0
index_j = 2
execute_on = TIMESTEP_END
[]
[strain_02]
type = RankTwoAux
variable = strain_02
rank_two_tensor = total_global_strain_t_points_0
index_i = 0
index_j = 2
execute_on = TIMESTEP_END
[]
[stress_20]
type = RankTwoAux
variable = stress_20
rank_two_tensor = global_stress_t_points_0
index_i = 2
index_j = 0
execute_on = TIMESTEP_END
[]
[strain_20]
type = RankTwoAux
variable = strain_20
rank_two_tensor = total_global_strain_t_points_0
index_i = 2
index_j = 0
execute_on = TIMESTEP_END
[]
[stress_12]
type = RankTwoAux
variable = stress_12
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 2
execute_on = TIMESTEP_END
[]
[strain_12]
type = RankTwoAux
variable = strain_12
rank_two_tensor = total_global_strain_t_points_0
index_i = 1
index_j = 2
execute_on = TIMESTEP_END
[]
[stress_21]
type = RankTwoAux
variable = stress_21
rank_two_tensor = global_stress_t_points_0
index_i = 2
index_j = 1
execute_on = TIMESTEP_END
[]
[strain_21]
type = RankTwoAux
variable = strain_21
rank_two_tensor = total_global_strain_t_points_0
index_i = 2
index_j = 1
execute_on = TIMESTEP_END
[]
[]
[BCs]
[BC_0]
type = ADDirichletBC
variable = disp_x
value = 0.0
boundary = '2' #left
[]
[BC_1]
type = ADDirichletBC
variable = disp_y
value = 0.0
boundary = 10 #left_side_mid
[]
[BC_2]
type = FunctionDirichletBC
variable = disp_x
boundary = '3'
function = displacement
[]
[]
[Functions]
[displacement]
type = PiecewiseLinear
x = '0.0 1.0'
y = '0.0 0.2'
[]
[]
[Materials]
[stress]
type = ADComputeShellStress
block = 1
through_thickness_order = SECOND
[]
[elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 210000
poissons_ratio = 0.3
block = 1
through_thickness_order = SECOND
[]
[strain]
type = ADComputeIncrementalShellStrain
block = 1
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
thickness = 0.1
through_thickness_order = SECOND
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
line_search = 'none'
nl_rel_tol = 1e-16
nl_abs_tol = 1e-16
dt = 1
dtmin = 1
end_time = 1
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/plate_bending2.i)
# Shell element verification test from Abaqus verification manual 1.3.13
# A 40 m x 20 m x 1 m plate that has E = 1000 Pa and Poisson's ratio = 0.3
# is subjected to the following boundary/loading conditions. A single shell
# element is used to model the plate.
# disp_z = 0 at vertices A (0, 0), B (40, 0) and D (20, 0).
# disp_x and disp_y are zero at all four vertices.
# F_z = -2.0 N at vertex C (40, 20).
# M_x = 20.0 Nm at vertices A and B (bottom boundary)
# M_x = -20.0 Nm at vertices C and D (top boundary)
# M_y = 10.0 Nm at vertices B and C (right boundary)
# M_y = -10.0 Nm at vertices A and D (left boundary)
# The disp_z at vertex C is -12.54 m using S4 elements in Abaqus.
# The solution obtained using Moose is -12.519 m with a relative error
# of 0.16%.
[Mesh]
[./gmg]
type = GeneratedMeshGenerator
dim = 2
nx = 1
ny = 1
xmin = 0.0
xmax = 40.0
ymin = 0.0
ymax = 20.0
[../]
[./c_node]
type = ExtraNodesetGenerator
input = gmg
new_boundary = 100
coord = '40.0 20.0'
[../]
[]
[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 = 'right top bottom left'
value = 0.0
[../]
[./simply_support_y]
type = DirichletBC
variable = disp_y
boundary = 'right top bottom left'
value = 0.0
[../]
[./simply_support_z]
type = DirichletBC
variable = disp_z
boundary = 'bottom left'
value = 0.0
[../]
[]
[NodalKernels]
[./force_C]
type = ConstantRate
variable = disp_z
boundary = 100
rate = -2.0
[../]
[./Mx_AB]
type = ConstantRate
variable = rot_x
boundary = bottom
rate = 20.0
[../]
[./Mx_CD]
type = ConstantRate
variable = rot_x
boundary = top
rate = -20.0
[../]
[./My_BC]
type = ConstantRate
variable = rot_y
boundary = right
rate = 10.0
[../]
[./My_AD]
type = ConstantRate
variable = rot_y
boundary = left
rate = -10.0
[../]
[]
[Preconditioning]
[./smp]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = NEWTON
line_search = 'none'
#nl_max_its = 2
nl_rel_tol = 1e-10
nl_abs_tol = 6e-6
dt = 1.0
dtmin = 1.0
end_time = 3
[]
[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 = 1e3
poissons_ratio = 0.3
block = 0
through_thickness_order = SECOND
[../]
[./strain]
type = ADComputeIncrementalShellStrain
block = '0'
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
thickness = 1.0
through_thickness_order = SECOND
[../]
[./stress]
type = ADComputeShellStress
block = 0
through_thickness_order = SECOND
[../]
[]
[Postprocessors]
[./disp_z2]
type = PointValue
point = '40.0 20.0 0.0'
variable = disp_z
[../]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/tank_shell_rotated.i)
# Test for Pressure on shell elements
# An inclined cylindrical tank (length:3m) with a wall thickness of t=0.03 m and a radius of 0.5m is subjected to an internal presure of p=40 MPa.
# The lower part of the cylinder is constrained in all directions
# Theorically, assuming a thin_walled cylinder t/r <0.1, the hoop stress is sigma_t=p*r/t
# Therefore, in-plane force in the circumference of the cylinder is F=sigma_t*t= p*r=0.5*40=20 MN (independent of material properties of the shell)
# Analytical solution for the radial displacement : u=p*r^2/(E*t)=0.00167 m
# We check the axial_force_1 at the upper part of the cylinder (far from the lower boundary to avoid boundary effects)
# The numerical modeling results in axial_force_1 =19.882 MPa (0.6% error) and radial displacement u=0.00165 (1.1% error)
[Mesh]
[gmg]
type = FileMeshGenerator
file = tank_shell_rotated.msh
[]
[]
[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
[]
[]
[AuxVariables]
[axial_force_1]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[axial_force_1]
type = ShellResultantsAux
variable = axial_force_1
stress_resultant = axial_force_1
thickness = 0.03
through_thickness_order = SECOND
execute_on = TIMESTEP_END
[]
[]
[BCs]
[simply_support_x]
type = DirichletBC
variable = disp_x
boundary = 'lower_circle'
value = 0.0
[]
[simply_support_y]
type = DirichletBC
variable = disp_y
boundary = 'lower_circle'
value = 0.0
[]
[simply_support_z]
type = DirichletBC
variable = disp_z
boundary = 'lower_circle'
value = 0.0
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options = '-snes_ksp_ew'
# best overall
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = ' lu mumps'
line_search = 'none'
nl_rel_tol = 1e-12
nl_abs_tol = 1e-13
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
[]
[load_x]
type = ADDistributedLoadShell
function = '40'
variable = disp_x
project_load_to_normal = true
displacements = 'disp_x disp_y disp_z'
[]
[load_y]
type = ADDistributedLoadShell
function = '40'
variable = disp_y
project_load_to_normal = true
displacements = 'disp_x disp_y disp_z'
[]
[load_z]
type = ADDistributedLoadShell
function = '40'
variable = disp_z
project_load_to_normal = true
displacements = 'disp_x disp_y disp_z'
[]
[]
[Materials]
[elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 2e5
poissons_ratio = 0.3
through_thickness_order = SECOND
[]
[strain]
type = ADComputeIncrementalShellStrain
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
thickness = 0.03
reference_first_local_direction = ' 1 0 1'
through_thickness_order = SECOND
[]
[stress]
type = ADComputeShellStress
through_thickness_order = SECOND
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD.i)
# Test that models bending of a cantilever beam using shell elements
# A cantilever beam of length 10 m (in Y direction) and cross-section
# 1 m x 0.1 m is modeled using 4 shell elements placed along the length
# (Figure 6a from Dvorkin and Bathe, 1984). All displacements and
# X rotations are fixed on the bottom boundary. E = 2100000 and v = 0.0.
# A load of 0.5 N (in the Z direction) is applied at each node on the top
# boundary resulting in a total load of 1 N.
# The analytical solution for displacement at tip using small strain/rotations # is PL^3/3EI + PL/AG = 1.90485714 m
# The FEM solution using 4 shell elements is 1.875095 m with a relative error
# of 1.5%.
# Similarly, the analytical solution for slope at tip is PL^2/2EI = 0.285714286
# The FEM solution is 0.2857143 and the relative error is 5e-6%.
# The stress_yy for the four elements at z = -0.57735 * (t/2) (first qp below mid-surface of shell) are:
# 3031.089 Pa, 2165.064 Pa, 1299.038 Pa and 433.0127 Pa.
# Note the above values are the average stresses in each element.
# Analytically, stress_yy decreases linearly from y = 0 to y = 10 m.
# The maximum value of stress_yy at y = 0 is Mz/I = PL * 0.57735*(t/2)/I = 3464.1 Pa
# Therefore, the analytical value of stress at z = -0.57735 * (t/2) at the mid-point
# of the four elements are:
# 3031.0875 Pa, 2165.0625 Pa, 1299.0375 Pa ,433.0125 Pa
# The relative error in stress_yy is in the order of 5e-5%.
# The stress_yz at z = -0.57735 * (t/2) at all four elements from the simulation is 10 Pa.
# The analytical solution for the shear stress is: V/2/I *((t^2)/4 - z^2), where the shear force (V)
# is 1 N at any y along the length of the beam. Therefore, the analytical shear stress at
# z = -0.57735 * (t/2) is 10 Pa at any location along the length of the beam.
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 4
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 10.0
[]
[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
[]
[]
[AuxVariables]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_yy]
type = RankTwoAux
variable = stress_yy
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 1
[]
[stress_yz]
type = RankTwoAux
variable = stress_yz
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 2
[]
[]
[BCs]
[fixy1]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[fixz1]
type = DirichletBC
variable = disp_z
boundary = 'bottom'
value = 0.0
[]
[fixr1]
type = DirichletBC
variable = rot_x
boundary = 'bottom'
value = 0.0
[]
[fixr2]
type = DirichletBC
variable = rot_y
boundary = 'bottom'
value = 0.0
[]
[fixx1]
type = DirichletBC
variable = disp_x
boundary = 'bottom'
value = 0.0
[]
[]
[NodalKernels]
[force_y2]
type = ConstantRate
variable = disp_z
boundary = 'top'
rate = 0.5
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
nl_max_its = 2
nl_rel_tol = 1e-10
nl_abs_tol = 5e-4
dt = 1
dtmin = 1
end_time = 1
[]
[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 = 2100000
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]
[disp_z_tip]
type = PointValue
point = '1.0 10.0 0.0'
variable = disp_z
[]
[rot_x_tip]
type = PointValue
point = '0.0 10.0 0.0'
variable = rot_x
[]
[stress_yy_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yy
[]
[stress_yy_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yy
[]
[stress_yy_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yy
[]
[stress_yy_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yy
[]
[stress_yz_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yz
[]
[stress_yz_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yz
[]
[stress_yz_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yz
[]
[stress_yz_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yz
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/large_strain_m_40_AD.i)
# Large strain/rotation test for shell elements
# A cantilever beam that is 40 m long (Y direction) with 1 m x 1 m
# cross-section is modeled using 5 shell elements placed along its
# length. The bottom boundary is fixed in all displacements and
# rotations. A load of 0.140625 N is applied at each node on the top
# boundary, resulting in a total load of 0.28125 N. E = 1800 Pa and
# v = 0.0.
# The reference solution for large deflection of this beam is based on
# K. E. Bisshopp and D.C. Drucker, Quaterly of Applied Mathematics,
# Vol 3, No. # 3, 1945.
# 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 = 25.2 m; relative error = 5 %
# disp_y = -9.42 m; relative error = 1.87 %
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 5
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 40.0
[]
[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]
[fixy1]
type = DirichletBC
variable = disp_y
boundary = bottom
value = 0.0
[]
[fixz1]
type = DirichletBC
variable = disp_z
boundary = bottom
value = 0.0
[]
[fixr1]
type = DirichletBC
variable = rot_x
boundary = bottom
value = 0.0
[]
[fixr2]
type = DirichletBC
variable = rot_y
boundary = bottom
value = 0.0
[]
[fixx1]
type = DirichletBC
variable = disp_x
boundary = bottom
value = 0.0
[]
[]
[NodalKernels]
[force_y2]
type = UserForcingFunctorNodalKernel
variable = disp_z
boundary = top
functor = force_y
[]
[]
[Functions]
[force_y]
type = PiecewiseLinear
x = '0.0 1.0 3.0'
y = '0.0 1.0 1.0'
scale_factor = 0.140625
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
petsc_options_iname = '-pc_type'
petsc_options_value = ' lu'
line_search = 'none'
nl_rel_tol = 1e-10
nl_abs_tol = 1e-14
dt = 0.1
dtmin = 0.1
end_time = 1.0
[]
[Kernels]
[solid_disp_x]
type = ADStressDivergenceShell
block = '0'
component = 0
variable = disp_x
through_thickness_order = SECOND
large_strain = true
[]
[solid_disp_y]
type = ADStressDivergenceShell
block = '0'
component = 1
variable = disp_y
through_thickness_order = SECOND
large_strain = true
[]
[solid_disp_z]
type = ADStressDivergenceShell
block = '0'
component = 2
variable = disp_z
through_thickness_order = SECOND
large_strain = true
[]
[solid_rot_x]
type = ADStressDivergenceShell
block = '0'
component = 3
variable = rot_x
through_thickness_order = SECOND
large_strain = true
[]
[solid_rot_y]
type = ADStressDivergenceShell
block = '0'
component = 4
variable = rot_y
through_thickness_order = SECOND
large_strain = true
[]
[]
[Materials]
[elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 1800
poissons_ratio = 0.0
block = 0
through_thickness_order = SECOND
[]
[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]
type = ADComputeShellStress
block = 0
through_thickness_order = SECOND
[]
[]
[Postprocessors]
[disp_z2]
type = PointValue
point = '1.0 40.0 0.0'
variable = disp_z
[]
[disp_y2]
type = PointValue
point = '1.0 40.0 0.0'
variable = disp_y
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/scordelis_lo_roof_shell.i)
# This model is a widely used benchmark model denoted the Scordelis-Lo roof.
# The maximum z-deformation is compared with the value given in "Proposed Standard Set of Problems to Test Finite Element Accuracy, Finite Elements in Analysis and Design, 1985".
# Based on the existing analytical Solutions, maximum deflection of the roof should be 0.3086
# The model results in a maximum deflection of 0.3090 (assuming a 15*15 structured mesh)
[Mesh]
[file]
type = FileMeshGenerator
file = scordelis_lo_roof_shell.msh
[]
[]
[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_y]
type = ADDirichletBC
variable = disp_x
boundary = 'back'
value = 0.0
[]
[simply_support_z]
type = ADDirichletBC
variable = disp_z
boundary = 'back'
value = 0.0
[]
[simply_support_x]
type = ADDirichletBC
variable = disp_y
boundary = 'front'
value = 0.0
[]
[simply_rot_x]
type = ADDirichletBC
variable = rot_x
boundary = 'front'
value = 0.0
[]
[]
[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-10
nl_abs_tol = 1e-10
dt = 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
[]
[self_weight]
type = ADDistributedLoadShell
function = '90'
variable = disp_z
displacements = 'disp_x disp_y disp_z'
[]
[]
[Materials]
[elasticity_tshell]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 4.32e8
poissons_ratio = 0.0
through_thickness_order = SECOND
[]
[strain_shell]
type = ADComputeIncrementalShellStrain
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
thickness = 0.25
through_thickness_order = SECOND
[]
[stress_shell]
type = ADComputeShellStress
through_thickness_order = SECOND
[]
[]
[Postprocessors]
[disp_z2]
type = PointValue
point = '-16.7 0 19.2'
variable = disp_z
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/clamped_plate_flat.i)
# Test for simply supported plate under uniform pressure
# One quarter of a 50 m x 50 m x 1m plate is modeled in this test.
# Pressure loading is applied on the top surface using nodal forces
# of magnitude -10 N on all nodes. This corresponds to a pressure (q) of
# -10.816 N/m^2.
# The FEM solution at (0,0), which is at the center of the full plate
# is -3.003319e-03 m (for a 5*5 mesh).
# The analytical solution for displacement at center of plate obtained
# using a thin plate assumption for a square plate is
# w = 16 q a^4/(D*pi^6) \sum_{m = 1,3,5, ..}^\inf \sum_{n = 1,3,5, ..}^\inf (-1)^{(m+n-2)/2}/(mn*(m^2+n^2)^2)
# The above solution is the Naviers series solution from the "Theory of plates
# and shells" by Timoshenko and Woinowsky-Krieger (1959).
# where a = 50 m, q = -10.816 N/m^2 and D = E/(12(1-v^2))
# The analytical series solution converges to 2.998535904e-03 m
# when the first 16 terms of the series are considered (i.e., until
# m & n = 7).
[Mesh]
[gmg]
type = FileMeshGenerator
file = clamped_plate_flat.msh
[]
[]
[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
[]
[symm_bottom_rot]
type = DirichletBC
variable = rot_x
boundary = 'bottom'
value = 0.0
[]
[simply_support_x]
type = DirichletBC
variable = disp_x
boundary = 'right top bottom left'
value = 0.0
[]
[simply_support_y]
type = DirichletBC
variable = disp_y
boundary = 'right top bottom left'
value = 0.0
[]
[simply_support_z]
type = DirichletBC
variable = disp_z
boundary = 'right top'
value = 0.0
[]
[]
[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-10
nl_abs_tol = 1e-10
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
[]
[load_z]
type = ADDistributedLoadShell
function = '10.816'
variable = disp_z
displacements = 'disp_x disp_y disp_z'
[]
[]
[Materials]
[elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 1e9
poissons_ratio = 0.3
through_thickness_order = SECOND
block = 'shell'
[]
[strain]
type = ADComputeIncrementalShellStrain
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
thickness = 1
through_thickness_order = SECOND
block = 'shell'
[]
[stress]
type = ADComputeShellStress
through_thickness_order = SECOND
block = 'shell'
[]
[]
[Postprocessors]
[disp_z2]
type = PointValue
point = '0.0 0.0 0.0'
variable = disp_z
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/plate_cantilever.i)
#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
[]
(modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment.i)
# Test that models bending of a cantilever beam using shell elements
# A cantilever beam of length 10 m (in Y direction) and cross-section
# 1 m x 0.1 m is modeled using 4 shell elements placed along the length
# (Figure 6a from Dvorkin and Bathe, 1984). All displacements and
# X rotations are fixed on the bottom boundary. E = 2100000 and v = 0.0.
# A load of 0.5 N (in the Z direction) is applied at each node on the top
# boundary resulting in a total load of 1 N.
# The analytical solution for displacement at tip using small strain/rotations # is PL^3/3EI + PL/AG = 1.90485714 m
# The FEM solution using 4 shell elements is 1.875095 m with a relative error
# of 1.5%.
# Similarly, the analytical solution for slope at tip is PL^2/2EI = 0.285714286
# The FEM solution is 0.2857143 and the relative error is 5e-6%.
# The stress_yy for the four elements at z = -0.57735 * (t/2) (first qp below mid-surface of shell) are:
# 3031.089 Pa, 2165.064 Pa, 1299.038 Pa and 433.0127 Pa.
# Note the above values are the average stresses in each element.
# Analytically, stress_yy decreases linearly from y = 0 to y = 10 m.
# The maximum value of stress_yy at y = 0 is Mz/I = PL * 0.57735*(t/2)/I = 3464.1 Pa
# Therefore, the analytical value of stress at z = -0.57735 * (t/2) at the mid-point
# of the four elements are:
# 3031.0875 Pa, 2165.0625 Pa, 1299.0375 Pa ,433.0125 Pa
# The relative error in stress_yy is in the order of 5e-5%.
# The stress_yz at z = -0.57735 * (t/2) at all four elements from the simulation is 10 Pa.
# The analytical solution for the shear stress is: V/2/I *((t^2)/4 - z^2), where the shear force (V)
# is 1 N at any y along the length of the beam. Therefore, the analytical shear stress at
# z = -0.57735 * (t/2) is 10 Pa at any location along the length of the beam.
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 4
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 10.0
[]
[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
[../]
[]
[AuxVariables]
[./stress_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz]
order = CONSTANT
family = MONOMIAL
[../]
# aux variables for dynamics
[./vel_x]
order = FIRST
family = LAGRANGE
[../]
[./vel_y]
order = FIRST
family = LAGRANGE
[../]
[./vel_z]
order = FIRST
family = LAGRANGE
[../]
[./accel_x]
order = FIRST
family = LAGRANGE
[../]
[./accel_y]
order = FIRST
family = LAGRANGE
[../]
[./accel_z]
order = FIRST
family = LAGRANGE
[../]
[./rot_vel_x]
order = FIRST
family = LAGRANGE
[../]
[./rot_vel_y]
order = FIRST
family = LAGRANGE
[../]
[./rot_accel_x]
order = FIRST
family = LAGRANGE
[../]
[./rot_accel_y]
order = FIRST
family = LAGRANGE
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
variable = stress_yy
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 1
[../]
[./stress_yz]
type = RankTwoAux
variable = stress_yz
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 2
[../]
# Kernels for dynamics
[./accel_x]
type = NewmarkAccelAux
variable = accel_x
displacement = disp_x
velocity = vel_x
beta = 0.25
execute_on = timestep_end
[../]
[./vel_x]
type = NewmarkVelAux
variable = vel_x
acceleration = accel_x
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_y]
type = NewmarkAccelAux
variable = accel_y
displacement = disp_y
velocity = vel_y
beta = 0.25
execute_on = timestep_end
[../]
[./vel_y]
type = NewmarkVelAux
variable = vel_y
acceleration = accel_y
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_z]
type = NewmarkAccelAux
variable = accel_z
displacement = disp_z
velocity = vel_z
beta = 0.25
execute_on = timestep_end
[../]
[./vel_z]
type = NewmarkVelAux
variable = vel_z
acceleration = accel_z
gamma = 0.5
execute_on = timestep_end
[../]
[./rot_accel_x]
type = NewmarkAccelAux
variable = rot_accel_x
displacement = rot_x
velocity = rot_vel_x
beta = 0.25
execute_on = timestep_end
[../]
[./rot_vel_x]
type = NewmarkVelAux
variable = rot_vel_x
acceleration = rot_accel_x
gamma = 0.5
execute_on = timestep_end
[../]
[./rot_accel_y]
type = NewmarkAccelAux
variable = rot_accel_y
displacement = rot_y
velocity = rot_vel_y
beta = 0.25
execute_on = timestep_end
[../]
[./rot_vel_y]
type = NewmarkVelAux
variable = rot_vel_y
acceleration = rot_accel_y
gamma = 0.5
execute_on = timestep_end
[../]
[]
[BCs]
[./fixy1]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[../]
[./fixz1]
type = DirichletBC
variable = disp_z
boundary = 'bottom'
value = 0.0
[../]
[./fixr1]
type = DirichletBC
variable = rot_x
boundary = 'bottom'
value = 0.0
[../]
[./fixr2]
type = DirichletBC
variable = rot_y
boundary = 'bottom'
value = 0.0
[../]
[./fixx1]
type = DirichletBC
variable = disp_x
boundary = 'bottom'
value = 0.0
[../]
[]
[Functions]
[./force_function]
type = PiecewiseLinear
x = '0.0 1.0'
y = '0.0 0.5'
[../]
[]
[NodalKernels]
[./force_y2]
type = UserForcingFunctorNodalKernel
variable = disp_z
boundary = 'top'
functor = force_function
[../]
[]
[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
[../]
[./inertial_force_x]
type = ADInertialForceShell
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 0
variable = disp_x
thickness = 0.1
[../]
[./inertial_force_y]
type = ADInertialForceShell
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 1
variable = disp_y
thickness = 0.1
[../]
[./inertial_force_z]
type = ADInertialForceShell
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 2
variable = disp_z
thickness = 0.1
[../]
[./inertial_force_rot_x]
type = ADInertialForceShell
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 3
variable = rot_x
thickness = 0.1
[../]
[./inertial_force_rot_y]
type = ADInertialForceShell
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 4
variable = rot_y
thickness = 0.1
[../]
[]
[Materials]
[./elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 2100000
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
[../]
[./density]
type = GenericConstantMaterial
block = 0
prop_names = 'density'
prop_values = '1.0'
[../]
[]
[Postprocessors]
[./disp_z_tip]
type = PointValue
point = '1.0 10.0 0.0'
variable = disp_z
[../]
[./rot_x_tip]
type = PointValue
point = '0.0 10.0 0.0'
variable = rot_x
[../]
[./stress_yy_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yy
[../]
[./stress_yy_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yy
[../]
[./stress_yy_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yy
[../]
[./stress_yy_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yy
[../]
[./stress_yz_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yz
[../]
[./stress_yz_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yz
[../]
[./stress_yz_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yz
[../]
[./stress_yz_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yz
[../]
[]
[Preconditioning]
[./smp]
type = SMP
full = true
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
petsc_options_iname = '-pc_type'
petsc_options_value = 'lu'
nl_max_its = 2
nl_rel_tol = 1e-10
nl_abs_tol = 5e-8
dt = 0.0005
dtmin = 0.0005
end_time = 1
[TimeIntegrator]
type = NewmarkBeta
beta = 0.25
gamma = 0.5
[]
[]
[Outputs]
csv = true
[]
(modules/solid_mechanics/test/tests/shell/static/pressure_error.i)
# Test for the stress and strain output for tapered shell elements.
# A tapered beam is represented with shell elements in XY plane
# having Young's Modulus of 210000 and poissons ratio of 0.3.
# The displacement in X direction is constrained in the left end and the
# displacement of center node of the left end is constrained in Y direction.
# A uniform pressure of 50 units is applied at the right end.
# The shell element strain calculation assumes four quadrature points.
# So, when a pressure is applied at the edge, two quadrature points
# are used in the strain calulator and an error is generated.
[Mesh]
[input]
type = FileMeshGenerator
file = taperedmesh.e
[]
[]
[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_10]
order = CONSTANT
family = MONOMIAL
[]
[stress_02]
order = CONSTANT
family = MONOMIAL
[]
[stress_20]
order = CONSTANT
family = MONOMIAL
[]
[stress_12]
order = CONSTANT
family = MONOMIAL
[]
[stress_21]
order = CONSTANT
family = MONOMIAL
[]
[strain_00]
order = CONSTANT
family = MONOMIAL
[]
[strain_11]
order = CONSTANT
family = MONOMIAL
[]
[strain_22]
order = CONSTANT
family = MONOMIAL
[]
[strain_01]
order = CONSTANT
family = MONOMIAL
[]
[strain_10]
order = CONSTANT
family = MONOMIAL
[]
[strain_02]
order = CONSTANT
family = MONOMIAL
[]
[strain_20]
order = CONSTANT
family = MONOMIAL
[]
[strain_12]
order = CONSTANT
family = MONOMIAL
[]
[strain_21]
order = CONSTANT
family = MONOMIAL
[]
[]
[Kernels]
[solid_disp_x]
type = ADStressDivergenceShell
block = 1
component = 0
variable = disp_x
through_thickness_order = SECOND
[]
[solid_disp_y]
type = ADStressDivergenceShell
block = 1
component = 1
variable = disp_y
through_thickness_order = SECOND
[]
[solid_disp_z]
type = ADStressDivergenceShell
block = 1
component = 2
variable = disp_z
through_thickness_order = SECOND
[]
[solid_rot_x]
type = ADStressDivergenceShell
block = 1
component = 3
variable = rot_x
through_thickness_order = SECOND
[]
[solid_rot_y]
type = ADStressDivergenceShell
block = 1
component = 4
variable = rot_y
through_thickness_order = SECOND
[]
[]
[AuxKernels]
[stress_00]
type = RankTwoAux
variable = stress_00
rank_two_tensor = global_stress_t_points_0
index_i = 0
index_j = 0
execute_on = TIMESTEP_END
[]
[strain_00]
type = RankTwoAux
variable = strain_00
rank_two_tensor = total_global_strain_t_points_0
index_i = 0
index_j = 0
execute_on = TIMESTEP_END
[]
[stress_11]
type = RankTwoAux
variable = stress_11
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 1
execute_on = TIMESTEP_END
[]
[strain_11]
type = RankTwoAux
variable = strain_11
rank_two_tensor = total_global_strain_t_points_0
index_i = 1
index_j = 1
execute_on = TIMESTEP_END
[]
[stress_22]
type = RankTwoAux
variable = stress_22
rank_two_tensor = global_stress_t_points_0
index_i = 2
index_j = 2
execute_on = TIMESTEP_END
[]
[strain_22]
type = RankTwoAux
variable = strain_22
rank_two_tensor = total_global_strain_t_points_0
index_i = 2
index_j = 2
execute_on = TIMESTEP_END
[]
[stress_01]
type = RankTwoAux
variable = stress_01
rank_two_tensor = global_stress_t_points_0
index_i = 0
index_j = 1
execute_on = TIMESTEP_END
[]
[strain_01]
type = RankTwoAux
variable = strain_01
rank_two_tensor = total_global_strain_t_points_0
index_i = 0
index_j = 1
execute_on = TIMESTEP_END
[]
[stress_10]
type = RankTwoAux
variable = stress_10
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 0
execute_on = TIMESTEP_END
[]
[strain_10]
type = RankTwoAux
variable = strain_10
rank_two_tensor = total_global_strain_t_points_0
index_i = 1
index_j = 0
execute_on = TIMESTEP_END
[]
[stress_02]
type = RankTwoAux
variable = stress_02
rank_two_tensor = global_stress_t_points_0
index_i = 0
index_j = 2
execute_on = TIMESTEP_END
[]
[strain_02]
type = RankTwoAux
variable = strain_02
rank_two_tensor = total_global_strain_t_points_0
index_i = 0
index_j = 2
execute_on = TIMESTEP_END
[]
[stress_20]
type = RankTwoAux
variable = stress_20
rank_two_tensor = global_stress_t_points_0
index_i = 2
index_j = 0
execute_on = TIMESTEP_END
[]
[strain_20]
type = RankTwoAux
variable = strain_20
rank_two_tensor = total_global_strain_t_points_0
index_i = 2
index_j = 0
execute_on = TIMESTEP_END
[]
[stress_12]
type = RankTwoAux
variable = stress_12
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 2
execute_on = TIMESTEP_END
[]
[strain_12]
type = RankTwoAux
variable = strain_12
rank_two_tensor = total_global_strain_t_points_0
index_i = 1
index_j = 2
execute_on = TIMESTEP_END
[]
[stress_21]
type = RankTwoAux
variable = stress_21
rank_two_tensor = global_stress_t_points_0
index_i = 2
index_j = 1
execute_on = TIMESTEP_END
[]
[strain_21]
type = RankTwoAux
variable = strain_21
rank_two_tensor = total_global_strain_t_points_0
index_i = 2
index_j = 1
execute_on = TIMESTEP_END
[]
[]
[BCs]
[BC_0]
type = ADDirichletBC
variable = disp_x
value = 0.0
boundary = '2' #left
[]
[BC_1]
type = ADDirichletBC
variable = disp_y
value = 0.0
boundary = 10 #left_side_mid
[]
[BC_2]
type = ADPressure
variable = disp_x
displacements = 'disp_x disp_y disp_z'
factor = -50
boundary = '3'
use_displaced_mesh = false
[]
[]
[Materials]
[stress]
type = ADComputeShellStress
block = 1
through_thickness_order = SECOND
[]
[elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 210000
poissons_ratio = 0.3
block = 1
through_thickness_order = SECOND
[]
[strain]
type = ADComputeIncrementalShellStrain
block = 1
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
thickness = 0.1
through_thickness_order = SECOND
[]
[]
[Executioner]
type = Steady
solve_type = NEWTON
nl_abs_tol = 1e-8
l_abs_tol = 1e-8
nl_max_its = 100
l_max_its = 100
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/inclined_straintest_local_stress.i)
# 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
[]
(modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm.i)
# Test for displacement of pinched cylinder
# Ref: Figure 10 and Table 6 from Dvorkin and Bathe, Eng. Comput., Vol. 1, 1984.
# A cylinder of radius 1 m and length 2 m (along Z axis) with clamped ends
# (at z = 0 and 2 m) is pinched at mid-length by placing point loads of 10 N
# at (1, 0, 1) and (-1, 0, 1). Due to the symmetry of the problem, only 1/8th
# of the cylinder needs to be modeled.
# The normalized series solution for the displacement at the loading point is
# w = Wc E t / P = 164.24; where Wc is the displacement in m, E is the Young's
# modulus, t is the thickness and P is the point load.
# For this problem, E = 1e6 Pa, L = 2 m, R = 1 m, t = 0.01 m, P = 10 N and
# Poisson's ratio = 0.3. This gives an analytic displacement of 0.16424 m.
# FEM results from different mesh discretizations are presented below. Only
# the 10x10 mesh is included as a test.
# As shown in the table below, the results from the MOOSE FEM analysis converge
# to the analytic solution and the convergence matches well with the results
# of Dvorkin and Bathe (1984).
# Mesh of 1/8 cylinder | FEM/analytical disp | FEM/analytical disp
# | (MOOSE implementation) | (Reported by Dvorkin)
#----------------------|-------------------------|-------------------------
# 10 x 10 | 0.82 | 0.83
# 20 x 20 | 0.95 | 0.96
# 40 x 40 | 0.99 | -
# 80 x 80 | 1.01 | -
[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
[]
# Note that the rotational DOFs are in the local coordinate system
# Also it isn't clear from the Dvorkin paper which DOFs should be fixed on the far
# end (boundary CD). If it were fully constrained we would need to fix disp_z and
# the rotations, but that makes it stiffer than the analytical solution.
[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
[]
[]
[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
[]
[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
[]
(modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment_free_orientation_inclined.i)
# Test to verify the fundamental natural frequency of a one element ADComputeShellStress
# BCs: Clamped on one end, free on others.
# Initial perturbation applied to edge of the beam. After that, the shell vibrates freely.
#
# Results have been compared for various thicknesses with the following approximate Results
# (Moose results were obtained with 8 elements along the length)
# Thickness = 0.1. Reference freq: 10.785 Hz, Moose freq: 10.612 Hz
# Thickness = 0.05. Reference freq: 5.393 Hz, Moose freq: 5.335 Hz
# Thickness = 0.025. Reference freq: 2.696 Hz, Moose freq: 2.660 Hz
#
# Reference values have been obtained from Robert Blevins, "Formulas for Dynamics, Acoustics and Vibration",
# Table 5.3 case 11. Formula looks like: f = lambda^2/(2*pi*a^2) * sqrt(E*h^2/(12*(1-nu*nu))), where lambda
# changes as a function of shell dimensions.
# This test uses one single element for speed reasons.
# Here, the shell, instead of being on the XY plane, is oriented at a 45 deg. angle
# with respect to the Y axis.
[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_yy]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz]
order = CONSTANT
family = MONOMIAL
[../]
# aux variables for dynamics
[./vel_x]
[../]
[./vel_y]
[../]
[./vel_z]
[../]
[./accel_x]
[../]
[./accel_y]
[../]
[./accel_z]
[../]
[./rot_vel_x]
[../]
[./rot_vel_y]
[../]
[./rot_accel_x]
[../]
[./rot_accel_y]
[../]
[]
[AuxKernels]
[./stress_yy]
type = RankTwoAux
variable = stress_yy
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 1
[../]
[./stress_yz]
type = RankTwoAux
variable = stress_yz
rank_two_tensor = global_stress_t_points_0
index_i = 1
index_j = 2
[../]
# Kernels for dynamics
[./accel_x]
type = NewmarkAccelAux
variable = accel_x
displacement = disp_x
velocity = vel_x
beta = 0.25
execute_on = timestep_end
[../]
[./vel_x]
type = NewmarkVelAux
variable = vel_x
acceleration = accel_x
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_y]
type = NewmarkAccelAux
variable = accel_y
displacement = disp_y
velocity = vel_y
beta = 0.25
execute_on = timestep_end
[../]
[./vel_y]
type = NewmarkVelAux
variable = vel_y
acceleration = accel_y
gamma = 0.5
execute_on = timestep_end
[../]
[./accel_z]
type = NewmarkAccelAux
variable = accel_z
displacement = disp_z
velocity = vel_z
beta = 0.25
execute_on = timestep_end
[../]
[./vel_z]
type = NewmarkVelAux
variable = vel_z
acceleration = accel_z
gamma = 0.5
execute_on = timestep_end
[../]
[./rot_accel_x]
type = NewmarkAccelAux
variable = rot_accel_x
displacement = rot_x
velocity = rot_vel_x
beta = 0.25
execute_on = timestep_end
[../]
[./rot_vel_x]
type = NewmarkVelAux
variable = rot_vel_x
acceleration = rot_accel_x
gamma = 0.5
execute_on = timestep_end
[../]
[./rot_accel_y]
type = NewmarkAccelAux
variable = rot_accel_y
displacement = rot_y
velocity = rot_vel_y
beta = 0.25
execute_on = timestep_end
[../]
[./rot_vel_y]
type = NewmarkVelAux
variable = rot_vel_y
acceleration = rot_accel_y
gamma = 0.5
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
[../]
[]
[Functions]
[./force_function]
type = PiecewiseLinear
x = '0.0 0.01 0.15 10.0'
y = '0.0 0.01 0.0 0.0'
[../]
[]
[NodalKernels]
[./force_y2]
type = UserForcingFunctorNodalKernel
variable = disp_z
boundary = '2'
functor = force_function
[../]
[]
[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
[../]
[./inertial_force_x]
type = ADInertialForceShell
use_displaced_mesh = true
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 0
variable = disp_x
thickness = 0.1
[../]
[./inertial_force_y]
type = ADInertialForceShell
use_displaced_mesh = true
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 1
variable = disp_y
thickness = 0.1
[../]
[./inertial_force_z]
type = ADInertialForceShell
use_displaced_mesh = true
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 2
variable = disp_z
thickness = 0.1
[../]
[./inertial_force_rot_x]
type = ADInertialForceShell
use_displaced_mesh = true
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 3
variable = rot_x
thickness = 0.1
[../]
[./inertial_force_rot_y]
type = ADInertialForceShell
use_displaced_mesh = true
eta = 0.0
block = 0
displacements = 'disp_x disp_y disp_z'
rotations = 'rot_x rot_y'
velocities = 'vel_x vel_y vel_z'
accelerations = 'accel_x accel_y accel_z'
rotational_velocities = 'rot_vel_x rot_vel_y'
rotational_accelerations = 'rot_accel_x rot_accel_y'
component = 4
variable = rot_y
thickness = 0.1
[../]
[]
[Materials]
[./elasticity]
type = ADComputeIsotropicElasticityTensorShell
youngs_modulus = 2100000
poissons_ratio = 0.3
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
[../]
[./density]
type = GenericConstantMaterial
block = 0
prop_names = 'density'
prop_values = '1.0'
[../]
[]
[Postprocessors]
[./disp_z_tip]
type = PointValue
point = '0.0 1.06 1.06'
variable = disp_z
[../]
[./rot_x_tip]
type = PointValue
point = '0.0 1.06 1.06'
variable = rot_x
[../]
[./stress_yy_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yy
[../]
[./stress_yy_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yy
[../]
[./stress_yy_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yy
[../]
[./stress_yy_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yy
[../]
[./stress_yz_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_yz
[../]
[./stress_yz_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_yz
[../]
[./stress_yz_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_yz
[../]
[./stress_yz_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_yz
[../]
[]
[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 = 0.005
dtmin = 0.005
timestep_tolerance = 2e-13
end_time = 0.5
[./TimeIntegrator]
type = NewmarkBeta
beta = 0.25
gamma = 0.5
[../]
[]
[Outputs]
perf_graph = true
csv = true
[]
(modules/solid_mechanics/test/tests/shell/static/straintest_shear.i)
# Test for the shear stress and strain output for 2D planar shell with uniform mesh.
# A cantilever beam of length 10 m and cross-section 1.5 m x 0.1 m having
# Young's Modulus of 5 N/mm^2 and poissons ratio of 0 is subjected to shear
# displacement of 0.05 m at the free end.
[Mesh]
type = GeneratedMesh
dim = 2
nx = 4
ny = 1
xmin = 0.0
xmax = 10
ymin = 0.0
ymax = 1.5
[]
[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
[]
[]
[AuxVariables]
[stress_xx]
order = CONSTANT
family = MONOMIAL
[]
[strain_xx]
order = CONSTANT
family = MONOMIAL
[]
[stress_yy]
order = CONSTANT
family = MONOMIAL
[]
[strain_yy]
order = CONSTANT
family = MONOMIAL
[]
[stress_xy]
order = CONSTANT
family = MONOMIAL
[]
[strain_xy]
order = CONSTANT
family = MONOMIAL
[]
[stress_yz]
order = CONSTANT
family = MONOMIAL
[]
[strain_yz]
order = CONSTANT
family = MONOMIAL
[]
[stress_xz]
order = CONSTANT
family = MONOMIAL
[]
[strain_xz]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_xx]
type = RankTwoAux
variable = stress_xx
selected_qp = 0
rank_two_tensor = global_stress_t_points_1
index_i = 0
index_j = 0
[]
[strain_xx]
type = RankTwoAux
variable = strain_xx
rank_two_tensor = total_global_strain_t_points_1
selected_qp = 0
index_i = 0
index_j = 0
[]
[stress_yy]
type = RankTwoAux
variable = stress_yy
rank_two_tensor = global_stress_t_points_1
selected_qp = 0
index_i = 1
index_j = 1
[]
[strain_yy]
type = RankTwoAux
variable = strain_yy
rank_two_tensor = total_global_strain_t_points_1
selected_qp = 0
index_i = 1
index_j = 1
[]
[stress_xy]
type = RankTwoAux
variable = stress_xy
rank_two_tensor = global_stress_t_points_1
selected_qp = 0
index_i = 0
index_j = 1
[]
[strain_xy]
type = RankTwoAux
variable = strain_xy
rank_two_tensor = total_global_strain_t_points_1
selected_qp = 0
index_i = 0
index_j = 1
[]
[stress_yz]
type = RankTwoAux
variable = stress_yz
rank_two_tensor = global_stress_t_points_1
selected_qp = 0
index_i = 1
index_j = 2
[]
[strain_yz]
type = RankTwoAux
variable = strain_yz
rank_two_tensor = total_global_strain_t_points_1
selected_qp = 0
index_i = 1
index_j = 2
[]
[stress_xz]
type = RankTwoAux
variable = stress_xz
rank_two_tensor = global_stress_t_points_1
selected_qp = 0
index_i = 0
index_j = 2
[]
[strain_xz]
type = RankTwoAux
variable = strain_yz
rank_two_tensor = total_global_strain_t_points_1
selected_qp = 0
index_i = 0
index_j = 2
[]
[]
[BCs]
[fixx]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[]
[fixy]
type = DirichletBC
variable = disp_y
boundary = left
value = 0.0
[]
[fixz]
type = DirichletBC
variable = disp_z
boundary = left
value = 0.0
[]
[fixr1]
type = DirichletBC
variable = rot_x
boundary = left
value = 0.0
[]
[fixr2]
type = DirichletBC
variable = rot_y
boundary = left
value = 0.0
[]
[disp]
type = FunctionDirichletBC
variable = disp_y
boundary = 'right'
function = displacement
[]
[]
[Functions]
[displacement]
type = PiecewiseLinear
x = '0.0 1.0'
y = '0.0 0.05'
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
line_search = 'none'
nl_rel_tol = 1e-12
nl_abs_tol = 1e-14
dt = 1
dtmin = 1
end_time = 1
[]
[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 = 4.0e6
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_xy_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_xy
[]
[strain_xy_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_xy
[]
[stress_xy_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_xy
[]
[strain_xy_el_1]
type = ElementalVariableValue
elementid = 1
variable = strain_xy
[]
[stress_xy_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_xy
[]
[strain_xy_el_2]
type = ElementalVariableValue
elementid = 2
variable = strain_xy
[]
[stress_xy_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_xy
[]
[strain_xy_el_3]
type = ElementalVariableValue
elementid = 3
variable = strain_xy
[]
[stress_xx_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_xx
[]
[strain_xx_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_xx
[]
[stress_xx_el_1]
type = ElementalVariableValue
elementid = 1
variable = stress_xx
[]
[strain_xx_el_1]
type = ElementalVariableValue
elementid = 1
variable = strain_xx
[]
[stress_xx_el_2]
type = ElementalVariableValue
elementid = 2
variable = stress_xx
[]
[strain_xx_el_2]
type = ElementalVariableValue
elementid = 2
variable = strain_xx
[]
[stress_xx_el_3]
type = ElementalVariableValue
elementid = 3
variable = stress_xx
[]
[strain_xx_el_3]
type = ElementalVariableValue
elementid = 3
variable = strain_xx
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_unstructured.i)
# Test for displacement of pinched cylinder (similar to pinch_cyl_symm.i)
# This variant of the test is run with an unstructured mesh
[Mesh]
[mesh]
type = FileMeshGenerator
file = pinched_cyl_10_10_unstructured.msh
[]
[block_100]
type = ParsedSubdomainMeshGenerator
input = mesh
combinatorial_geometry = 'x > -1.1 & x < 1.1 & y > -1.1 & y < 1.1 & z > -0.1 & z < 2.1'
block_id = 100
[]
[nodeset_1]
type = BoundingBoxNodeSetGenerator
input = block_100
top_right = '1.1 1.1 0'
bottom_left = '-1.1 -1.1 0'
new_boundary = 'CD' #CD
[]
[nodeset_2]
type = BoundingBoxNodeSetGenerator
input = nodeset_1
top_right = '1.1 1.1 1.0'
bottom_left = '-1.1 -1.1 1.0'
new_boundary = 'AB' #AB
[]
[nodeset_3]
type = BoundingBoxNodeSetGenerator
input = nodeset_2
top_right = '0.02 1.1 1.0'
bottom_left = '-0.1 0.98 0.0'
new_boundary = 'AD' #AD
[]
[nodeset_4]
type = BoundingBoxNodeSetGenerator
input = nodeset_3
top_right = '1.1 0.02 1.0'
bottom_left = '0.98 -0.1 0.0'
new_boundary = 'BC' #BC
[]
[]
[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
[]
[]
[AuxVariables]
[stress_xx0]
order = CONSTANT
family = MONOMIAL
[]
[stress_xx1]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_xx0]
type = RankTwoAux
variable = stress_xx0
rank_two_tensor = global_stress_t_points_0
index_i = 0
index_j = 0
execute_on = TIMESTEP_END
[]
[stress_xx1]
type = RankTwoAux
variable = stress_xx1
rank_two_tensor = global_stress_t_points_1
index_i = 0
index_j = 0
execute_on = TIMESTEP_END
[]
[]
[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 = 'AD BC'
value = 0.0
[]
[simply_support_rot_y]
type = DirichletBC
variable = rot_y
boundary = 'AB'
value = 0.0
[]
[]
[DiracKernels]
[point1]
type = ConstantPointSource
variable = disp_x
point = '1 0 1'
value = -2.5 # P = 10
[]
[]
[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-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 = '0 0 1'
[]
[stress]
type = ADComputeShellStress
block = '100'
through_thickness_order = SECOND
[]
[]
[Postprocessors]
[disp_z2]
type = PointValue
point = '1 0 1'
variable = disp_x
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/straintest.i)
# Test for the axial stress and strain output for single shell element
# for 2D planar shell with uniform mesh.
# A single shell 1 mm x 1 mm element having Young's Modulus of 5 N/mm^2
# and poissons ratio of 0 is fixed at the left end and
# an axial displacement of 0.2 mm is applied at the right.
# Theoretical value of axial stress and strain are 1 N/mm^2 and 0.2.
[Mesh]
type = GeneratedMesh
dim = 2
nx = 1
ny = 1
xmin = 0.0
xmax = 1.0
ymin = 0.0
ymax = 1.0
[]
[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
[]
[]
[AuxVariables]
[stress_xx]
order = CONSTANT
family = MONOMIAL
[]
[strain_xx]
order = CONSTANT
family = MONOMIAL
[]
[]
[AuxKernels]
[stress_xx]
type = RankTwoAux
variable = stress_xx
rank_two_tensor = global_stress_t_points_1
index_i = 0
index_j = 0
[]
[strain_xx]
type = RankTwoAux
variable = strain_xx
rank_two_tensor = total_global_strain_t_points_1
index_i = 0
index_j = 0
[]
[]
[BCs]
[fixx]
type = DirichletBC
variable = disp_x
boundary = left
value = 0.0
[]
[fixy]
type = DirichletBC
variable = disp_y
boundary = left
value = 0.0
[]
[fixz]
type = DirichletBC
variable = disp_z
boundary = left
value = 0.0
[]
[fixr1]
type = DirichletBC
variable = rot_x
boundary = left
value = 0.0
[]
[fixr2]
type = DirichletBC
variable = rot_y
boundary = left
value = 0.0
[]
[disp]
type = FunctionDirichletBC
variable = disp_x
boundary = 'right'
function = displacement
[]
[]
[Functions]
[displacement]
type = PiecewiseLinear
x = '0.0 1.0'
y = '0.0 0.2'
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
automatic_scaling = true
line_search = 'none'
nl_rel_tol = 1e-12
nl_abs_tol = 1e-14
dt = 1
dtmin = 1
end_time = 1
[]
[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.0
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]
[disp_x]
type = PointValue
point = '0.5 0.0 0.0'
variable = disp_z
[]
[stress_xx_el_0]
type = ElementalVariableValue
elementid = 0
variable = stress_xx
[]
[strain_xx_el_0]
type = ElementalVariableValue
elementid = 0
variable = strain_xx
[]
[]
[Outputs]
exodus = true
[]
(modules/solid_mechanics/test/tests/shell/static/plate_concentrated_loads.i)
# 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<x<3 =10000
# Numerical model: out of plane shear force (q13) for 0<x<3 =10000
# Analytical solution: out of plane shear force (q13) for 3<x<6 =0
# Numerical model: out of plane shear force (q13) at for 3<x<6 =0
# Analytical solution: out of plane shear force (q13) for 6<x<9 =-10000
# Numerical model: out of plane shear force (q13) for 6<x<9 =-10000
[Mesh]
[gmg]
type = FileMeshGenerator
file = Plate_Concentrated_Loads.msh
[]
[p1]
type = BoundingBoxNodeSetGenerator
input = gmg
bottom_left = '2.99 0.0 -0.1'
top_right = '3.1 0.0 1.1'
new_boundary = 100
[]
[p2]
type = BoundingBoxNodeSetGenerator
input = p1
bottom_left = '5.99 0.0 -0.1'
top_right = '6.1 0.0 1.1'
new_boundary = 101
[]
[]
[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 = 'right left'
value = 0.0
[]
[simply_support_y]
type = DirichletBC
variable = disp_y
boundary = 'right left'
value = 0.0
[]
[]
[NodalKernels]
[force_p1]
type = ConstantRate
variable = disp_y
boundary = 100
rate = -2500 # applied to the four nodes at x=3
[]
[force_p2]
type = ConstantRate
variable = disp_y
boundary = 101
rate = -2500 # applied to the four nodes at x=6
[]
[]
[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-10
nl_abs_tol = 1e-10
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 = 200e9
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
[]
[]
[AuxVariables]
[moment_22]
order = CONSTANT
family = MONOMIAL
[]
[shear_13]
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
[]
[shear_13]
type = ShellResultantsAux
variable = shear_13
stress_resultant = shear_force_02
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
[]
(modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm_local_stress.i)
# 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
[]