Higher order Elements

Both are higher-order rectangular (or quadrilateral) elements built on a $\xi$-$\eta$ parent square ($-1 \le \xi, \eta \le 1$).

Comparison

 9-node Lagrange        8-node serendipity
 4--7--3                4--7--3
 |     |                |     |
 8  9  6                8     6
 |     |                |     |
 1--5--2                1--5--2

Polynomial terms (Pascal's triangle)

• 9-node Lagrange: $1,\ \xi,\ \eta,\ \xi^2,\ \xi\eta,\ \eta^2,\ \xi^2\eta,\ \xi\eta^2,\ \xi^2\eta^2$.
• 8-node serendipity: the same list without $\xi^2\eta^2$.

The 9-node element is a complete biquadratic. The 8-node element is a complete quadratic with two additional cubic terms, $\xi^2\eta$ and $\xi\eta^2$, which keep the polynomial symmetric.

A serendipity element is preferred in most practice because the interior node adds degrees of freedom but little accuracy.

Each shape function must equal 1 at its own node and 0 at the other seven nodes, and the sum must be 1.

Mid-side node 5 (0, -1)

$N_5$ must vanish at the other two nodes of the edge $\eta = -1$ ($\xi = \pm1$) and on the edge $\eta = +1$, and at nodes 6 and 8 ($\xi = \pm1$). Use a quadratic in $\xi$ times a linear function in $\eta$:

$$N_5 = c\,(1 - \xi^2)(1 - \eta)$$

At node 5, $N_5 = c(1)(2) = 1$, so $c = \dfrac{1}{2}$:

$$N_5 = \tfrac{1}{2}(1 - \xi^2)(1 - \eta)$$

Corner node 1 (-1, -1)

Start with the bilinear function $\tfrac{1}{4}(1 - \xi)(1 - \eta)$. It is 1 at node 1 and 0 at nodes 2, 3, 4, but it is non-zero at mid-side nodes 5 and 8 (value $\tfrac{1}{2}$ each). Subtract $\tfrac{1}{2}N_5$ and $\tfrac{1}{2}N_8$:

$$N_1 = \tfrac{1}{4}(1 - \xi)(1 - \eta) - \tfrac{1}{2}N_5 - \tfrac{1}{2}N_8$$

With $N_8 = \tfrac{1}{2}(1 - \xi)(1 - \eta^2)$, simplifying gives

$$N_1 = \tfrac{1}{4}(1 - \xi)(1 - \eta)(-\xi - \eta - 1)$$

Check: at node 1 ($\xi = \eta = -1$): $\tfrac{1}{4}(2)(2)(1 + 1 - 1) = 1$. At node 2 ($\xi = 1$): factor $(1 - \xi) = 0$. At node 5 ($\xi = 0$, $\eta = -1$): factor $(1 - \eta)(-\xi - \eta - 1) = 2(0) = 0$. At node 8 similarly zero.

General formulae ($\xi_i$, $\eta_i$ = nodal values)

• Corner nodes ($i = 1$ to 4):

$$N_i = \tfrac{1}{4}(1 + \xi\xi_i)(1 + \eta\eta_i)(\xi\xi_i + \eta\eta_i - 1)$$

• Mid-side nodes with $\xi_i = 0$ (nodes 5, 7):

$$N_i = \tfrac{1}{2}(1 - \xi^2)(1 + \eta\eta_i)$$

• Mid-side nodes with $\eta_i = 0$ (nodes 6, 8):

$$N_i = \tfrac{1}{2}(1 + \xi\xi_i)(1 - \eta^2)$$

The edge shape functions are quadratic along each edge, so adjacent elements share three nodes per edge and compatibility is maintained.

Parametric mapping and isoparametric element

An arbitrary quadrilateral in the $x$-$y$ plane is mapped from a square parent element with natural coordinates $-1 \le \xi, \eta \le 1$. The same shape functions are used for the geometry and for the displacement:

$$x = \sum N_ix_i,\quad y = \sum N_iy_i,\quad u = \sum N_iu_i,\quad v = \sum N_iv_i$$

Equal order for geometry and displacement is the isoparametric concept. Derivatives need the chain rule, so the Jacobian matrix is required:

$$\begin{Bmatrix} \partial N/\partial\xi \\ \partial N/\partial\eta \end{Bmatrix} = [J]\begin{Bmatrix} \partial N/\partial x \\ \partial N/\partial y \end{Bmatrix},\qquad [J] = \begin{bmatrix} \partial x/\partial\xi & \partial y/\partial\xi \\ \partial x/\partial\eta & \partial y/\partial\eta \end{bmatrix},\quad dA = \det[J]\,d\xi\,d\eta$$

Shape functions at $\xi = \eta = 0.5$

$N_i = \tfrac{1}{4}(1 + \xi\xi_i)(1 + \eta\eta_i)$:

$$N_1 = \tfrac{1}{4}(0.5)(0.5) = 0.0625,\ N_2 = 0.1875,\ N_3 = 0.5625,\ N_4 = 0.1875$$

(a) Global coordinates

$$\begin{aligned} x &= 0.0625(0) + 0.1875(4) + 0.5625(5) + 0.1875(1) = 3.75\ \text{mm} \\ y &= 0 + 0 + 0.5625(3) + 0.1875(3) = 2.25\ \text{mm} \end{aligned}$$

(b) Jacobian

Derivatives with respect to the natural coordinates:

$$\frac{\partial N_i}{\partial\xi} = \tfrac{1}{4}\xi_i(1 + \eta\eta_i),\qquad \frac{\partial N_i}{\partial\eta} = \tfrac{1}{4}\eta_i(1 + \xi\xi_i)$$ $$\begin{aligned} J_{11} &= \sum \frac{\partial N_i}{\partial\xi}x_i = 0.125(4) + 0.375(5) - 0.375(1) = 2.0 \\ J_{12} &= \sum \frac{\partial N_i}{\partial\xi}y_i = 0.375(3) - 0.375(3) = 0 \\ J_{21} &= \sum \frac{\partial N_i}{\partial\eta}x_i = -0.375(4) + 0.375(5) + 0.125(1) = 0.5 \\ J_{22} &= \sum \frac{\partial N_i}{\partial\eta}y_i = 0.375(3) + 0.125(3) = 1.5 \end{aligned}$$ $$[J] = \begin{bmatrix} 2.0 & 0 \\ 0.5 & 1.5 \end{bmatrix},\qquad \det[J] = 3.0$$

(The element is a parallelogram of area $4 \times 3 = 12$ mm$^2$ = $4\det J$, because the parent area is 4, so $J$ is constant.)

(c) Global derivatives

$$[J]^{-1} = \frac{1}{3}\begin{bmatrix} 1.5 & 0 \\ -0.5 & 2.0 \end{bmatrix}$$

Each column sums to zero, as it must because $\sum N_i = 1$.

Answer: $(x, y) = (3.75,\ 2.25)$ mm; $[J] = \begin{bmatrix} 2 & 0 \\ 0.5 & 1.5 \end{bmatrix}$, $\det J = 3.0$; derivatives as tabulated.

Gauss quadrature

In isoparametric elements, the stiffness integral $\int [B]^T[D][B]\det J\,d\xi\,d\eta$ is rational in $\xi$ and $\eta$, so it is evaluated numerically over $-1 \le \xi \le 1$:

$$\int_{-1}^{1} f(\xi)\,d\xi \approx \sum_{i=1}^{n} w_i\,f(\xi_i)$$

An $n$-point rule has sampling points $\xi_i$ and weights $w_i$ chosen to integrate exactly a polynomial of degree up to $2n - 1$.

For two dimensions, the rule is applied in both directions, $\sum\sum w_iw_jf(\xi_i, \eta_j)$ ($2\times2$ gives 4 points for a 4-node element).

Numerical example

Change of limits: $x = 1 + \xi$, so $dx = d\xi$ and $f(\xi) = (1 + \xi)^4 + 3(1 + \xi)$.

Exact value:

$$I = \left[\frac{x^5}{5} + \frac{3x^2}{2}\right]_0^2 = \frac{32}{5} + 6 = 12.4$$

2-point rule: $\xi = \pm0.57735$

$$I_2 = 1.29986 + 10.92236 = 12.222$$

The error is $-1.4\%$; the rule is not exact because $x^4$ is of degree 4 > 3.

3-point rule: $\xi = -0.77460,\ 0,\ +0.77460$

$$I_3 = \tfrac{5}{9}(0.67879 + 15.24121) + \tfrac{8}{9}(4.0) = 8.8444 + 3.5556 = 12.400$$

This equals the exact value, since the polynomial is of degree 4 ≤ 5.

Answer: 2-point: $I = 12.222$ (error 1.4 %); 3-point: $I = 12.400$ (exact).