One-Point Quadrature Shell Element

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 ¹ 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 ² which free either drilling rotation or rigid motion should be used. This full projection is only used for QEPH element.

Constant part

with:

${\text{Δ}}_I=\left[t_I-\left(t_I{x}^I\right)b_{xI}-\left(t_I{y}^I\right)b_{yI}\right];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\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;\varphi =\xi \eta ;$

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

with

$<b_{xI}^c>=4z_{\gamma }/A^2<b_{y2};b_{y1};b_{y4};b_{y3}>$

$<b_{yI}^c>=4z_{\gamma }/A^2<b_{x4};b_{x3};b_{x2};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}\left({\overline{f}}_{ij}=\left(f_i+f_j\right)/2\right),$ 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 ⁴, ⁶though this part has been obtained in a very different manner. Letellier has given the same result in his thesis ⁵, 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$

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).

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.

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$

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{Δ}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{Δ}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.

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=\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.

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 ⁶, and it has been shown that this choice has little effect on the result of the "Cylindrical panel test". In ⁷, 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.