# 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 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, i.e.: 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

with:

${\text{\Delta}}_{I}=\left[{t}_{I}-\left({t}_{I}{x}^{I}\right)\text{\hspace{0.17em}}{b}_{xI}-\left({t}_{I}{y}^{I}\right){b}_{yI}\right]\text{\hspace{0.17em}};\text{\hspace{1em}}{t}_{I}=\left(1,1,1,1\right)$

${b}_{xI}=\left({y}_{24}{y}_{31}{y}_{42}{y}_{13}\right)/A;{b}_{yI}=\left({x}_{42}{x}_{13}{x}_{24}{x}_{31}\right)/A\text{\hspace{1em}}\left(fij=\left(fi-fj\right)/2\right)$

${\gamma}_{I}=\left[{h}_{I}-\left({h}_{J}{x}^{J}\right){b}_{xI}-\left({h}_{J}{y}^{J}\right){b}_{yI}\right]/4\text{\hspace{0.17em}};\text{\hspace{0.17em}}\varphi =\xi \eta \text{\hspace{1em}};$

$A$ is the area of element.

The derivation of the shape functions is given by:

${N}_{I,\alpha}={b}_{\alpha I}+{\gamma}_{I}{\varphi}_{,\alpha}\left(\alpha =x,y\right)$

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.

with:

$\lfloor {B}_{I}\rfloor =\left[{B}_{I}^{m}\right]+z\left[{B}_{I}^{b}\right]\text{\hspace{0.17em}};\text{\hspace{1em}}<{v}_{n}^{I}>=<{v}_{x}^{I}{v}_{y}^{I}{v}_{z}^{I}{\varpi}_{x}^{I}{\varpi}_{y}^{I}>.$with

$<{b}_{xI}^{c}>=4{z}_{\gamma}/{A}^{2}<{b}_{y2}\text{\hspace{0.17em}};\text{\hspace{0.17em}}{b}_{y1}\text{\hspace{0.17em}};\text{\hspace{0.17em}}{b}_{y4}\text{\hspace{0.17em}};\text{\hspace{0.17em}}{b}_{y3}>$

$<{b}_{yI}^{c}>=4{z}_{\gamma}/{A}^{2}<{b}_{x4}\text{\hspace{0.17em}};\text{\hspace{0.17em}}{b}_{x3}\text{\hspace{0.17em}};\text{\hspace{0.17em}}{b}_{x2}\text{\hspace{0.17em}};\text{\hspace{0.17em}}{b}_{x1}>$

The parameter ${z}_{\gamma}={\gamma}_{I}{z}^{I}$ is a measure of the warping of the element.

The first term $2H\left[{\left({B}_{I}^{m}\right)}^{0}\right]$ is neglectable. You have verified that the order of $2H\left[{\left({B}_{I}^{m}\right)}^{0}\right]$ is $\epsilon $ times the second term of $\left[{\left({B}^{b}\right)}^{0}\right]$ with $\epsilon =2{\overline{x}}_{34}/{\overline{y}}_{34}\text{\hspace{1em}}\left({\overline{f}}_{ij}=\left({f}_{i}+{f}_{j}\right)/2\right)\text{\hspace{0.17em}},$ 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}, ^{6}though 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.

and:

${\overline{J}}_{0}=\left[\left({\xi}_{I}{x}^{I}\right)\left({\eta}_{I}{y}^{I}\right)-\left({\eta}_{I}{x}^{I}\right)\left({\xi}_{I}{y}^{I}\right)\right]/16=A/4$

${\overline{J}}_{1}=\left[\left({\xi}_{I}{x}^{I}\right)\left({h}_{I}{y}^{I}\right)-\left({h}_{I}{x}^{I}\right)\left({\xi}_{I}{y}^{I}\right)\right]/16$

