One-Point Quadrature Shell Element

In this section, a one-point quadrature shell element formulation will be developed from the general formulation described in the previous section. It is based on the Physical Stabilization method which explicitly computes the stabilization terms in making some simplifications.

The following formulations will be written in the local coordinate system [ t 1 t 2 n ] (the circumflex in the co-rotational system notation has been omitted for convenience).

Kinematic Approximation

The velocity interpolation using the nodal tangent vectors ( t 1 I , t 2 I ) complicates the strain computation, especially for transverse shear which is used mainly as a penalty function. To be consistent with the one-point quadrature approach, the kinematic approximation is performed by:(1) v i = N I [ v i I + z ( δ i 2 ω ¯ 1 I + δ i 1 ω ¯ 2 I ) ]
Where, ω ¯ i I ( i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamyAaaaa@36E4@ =1,2) is the nodal rotation velocity around t i I . ω ¯ i I can be computed by a projection scheme by:(2) ω ¯ I = ω I ( ω I n I ) n I = P ( n I ) ω I

The projection consists in eliminating the nodal drilling rotations in order to reinforce Mindlin's kinematic condition at the nodes. It has been pointed out 1 that without this projection, the element is too flexible and cannot pass the Twisted Beam test.

This projection has a drawback of changing the invariance property to rigid body motion, that is, a warped element being invariant to rigid body rotation will now strains under to rigid body rotation if the drill projection is applied. To overcome this problem, a full projection proposed by Belytschko-Leviathan 2 which free either drilling rotation or rigid motion should be used. This full projection is only used for QEPH element.

In-plane Strain-rate Construction

Constant part

It is useful to write the shape functions in Belytschko-Bachrach's mixed form:(3) N I ( x , y , ξ η ) = Δ I + b x I x + b y I y + γ I ϕ

with:

Δ I = [ t I ( t I x I ) b x I ( t I y I ) b y I ] ; t I = ( 1 , 1 , 1 , 1 )

b x I = ( y 24 y 31 y 42 y 13 ) / A ; b y I = ( x 42 x 13 x 24 x 31 ) / A ( f i j = ( f i f j ) / 2 )

γ I = [ h I ( h J x J ) b x I ( h J y J ) b y I ] / 4 ; ϕ = ξ η ;

A is the area of element.

The derivation of the shape functions is given by:

N I , α = b α I + γ I ϕ , α ( α = x , y )

Where,(4) ϕ , x = η ξ , x + ξ η , x = J 11 1 η + J 12 1 ξ ϕ , y = η ξ , y + ξ η , y = J 21 1 η + J 22 1 ξ

The advantage of this shape function form is that a linear field expressed with Cartesian coordinates and a bilinear field expressed with Natural coordinates is decomposed so that the constant part is directly formulated with the Cartesian coordinates, and the non-constant part is to be approached separately.

The in-plane rate-of-deformation (decomposed on membrane and bending) is given by:(5) { D } = { D m } + z { D b } = [ B I ] { v n I }

with:

B I = [ B I m ] + z [ B I b ] ; < v n I > = < v x I v y I v z I ϖ x I ϖ y I > .
The development of the general formulations leads to the constant part, denoted by superscript 0, of the matrix [ B I ]:(6) [ ( B I m ) 0 ] = [ b x I 0 0 0 b y I 0 b y I b x I 0 ] [ ( B b ) 0 ] = 2 H [ ( B I m ) 0 ] + [ b x I c 0 0 0 b x I 0 b y I c 0 b y I 0 b y I c b x I c 0 b x I b y I ]

with

< b x I c > = 4 z γ / A 2 < b y 2 ; b y 1 ; b y 4 ; b y 3 >

< b y I c > = 4 z γ / A 2 < b x 4 ; b x 3 ; b x 2 ; b x 1 >

The parameter z γ = γ I z I is a measure of the warping of the element.

