General Degenerated 4-Node Shell Formulations

The following formulations of degenerated quadrilateral shells are based on the successful full integration element MITC4 developed by Dvorkin-Bathe ¹ and Q4

γ

24 developed by Batoz and Dhatt ²; they are suitable for both thin and thick shells and are applicable to linear and nonlinear problems. Their main feature is that a classical displacement method is used to interpolate the in-plane strains (membrane, bending), and a mix/collocation (or assumed strain) method is used to interpolate the out-plane strains (transverse shear). Certain conditions are also specified:

They are based on the Reissner-Mindlin model,
In-plane strains are linear, out-plane strains(transverse shear) are constant throughout the thickness,
Thickness is constant in the element (the normal and the fiber directions are coincident),
5 DOF in the local system (that is, the nodal normal vectors are not constant from one element to another).

Notational Conventions

A bold letter denotes a vector or a tensor.
An upper case index denotes a node number; a lower case index denotes a component of vector or tensor.
The Einstein convention applies only for the repeated index where one is subscript and another is superscript, e.g.:
$N_{I} x^{I} = \sum N_{I} x^{I}$
{} denotes a vector and [ ] denotes a matrix.

Geometry and Kinematics

The geometry of the 4-node degenerated shell element, as shown in 図 1, is defined by its mid-surface with coordinates denoted by

X p

interpolated by the node coordinates

X^{I}

(

I

=1,4):(1)

X p (ξ, η) = N_{I} X^{I}

Where,

N_{I} (ξ, η)

are the bilinear isoparametric shape functions, given by:(2)

N_{I} (ξ, η) = (1 + ξ_{I} ξ) (1 + η_{I} η) / 4

A generic point

q

within the shell is derived from point p on the mid-surface and its coordinate along the normal (fiber):(3)

X q (ξ, η, ζ) = X p + z \cdot n

with $z = \frac{a ζ}{2}$

Where,

$a$: Shell thickness

The transformation between the Cartesian system and the Natural system is given by the differential relation (in matrix form):(4)