${\overline{J}}_{2}=\left[\left({h}_{I}{x}^{I}\right)\left({\eta}_{I}{y}^{I}\right)-\left({\eta}_{I}{x}^{I}\right)\left({h}_{I}{y}^{I}\right)\right]/16$^{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 $\left[{\left({B}_{I}^{m}\right)}^{H}\right]$ is chosen:

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 $${\left({B}_{I}^{m}\right)}_{i3}$$ do not vanish for a warped element due to the tangent vectors $${t}_{i}$$ $(\xi ,\eta )$ which differ from $${t}_{i}$$ (0,0).

## Out-plane Strain-rate Construction

The straightforward form of [${B}_{\mathrm{Ic}}$] is obtained using one additional simplification:

${\xi}_{,x}\left(\xi ,\eta \right)={\xi}_{,x}\left(0,0\right){\eta}_{,x}\left(\xi ,\eta \right)={\eta}_{,x}\left(0,0\right)$

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:

$\left\{{f}_{I}^{\mathrm{int}}\right\}={{\displaystyle {\int}_{ve}\left[{B}_{I}\right]}}^{t}\left[C\right]\left[{B}_{J}\right]\left\{{v}_{n}^{J}\right\}d{v}^{e}$

- 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:

$\left\{{f}_{I}^{\mathrm{int}}\right\}=\left\{{\left({f}_{I}^{\mathrm{int}}\right)}^{0}\right\}+\left\{{\left({f}_{I}^{\mathrm{int}}\right)}^{H}\right\}$

With $\left\{{\left({f}_{I}^{\mathrm{int}}\right)}^{0}\right\}$ the constant part being computed with one-point quadrature, and

$\left\{{\left({f}_{I}^{\mathrm{int}}\right)}^{H}\right\}={{\displaystyle {\int}_{ve}\left[{\left({B}_{I}\right)}^{H}\right]}}^{t}\left[C\right]\left[{\left({B}_{J}\right)}^{H}\right]\left\{{v}_{n}^{J}\right\}d{v}^{e}$

These can be evaluated explicitly.

The rate-of-deformations will be written explicitly.

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

$\dot{\sigma}=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 $\dot{q}$:

${\left(q\right)}_{n+1}={\left(q\right)}_{n}+{\dot{q}}_{n+\frac{1}{2}}\text{\Delta}t$

Noting that if $\varphi {,}_{\alpha}$ 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:

$\left\{{f}_{I}^{\mathrm{int}}\right\}={\displaystyle {\int}_{ve}{\left[{B}_{I}\right]}^{t}\left\{\sigma \right\}}d{v}^{e}$

The constitutive relation is written by either a tangent form: $\dot{\sigma}={C}^{t}:D$, or a projection form: $\sigma =P\left({\sigma}^{e},\mathrm{...}\right)$

Where, ${C}^{t}\left({\sigma}^{},\mathrm{...}\right)$ is the history-dependent tangent tensor; $\left\{{\sigma}_{n+1}^{e}\right\}=\left\{{\sigma}_{n}\right\}+\left[C\right]\left\{{D}_{n+1/2}\right\}\text{\Delta}t$ is the trial stress, and the function $$P$$ consists of projecting the trial stress on the updated yield surface.

The decomposed form of the last equation is written:

$\left\{{f}_{I}^{\mathrm{int}}\right\}=\left\{{\left({f}_{I}^{\mathrm{int}}\right)}^{0}\right\}+\left\{{\left({f}_{I}^{\mathrm{int}}\right)}^{H}\right\}$

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

- Keep the same orthogonalities as in the elastic case, and
- Use the unidimensional tangent modulus ${E}_{t}\left(\zeta \right)$ to evaluate the non-constant rate-stress, that
is,
(19) ${\left\{\dot{\sigma}\right\}}^{H}=\left({E}_{t}\left(\zeta \right)/E\right)\left[C\right]{\left\{D\right\}}^{H}$ Where,- $$E$$
- Young's modulus
- $\left[C\right]$
- Matrix of elastic moduli

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

Where, ${E}_{t}\left(\zeta \right)$ is obtained by the constant stress incremental computation along the thickness and ${E}_{t}\left(\zeta \right)=E$ in the elastic zone, and ${\lambda}_{m}={\overline{E}}_{t}/E;{\overline{E}}_{t}$ is the average value of ${E}_{t}\left(\zeta \right)$ and ${\lambda}_{b}=\mathrm{min}{E}_{t}\left(\zeta \right)/E$.

### For the QPH shell, $${\lambda}_{m}$$= 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.

- Take the elastic matrix $\left[C\right]$ directly:
$\left\{\dot{\sigma}\right\}={\left\{\dot{\sigma}\right\}}^{0}+\left[C\right]\left\{{\left(\dot{\epsilon}-{\dot{\epsilon}}^{P}\right)}^{H}\right\}={\left\{\dot{\sigma}\right\}}^{0}+\left[C\right]\left\{{\left(\dot{\epsilon}\right)}^{H}\right\}\to \left\{{\left({\dot{\epsilon}}^{P}\right)}^{H}\right\}=0$

which means that the plastic rate-of-deformation $\left\{{\dot{\epsilon}}^{P}\right\}$ is constant within the element; or

- Take the tangent matrix [$${C}^{t}$$] ($\xi =\eta =0\text{\hspace{1em}},\zeta =$ 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 $\left[C\right]$ 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.