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@+=feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKfMBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhiov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake 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.


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


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


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


図 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)
qiαt=I=18ΓIαviI
Where,
Γ
Non-orthogonal hourglass mode shape vector
ν
Node velocity vector
i MathType@MTEF@5@5@+=feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKfMBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhiov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaakeaacaWGPbaaaa@39C5@
Direction index, running from 1 to 3
I MathType@MTEF@5@5@+=feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKfMBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhiov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake 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)
fiIhgr=14ρch(Ω3)2αqiα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)
viIHour=viIviILin
The linear portion of the velocity field can be expanded to give:(4)
viIHour=viI(viI¯+viIxj(xjxj¯))
Decomposition on the hourglass vectors base gives 1:(5)
qiαt=ΓIαviIHour=(viIvilxjxj)ΓIα
Where,
qiαt
Hourglass modal velocities
ΓIα
Hourglass vectors base
Remembering that vixj=ΦjxjviJ and factorizing 式 5 gives:(6)
qiαt=viI(ΓIαΦjxjxjΓIα)
(7)
γIα=ΓIαΦjxjxjΓJα

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

Physical Hourglass Formulation

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

The constant part {(fIint)0}=Ω([BI]0)t[C]j=18[BJ]0{vJ}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 xiξj=0;(ij) (that is the Jacobian matrix of Strain Rate, 式 1 is diagonal), you have:(10)
(fiIint)H=α=14QiαγIα
with 12 generalized hourglass stress rates Q.iα calculated by:(11)
Q.ii=μ[(Hjj+Hkk)q˙ii+Hijq˙jj+Hikq˙kk]Q.jj=μ[11νHiiq˙ij+ν¯Hijq˙ji]Q.i4=2μ1+ν3Hiiq˙j4
and(12)
Hii=Ω(ϕjxi)2dΩ=Ω(ϕkxi)2dΩ=3Ω(ϕ4xi)2dΩHij=ΩϕixjϕjxidΩ

Where, i MathType@MTEF@5@5@+=feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKfMBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhiov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGPbaaaa@39C5@ , j MathType@MTEF@5@5@+=feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKfMBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhiov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGPbaaaa@39C5@ , k MathType@MTEF@5@5@+=feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKfMBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhiov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGPbaaaa@39C5@ are permuted between 1 to 3 and q˙iα has the same definition than in 式 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=σeq2(ξ,η,ζ)σy2=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(σ¯eq);σ¯eq=1ΩΩσeqdΩ
  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+1trH=(σi)nH+[C]{ε˙}HΔt

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

    (σi)n+1H=P((σi)n+1trH,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.