{d X q} = ([F_{0}] + z [{n,_{ξ}} {n,_{η}} {0}] {d ξ} = [F] {d ξ}

with $[F_{0}] = [{g_{1}} {g_{2}} \frac{a}{2} {n}]; g_{1} (ξ, η) = X_{p, ξ} = N_{I, ξ} X^{I}, g_{2} (ξ, η) = X_{p, η} = N_{I, η} X^{I}$

$F$ is the gradient tensor which is related to the Jacobian tensor $[J] = {[F]}^{t}$ .

With 5 DOF at each node I (three translational velocities

v^{I} \to v_{i}^{I}

and two rotational velocities

(ω^{I} \to ω_{i}^{I})

, the velocity interpolation is given by the Mindlin model:(5)

v_{q} (ξ, η, ζ) = v_{p} + z β = N_{I} (ξ, η) (v^{I} + z β^{I})

β = ω \times n

Where, $β$ and $ω$ are the rotational velocity vectors of the normal: $β = β_{1} t_{1} + β_{2} t_{2} = ω_{2} t_{1} - ω_{1} t_{2}$

and $(t_{1}, t_{2}, n)$ is base of the local coordinate system.

式 5 can be written also by:

$v_{i} = N_{I} v_{i}^{I} + z (- N_{I} ω_{1}^{I} t_{2 i}^{I} + N_{I} ω_{2}^{I} t_{1 i}^{I}); i = 1, 3$

This velocity interpolation is expressed in the global system, but $ω^{I}$ must be defined first in the local nodal coordinate system to ensure Mindlin's kinematic condition.

Strain-Rate Construction

The in-plane rate-of-deformation is interpolated by the usual displacement method.

The rate-of-deformation tensor (or velocity-strains)

D = (L_{t} + L_{t}^{T}) / 2

is defined by the velocity gradient tensor

L

:(6)

{d v_{q}} = [L_{ζ}] {d ξ} = [L_{ζ}] {[F]}^{- 1} {d X q} = [L] {d X q} = [Q] [L_{t}] {[Q]}^{T} {d X q}

with $[Q] = [{t_{1}} {t_{2}} {n}]$ .

The Reissner-Mindlin conditions $ε_{ζ} = 0$ and $σ_{ζ} = 0$ requires that the strain and stress tensors are computed in the local coordinate system (at each quadrature point).

After the linearization of $L_{t}$ with respect to $z$ , the in-plane rate-of-deformation terms are given by:

${[L_{t}]}_{2 \times 2} = [L_{s 0}] + z [L_{s 1}]$

with the membrane terms:

$[L_{s 0}] = [C_{1}] [C_{0}]$

the bending terms:(7)

[L_{s 1}] = 2 H [L_{s 0}] + [C_{1}] [\begin{matrix} g^{2} \cdot n,_{η} & - g^{1} \cdot n,_{η} \\ - g^{2} \cdot n,_{ξ} & g^{1} \cdot n,_{ξ} \end{matrix}] [C_{0}] + [\begin{matrix} t_{1} \cdot β,_{ξ} & t_{1} \cdot β,_{η} \\ t_{1} \cdot β,_{ξ} & t_{1} \cdot β,_{η} \end{matrix}]

Where, the contravariant vectors

g^{α} (α = 1, 2)

, dual to

g_{α}

, satisfy the orthogonality condition:

g^{α}_{\cdot}

g_{β} = δ_{β}^{α}

(Kronecker delta symbol);

H (ξ, η)

is the average curvature:

2 H = - (g^{1}_{\cdot} n,_{ξ} + g^{2}_{\cdot} n,_{η})

(8)

[C_{0}] = [\begin{matrix} t_{1} \cdot g^{1} & t_{2} \cdot g^{1} \\ t_{1} \cdot g^{2} & t_{2} \cdot g^{2} \end{matrix}] and [C_{1}] = [\begin{matrix} t_{1} \cdot v_{p},_{ξ} & t_{1} \cdot v_{p},_{η} \\ t_{2} \cdot v_{p},_{ξ} & t_{2} \cdot v_{p},_{η} \end{matrix}]

The curvature-translation coupling is presented in the bending terms for a warped element (the first two terms in the last equation.)

The out-plane rate-of-deformation (transverse shear) is interpolated by the "assumed strain" method, which is based on the Hu-Washizu variation principle.

If the out-plane rate-of-deformation is interpolated in the same manner for a full integration scheme, it will lead to shear locking. It is known that the transverse shear strains energy cannot vanish when it is subjected to a constant bending moment. Dvorkin-Bathe's ¹ mix/collocation method has been proved very efficient in overcoming this problem. This method consists in interpolating the transverse shear from the values of the covariant components of the transverse shear strains at 4 mid-side points. That is:(9)

γ_{ξ} = [(1 - η) γ_{ξ}^{A 1} + (1 + η) γ_{ξ}^{A 2}] / 2

(10)

γ_{η} = [(1 - ξ) γ_{η}^{B 1} + (1 + ξ) γ_{η}^{B 2}] / 2

Where,

γ_{ξ}^{A α}, γ_{η}^{B α}

are the values of the covariant components at 4 mid-side points which vanish under a constant bending moment (図 2).

Special Case for One-point Quadrature and the Difficulties in Stabilization

The formulations described above are general for both the full integration and reduced integration schemes. For a one-point quadrature element, you have the following particularities:

The quadrature point is often chosen at

(ξ = 0, η = 0)

. The derivatives of the shape functions are:(11)

N_{I, ξ} = (ξ_{I} + h_{I} η) / 4

N_{I, η} = (η_{I} + h_{I} ξ) / 4

Where, $h_{I} = (\begin{matrix} 1 & - 1 & 1 & - 1 \end{matrix})$ .

This implies that all the terms computed at the quadrature point are the constant parts with respect to $(ξ, η)$ , and the stabilizing terms (hourglass) are the non-constant parts.

The constant parts can be derived directly from the general formulations at the quadrature point without difficulty. The difficulties in stabilization lie in correctly computing the internal force vector (or stiffness matrices):(12)

f^{int} = \int_{Ω} {\bar{B}}^{t} \cdot σ d Ω = f^{int} (ξ = η = 0) + f_{s t a b}^{int}

It would be ideal if the integration term $f_{s t a b}^{int}$ could be evaluated explicitly. But such is not the case, and the main obstacles are the following:

For a non-coplanar element, the normal varies at each point so that it is difficult to write the non-constant part of strains explicitly. For a physically nonlinear problem, the non-constant part of stress is not generally in an explicit form. Thus, simplification becomes necessary.

¹ Dvorkin E. and Bathe K.J.「A continuum mechanics four-node shell element for 35 general nonlinear analysis」, Engrg Comput, 1:77-88, 1984.

² Batoz J.L. and Dhatt G., 「Modeling of Structures by finite element」, volume 3, Hermes, 1992.