The first term 2 H [ ( B I m ) 0 ] is neglectable. You have verified that the order of 2 H [ ( B I m ) 0 ] is ε times the second term of [ ( B b ) 0 ] with ε = 2 x ¯ 34 / y ¯ 34 ( f ¯ i j = ( f i + f j ) / 2 ) , which vanishes when the element is rectangular. Thus, this term is not used in the program.

The constant part of the in-plane rate-of-deformation formulation without the H term is consistent with the result of Belytschko's family shell element 4, 6though this part has been obtained in a very different manner. Letellier has given the same result in his thesis 5, and studies were also made of the quadratic terms with respect to z .

Non-constant Part

The main simplification for the non-constant part formulation, in order to overcome the difficulties described above, is:

The element is considered to be flat.

In this case, the Jacobian matrix is written as:(7) [ J ] = [ F 0 ] t = [ x , ξ y , ξ 0 x , η y , η 0 0 0 a J ¯ 0 / 2 J ¯ ] = [ [ J ¯ ] 0 0 a J ¯ 0 / 2 J ¯ ]
with the determinant J ¯ of the in-plane Jacobian:(8) J ¯ = det [ J ¯ ] = J ¯ 0 + J ¯ 1 ξ + J ¯ 2 η

and:

J ¯ 0 = [ ( ξ I x I ) ( η I y I ) ( η I x I ) ( ξ I y I ) ] / 16 = A / 4

J ¯ 1 = [ ( ξ I x I ) ( h I y I ) ( h I x I ) ( ξ I y I ) ] / 16

