Hourglass Modes

Hourglass modes are element distortions that have zero strain energy. Thus, no stresses are created within the element. There are 12 hourglass modes for a brick element, 4 modes for each of the 3 coordinate directions. Γ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacqqHtoWraaa@3A3F@ represents the hourglass mode vector, as defined by Flanagan-Belytschko. 1 They produce linear strain modes, which cannot be accounted for using a standard one point integration scheme.


Figure 1.
Γ 1 = ( + 1 , 1 , + 1 , 1 , + 1 , 1 , + 1 , 1 )


Figure 2.
Γ 2 = ( + 1 , + 1 , 1 , 1 , 1 , 1 , + 1 , + 1 )


Figure 3.
Γ 3 = ( + 1 , 1 , 1 , + 1 , 1 , + 1 , + 1 , 1 )


Figure 4.

Γ 4 = ( + 1 , 1 , + 1 , 1 , 1 , + 1 , 1 , + 1 )

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

Kosloff and Frasier Formulation

The Kosloff-Frasier hourglass formulation 2 uses a simplified hourglass vector. The hourglass velocity rates are defined as:(1)
q i α t = I = 1 8 Γ I α v i I
Where,
Γ
Non-orthogonal hourglass mode shape vector
ν
Node velocity vector
i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGPbaaaa@39C5@
Direction index, running from 1 to 3
I MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGPbaaaa@39C5@
Node index, from 1 to 8
α
Hourglass mode index, from 1 to 4

This vector is not perfectly orthogonal to the rigid body and deformation modes.

All hourglass formulations except the physical stabilization formulation for solid elements in Radioss use a viscous damping technique. This allows the hourglass resisting forces to be given by:(2)
f i I h g r = 1 4 ρ c h ( Ω 3 ) 2 α q i α t Γ I α
Where,
ρ
Material density
c
Sound speed
h
Dimensional scaling coefficient defined in the input
Ω
Volume

Flanagan-Belytschko Formulation

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:(3)
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:(4)
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 vectors base gives 1:(5)
q i α t = Γ I α v i I H o u r = ( v i I v i l 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 Equation 5 gives:(6)
q i α t = v i I ( Γ I α Φ j x j x j Γ I α )
(7)
γ I α = Γ I α Φ j x j x j Γ J α

is the hourglass shape vector used in place of Γ I α in Equation 2.

Physical Hourglass Formulation

You also try to decompose the internal force vector as:(8)
{ f I int } = { ( f I int ) 0 } + { ( f I int ) H }
In elastic case, you have:(9)
{ f I int } = Ω [ B I ] t [ C ] j = 1 8 [ B J ] { v J } d Ω = Ω ( [ B I ] 0 + [ B ¯ I ] H ) t [ C ] j = 1 8 ( [ B J ] 0 + [ B ¯ J ] H ) { v J } d Ω

The constant part { ( f I int ) 0 } = Ω ( [ B I ] 0 ) t [ C ] j = 1 8 [ B J ] 0 { v J } d Ω is evaluated at the quadrature point just like other one-point integration formulations mentioned before, and the non-constant part (Hourglass) will be calculated as:

Taking the simplification of x i ξ j = 0 ; ( i j ) (that is the Jacobian matrix of Strain Rate, Equation 1 is diagonal), you have:(10)
( f i I int ) H = α = 1 4 Q i α γ I α
with 12 generalized hourglass stress rates Q . i α calculated by:(11)
Q . i i = μ [ ( H j j + H k k ) q ˙ i i + H i j q ˙ j j + H i k q ˙ k k ] Q . j j = μ [ 1 1 ν H i i q ˙ i j + ν ¯ H i j q ˙ j i ] Q . i 4 = 2 μ 1 + ν 3 H i i q ˙ j 4
and(12)
H ii = Ω ( ϕ j x i ) 2 dΩ= Ω ( ϕ k x i ) 2 dΩ=3 Ω ( ϕ 4 x i ) 2 dΩ H ij = Ω ϕ i x j ϕ j x i dΩ

Where, i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGPbaaaa@39C5@ , j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGPbaaaa@39C5@ , k MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGPbaaaa@39C5@ are permuted between 1 to 3 and q ˙ i α has the same definition than in Equation 6.

Extension to nonlinear materials has been done simply by replacing shear modulus μ by its effective tangent values which is evaluated at the quadrature point. For the usual elastoplastic materials, use a more sophistic procedure which is described in Advanced Elasto-plastic Hourglass Control.

Advanced Elasto-plastic Hourglass Control

With one-point integration formulation, if the non-constant part follows exactly the state of constant part for the case of elasto-plastic calculation, the plasticity will be under-estimated due to the fact that the constant equivalent stress is often the smallest one in the element and element will be stiffer. Therefore, defining a yield criterion for the non-constant part seems to be a good idea to overcome this drawback.
Plastic yield criterion
The von Mises type of criterion is written by:(13)
f = σ e q 2 ( ξ , η , ζ ) σ y 2 = 0
for any point in the solid element, where σ y is evaluated at the quadrature point.
As only one criterion is used for the non-constant part, two choices are possible:
  1. taking the mean value, that is, f = f ( σ ¯ e q ) ; σ ¯ e q = 1 Ω Ω σ e q d Ω
  2. taking the value by some representative points, for example: eight Gauss points
The second choice has been used in this element.
Elastro-plastic hourglass stress calculation
The incremental hourglass stress is computed by:
  • Elastic increment

    ( σ i ) n + 1 t r H = ( σ i ) n H + [ C ] { ε ˙ } H Δ t

  • Check the yield criterion
  • If f 0 , the hourglass stress correction will be done by un radial return

    ( σ i ) n + 1 H = P ( ( σ i ) n + 1 t r H , f )

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 Kosloff D. and Frazier G., “Treatment of hourglass pattern in low order finite element code”, International Journal for Numerical and Analytical Methods in Geomechanics, 1978.