- 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
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>
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.
- 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
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.
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
Tagging 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>
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>
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
Input Files
- (modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment_free_orientation_inclined_hht.i)
- (modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD.i)
- (modules/solid_mechanics/test/tests/shell/static/inclined_straintest.i)
- (modules/solid_mechanics/test/tests/shell/static/pinched_cylinder_symm.i)
- (modules/solid_mechanics/test/tests/shell/static/straintest.i)
- (modules/solid_mechanics/test/tests/shell/static/straintest_shear.i)
- (modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment.i)
- (modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment_free_orientation_inclined.i)
- (modules/solid_mechanics/test/tests/shell/static/large_strain_m_40_AD.i)
- (modules/solid_mechanics/test/tests/shell/static/beam_bending_moment_AD_2.i)
- (modules/solid_mechanics/test/tests/shell/dynamics/shell_dynamics_bending_moment_free.i)
- (modules/solid_mechanics/test/tests/shell/static/tapered.i)
- (modules/solid_mechanics/test/tests/shell/static/pressure_error.i)
- (modules/solid_mechanics/test/tests/shell/static/plate_bending2.i)
- (modules/solid_mechanics/test/tests/shell/static/finite_straintest.i)
- (modules/solid_mechanics/test/tests/shell/static/plate_bending.i)
(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 = UserForcingFunctionNodalKernel
variable = disp_z
boundary = '2'
function = 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/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/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/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. FEM results from different mesh discretizations are
# presented below. Only the 10x10 mesh is included as a test.
# Mesh of 1/8 cylinder | FEM/analytical (Moose) | FEM/analytical (Dvorkin)
# |ratio of normalized disp.| ratio of normalized disp.
#----------------------|-------------------------|-------------------------
# 10 x 10 | 0.806 | 0.83
# 20 x 20 | 1.06 | 0.96
# 40 x 40 | 0.95 | -
# 80 x 160 | 0.96 | -
# The results from FEM analysis matches well with the series solution and with
# the solution presented by Dvorkin and Bathe (1984).
[Mesh]
[./mesh]
type = FileMeshGenerator
file = pinched_cyl_10_10.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
[../]
[]
[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 = 'CD AB'
value = 0.0
[../]
[./simply_support_rot_x]
type = DirichletBC
variable = rot_x
boundary = 'CD BC'
value = 0.0
[../]
[./simply_support_rot_y]
type = DirichletBC
variable = rot_y
boundary = 'CD AD'
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
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_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/straintest_shear.i)
# Test for the shear stress and strain output for 2D planar shell with uniform mesh.
# A cantiliver 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/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 = UserForcingFunctionNodalKernel
variable = disp_z
boundary = 'top'
function = 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/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 = UserForcingFunctionNodalKernel
variable = disp_z
boundary = '2'
function = 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/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 = 24.85069 m; relative error = 3.54 %
# disp_y = -9.125937 m; relative error = 5.19 %
[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 = UserForcingFunctionNodalKernel
variable = disp_z
boundary = top
function = 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
line_search = 'none'
nl_rel_tol = 1e-12
nl_abs_tol = 1e-14
dt = 0.1
dtmin = 0.1
end_time = 3.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/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 = '0.0 0.0 10.0'
variable = rot_y
[../]
[./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/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 = UserForcingFunctionNodalKernel
variable = disp_z
boundary = 'top'
function = 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/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/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/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/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/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
[]