J ¯ 2 = [ ( h I x I ) ( η I y I ) ( η I x I ) ( h I y I ) ] / 16
The inverse of the in-plane Jacobian matrix can be expressed explicitly:(9) [ J ] = [ ξ , x η , x ξ , y η , y ] = 1 J [ y , η y , ξ x , η x , ξ ] = 1 4 J ¯ [ ( η I + h I ξ ) y I ( ξ I + h I η ) y I ( η I + h I ξ ) x I ( ξ I + h I η ) x I ]
You can now write the non-constant part, denoted by superscript H , of the matrix [ B I ] for in-plane rate-of-deformation:(10) [ ( B I m ) H ] = [ γ I ϕ , x 0 0 0 γ I ϕ , y 0 γ I ϕ , y γ I ϕ , x 0 ] ; [ ( B I b ) H ] = [ 0 γ I ϕ , x γ I ϕ , y 0 γ I ϕ , x γ I ϕ , y ]
It is shown in 3 that the non-constant part of membrane rate-strain does not vanish when a warped element undergoes a rigid body rotation. Thus, a modified matrix [ ( B I m ) H ] is chosen:(11) [ ( 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 correction term added at rotational positions, which couples translations to curvatures:(12) [ ( 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 the bending terms not in the membrane terms, yet the normal translation components in ( B I m ) i3 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aadaqadaqaaiaadkeadaqhaaWcbaGaamysaaqaaiaad2gaaaaakiaa wIcacaGLPaaadaWgaaWcbaGaamyAaiaaiodaaeqaaaaa@3EF5@ do not vanish for a warped element due to the tangent vectors t i ( ξ , η ) which differ from t i (0,0).

Out-plane Strain-rate Construction

The out-plane rate-of-deformation (transverse shear) is interpolated by the Dvorkin-Bathe method, whose closed form is given by Belytschko-Leviathan:(13) { D x z D y z } = [ B I c ] { v z I ϖ x I ϖ y I }
Where, [ B I c ] ( ζ , η ) = [ B I c η ] ( η ) + [ B I c ξ ] ( ξ ) (14) ( B I c η ) 0 = 1 16 A [ 4 η I y I ξ I ( ξ J + h J η I ) η k y k y J ( ξ J + h J η I ) η k y k x J 4 η I x I ξ I ( ξ J + h J η I ) η k x k y J ( ξ J + h J η I ) η k x k x J ] ( B I c ξ ) 0 = 1 16 A [ 4 ξ I y I η I ( η J + h J ξ I ) ξ k y k y J ( η J + h J ξ I ) ξ k y k x J 4 ξ I y I η I ( η J + h J ξ I ) ξ k x k y J ( η J + h J ξ I ) ξ k x k x J ] and ( B I c η ) H = η I ( B I c η ) 0 η ; ( B I c ξ ) H = ξ I ( B I c ξ ) 0 ξ

The straightforward form of [ B Ic ] is obtained using one additional simplification:

ξ , x ( ξ , η ) = ξ , x ( 0 , 0 ) η , x ( ξ , η ) = η , x ( 0 , 0 )

which is true for a parallelogram element. Although this simplification is not necessary, it is justified by the fact that the transverse shear terms serve mainly as a penalty function.

Explicit Integration of the Nodal Internal Force Vector

Elastic Case

The elementary nodal force vector is computed by:

{ f I int } = v e [ B I ] t [ C ] [ B J ] { v n J } d v e

Taking advantage of substantial orthogonality between the:
  • Constant in-plane fields along with the non-constant ones
  • Membrane and the bending
  • In-plane fields and the out-plane fields
  • Decomposed non-constant out-plane fields

Resulting in:

{ f I int } = { ( f I int ) 0 } + { ( f I int ) H }

With { ( f I int ) 0 } the constant part being computed with one-point quadrature, and

{ ( f I int ) H } = v e [ ( B I ) H ] t [ C ] [ ( B J ) H ] { v n J } d v e

It can be shown in the last equation that only the following scalar functions need to be integrated:(15) H xx = A e ϕ ,x ϕ ,x d A e , H xy = A e ϕ ,x ϕ ,y d A e , H yy = A e ϕ ,y ϕ ,y d A e

These can be evaluated explicitly.

Defining 6 hourglass generalized rate-of-deformation q ˙ by:(16) membrane : { q ˙ x m = γ I v x I + z 1 b x I v z I q ˙ y m = γ I v y I + z 1 b y I v z I (17) bending : { q ˙ x b = γ I ϖ y I q ˙ y b = γ I ϖ x I (18) shear : { q ˙ x s = 4 h I v z I [ ( ξ I y I ) η I + ( h I y I ) t I ] ϖ x I + [ ( ξ I x I ) η I + ( h I x I ) t I ] ϖ y I q ˙ y s = 4 h I v z I [ ( η I y I ) ξ I + ( h I y I ) t I ] ϖ x I + [ ( η I x I ) ξ I + ( h I x I ) t I ] ϖ y I

The rate-of-deformations will be written explicitly.

The rate form of the constitutive relation is written as (stress plane for in-plane terms):

σ ˙ = C : D

With the assumption: the spin is constant within the element, the objectivity principle will be satisfied. The incremental computation is performed with the hourglass generalized rate-of-deformation q ˙ :

( q ) n+1 = ( q ) n + q ˙ n+ 1 2 Δt

Noting that if ϕ , α is considered as constant over a time step, it is equivalent to the incremental stress computation.

Physical Nonlinear Case

Now consider an elastoplastic problem.

The elementary nodal internal force vector is now computed by:

{ f I int } = v e [ B I ] t { σ } d v e

The constitutive relation is written by either a tangent form: σ ˙ = C t : D , or a projection form: σ = P ( σ e , ... )

Where, C t ( σ , ... ) is the history-dependent tangent tensor; { σ n + 1 e } = { σ n } + [ C ] { D n + 1 / 2 } Δ t is the trial stress, and the function P MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGqbaaaa@39AC@ consists of projecting the trial stress on the updated yield surface.

The decomposed form of the last equation is written:

{ f I int } = { ( f I int ) 0 } + { ( f I int ) H }

The constant part is computed by integration at each integration point through the thickness.

The stabilization part will be approached by relying on two hypotheses:
  • Keep the same orthogonalities as in the elastic case, and
  • Use the unidimensional tangent modulus E t ( ζ ) to evaluate the non-constant rate-stress, that is,(19) { σ ˙ } H = ( E t ( ζ ) / E ) [ C ] { D } H
    Where,
    E MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake aacaWGfbaaaa@39A1@
    Young's modulus
    [ C ]
    Matrix of elastic moduli

Thus, the elastic case easily extends to the nonlinear case.

The incremental computation with the hourglass generalized rate-of-deformation q ˙ becomes:(20) m e m b r a n e : ( q m ) n + 1 = ( q m ) n + λ m q ˙ n + 1 / 2 m Δ t (21) b e n d i n g : ( q b ) n + 1 = ( q b ) n + λ b q ˙ n + 1 / 2 b Δ t (22) s h e a r : ( q s ) n + 1 = ( q s ) n + q ˙ n + 1 / 2 s Δ t

Where, E t ( ζ ) is obtained by the constant stress incremental computation along the thickness and E t ( ζ ) = E in the elastic zone, and λ m = E ¯ t / E ; E ¯ t is the average value of E t ( ζ ) and λ b = min E t ( ζ ) / E .

For the QPH shell, λ m MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbwvMCKf MBHbqefqvATv2CG4uz3bIuV1wyUbqedmvETj2BSbqefm0B1jxALjhi ov2DaebbnrfifHhDYfgasaacH8srps0lbbf9q8WrFfeuY=Hhbbf9v8 qqaqFr0xc9pk0xbba9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9 q8qqQ8frFve9Fve9Ff0dmeaacaGacmGadaWaaiqacaabaiaafaaake GabaaVqiabeU7aSnaaBaaaleaacaWGTbaabeaaaaa@3C78@ = 1

The key orthogonalities has been maintained without any significant deterioration in performance, although the first two orthogonalities might have been slightly violated. In fact, it is simply due to these orthogonalities that a one-point quadrature element dramatically reduces the computation cost; otherwise you return to the full integration scheme.

Most of the physical stabilization elements have incorporated the following assumption.

The material response is constant within the element.

There are two alternatives to this assumption:
  • Take the elastic matrix [ C ] directly:

    { σ ˙ } = { σ ˙ } 0 + [ C ] { ( ε ˙ ε ˙ P ) H } = { σ ˙ } 0 + [ C ] { ( ε ˙ ) H } { ( ε ˙ P ) H } = 0

    which means that the plastic rate-of-deformation { ε ˙ P } is constant within the element; or

  • Take the tangent matrix [ C t ] ( ξ = η = 0 , ζ = constant).

Since the components of [ C t ] are generally functions of the updated stress (precisely, the stress deviator for an elastoplastic problem with an associative flow rule, which means that the stress is constant within the element.

Neither of these alternatives is theoretically perfect. Note that the [ C ] option results in a contradiction with the stress computation (which yields different results for the constant part and the non-constant part); it is more expensive and the tangent form is not generally used for constant stress computations within an explicit program. Hence, the approximation based on the above assumption is not necessary.

The choice of the moduli for the nonlinear case has not been studied for Belytschko-Leviathan's element 6, and it has been shown that this choice has little effect on the result of the "Cylindrical panel test". In 7, the elastic tangent matrix has been used for the evaluation of the stabilizing forces.

QPH,QPPS have shown often stiffer behaviors and sometimes have certain numerical problems in crash simulations.

1 Belytschko T. and Leviathan I., 「Projection schemes for one-point quadrature shell elements」, Computer Methods in Applied Mechanics and Engineering, 115:277-286, 1993.
2 Belytschko Ted and Leviathan Itai, 「Projection schemes for one-point quadrature shell elements」, Computer Methods in Applied Mechanics and Engineering, Vol. 115, 227-286, 1994.
3 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.
4 Belytschko T., Wong B.L. and Chiang H.Y.「Advances in one-point quadrature shell elements」, Computer Methods in Applied Mechanics and Engineering, 96:93-107, 1989.
5 Letellier Antoine, 「Contribution to the modeling of impacts of birds on the blades of aircraft engines」, PhD thesis, University d'Evry, 1996.
6 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.
7 Zhu Y. and Zacharia T., 「A new one-point quadrature, quadrilateral shell element with drilling degree of freedom」, Computer Methods in Applied Mechanics and Engineering, 136:165-203, 1996.