Hourglass Resistance

To correct this phenomenon, it is necessary to introduce anti-hourglass forces and moments. Two possible formulations are presented hereafter.

Flanagan-Belytschko Formulation

Ishell=11

In the Kosloff-Frasier formulation seen in Kosloff and Frasier Formulation, the hourglass base vector Γ I α is not perfectly orthogonal to the rigid body and deformation modes that are taken into account by the one point integration scheme. The mean stress/strain formulation of a one point integration scheme only considers a fully linear velocity field, so that the physical element modes generally contribute to the hourglass energy. To avoid this, the idea in the Flanagan-Belytschko formulation is to build an hourglass velocity field which always remains orthogonal to the physical element modes. This can be written as:(1)
v i I H o u r = v i I v i I L i n
The linear portion of the velocity field can be expanded to give:(2)
v i I H o u r = v i I ( v i I ¯ + v i I x j ( x j x j ¯ ) )
Decomposition on the hourglass base vectors gives 1:(3)
q i α t = Γ I α v i I H o u r = ( v i I v i I x j x j ) Γ I α
Where,
q i α t
Hourglass modal velocities
Γ I α
Hourglass vectors, base
Remembering that v i x j = Φ j x j v i J and factorizing Hourglass Modes, 式 5 gives:(4)
q i α t = v i I ( Γ I α Φ J x j x j Γ J α )
(5)
γ I α = Γ I α Φ J x j x j Γ J α
is the hourglass shape vector used in place of Γ I α in Hourglass Modes, 式 2.


図 1. Flanagan Belytschko Hourglass Formulation


Elastoplastic Hourglass Forces

Ishell=3

The same formulation as elastic hourglass forces is used (Hourglass Elastic Stiffness Forces and Flanagan et al. 1) but the forces are bounded with a maximum force depending on the current element mean yield stress. The hourglass forces are defined as:(6)
f i I h g r = f i h g r Γ I
For in plane mode:(7)
f i h g r ( t + Δ t ) = f i h g r ( t ) + 1 8 h m E t q i t Δ t
(8)
f i h g r ( t + Δ t ) = min ( f i h g r ( t + Δ t ) , 1 2 h m σ y t A )
For out of plane mode:(9)
f i h g r ( t + Δ t ) = f i h g r ( t ) + 1 40 h f E t 3 q i t Δ t
(10)
f i h g r ( t + Δ t ) = min ( f i h g r ( t + Δ t ) , 1 4 h f σ y t 2 )
Where,
t MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyyaaaa@36DC@
Element thickness
σ y
Yield stress
A MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyyaaaa@36DC@
Element area

Physical Hourglass Forces

Ishell=22, 24

The hourglass forces are given by:(11)
f int h g r = Ω e B H t C B H v d Ω e
For in-plane membrane rate-of-deformation, with Φ = ξ η and γ I defined in Bilinear Shape Functions, 式 3:(12)
[ ( B I m ) H ]=[ γ I ϕ,x 0 0 0 γ I ϕ,y 0 γ I ϕ,y γ I ϕ,x 0 ]
For bending:(13)
[ ( B I b ) H ]=[ 0 γ I ϕ,x γ I ϕ,y 0 γ I ϕ,x γ I ϕ,y ]
It is shown in 2 that the non-constant part of the membrane strain rate does not vanish when a warped element undergoes a rigid body rotation. Thus, a modified matrix [ ( B I m ) H ] is chosen using z γ = γ I z I as a measure of the warping:(14)
[ ( B I m ) H ] = [ γ I ϕ , x 0 z γ b x I ϕ , x 0 γ I ϕ , y z γ b y I ϕ , y γ I ϕ , y γ I ϕ , x z γ ( b x I ϕ , y + b y I ϕ , x ) ]
This matrix is different from the Belytschko-Leviathan 3 correction term added at rotational positions, which couples translations to curvatures as:(15)
[ ( B I m ) H ] = [ γ I ϕ , x 0 0 0 1 4 z γ ϕ , x 0 γ I ϕ , y 0 1 4 z γ ϕ , y 0 γ I ϕ , y γ I ϕ , x 0 1 4 z γ ϕ , x 1 4 z γ ϕ , y ]

This will lead to membrane locking (the membrane strain will not vanish under a constant bending loading). According to the general formulation, the coupling is presented in terms of bending and not in terms of membrane, yet the normal translation components in ( B I m ) do not vanish for a warped element due to the tangent vectors t i ( ξ , η ) which differ from t i ( 0 , 0 ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWG0bWaaSbaaSqaaiaadMgaaeqaaOWaaeWaaeaacaaIWaGaaiil aiaaicdaaiaawIcacaGLPaaaaaa@3EA1@ .

Fully Integrated Formulation

Ishell=12

The element is based on the Q4 γ 24 shell element developed in 4 by Batoz and Dhatt. The element has 4 nodes with 5 local degrees-of-freedom per node. Its formulation is based on the Cartesian shell approach where the middle surface is curved. The shell surface is fully integrated with four Gauss points. Due to an in-plane reduced integration for shear, the element shear locking problems are avoided. The element without hourglass deformations is based on Mindlin-Reissner plate theory where the transversal shear deformation is taken into account in the expression of the internal energy. Consult the reference for more details.

Shell Membrane Damping

The shell membrane damping, dm, is only used for LAWS 25, 27, 19, 32 and 36. The Shell membrane damping factor is a factor on the numerical VISCOSITY and not a physical viscosity. Its effect is shown in the formula of the calculation of forces in a shell element:

d m = d m read in D00 input (Shell membrane damping factor parameter) then:(16)
d m = 2 d m D 00 ρ 0 c A R E A
Effect in the force vector ( F MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyyaaaa@36DC@ ) calculation:(17)
F 1 n e w = F 1 o l d + d m ( ε ˙ 11 + ε ˙ 22 2 )
(18)
F 2 n e w = F 2 o l d + d m ( ε ˙ 22 + ε ˙ 11 2 )
(19)
F 3 n e w = F 3 o l d + d m ε ˙ 12 3
Where,
ρ 0
Density
A R E A MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyqaiaadk facaWGfbGaamyqaaaa@3923@
Area of the shell element surface
dt
Time step
c MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyyaaaa@36DC@
Sound speed

In order to calibrate the dm value so that it represents the physical viscosity, one should obtain the same size for all shell elements (Cf. A R E A MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyqaiaadk facaWGfbGaamyqaaaa@3923@ factor), then scale the physical viscosity value to the element size.

1 Flanagan D. and Belytschko T., 「A Uniform Strain Hexahedron and Quadrilateral with Orthogonal Hourglass Control」, Int. Journal Num.Methods in Engineering, 17 679-706, 1981.
2 Belytschko T., Lin J.L. and Tsay C.S., 「Explicit algorithms for the nonlinear dynamics of shells」, Computer Methods in Applied Mechanics and Engineering, 42:225-251, 1984.
3 Belytschko T. and Leviathan I., 「Physical stabilization of the 4-node shell element with one-point quadrature」, Computer Methods in Applied Mechanics and Engineering, 113:321-350, 1992.
4 Batoz J.L. and Dhatt G., 「Modeling of Structures by finite element」, volume 3, Hermes, 1992.