Interpolation Functions

For the FE solution to converge to the exact one as the mesh is refined, the assumed displacement function must satisfy these conditions.

Requirements

1. Completeness. The function must be able to represent rigid-body motion (constant terms) and constant strain states. For example, a bar needs $u = a_1 + a_2x$: $a_1$ is rigid-body translation and $a_2$ is constant strain. As elements shrink, the strain in each becomes nearly constant, so this is essential.
2. Compatibility (conformity). Displacements must be continuous across element boundaries so that no gaps or overlaps develop. For second-order problems (bars, plane elasticity) this means $C^0$ continuity; for beams the slope must also be continuous ($C^1$).
3. Rigid-body modes. An element should have zero strain energy under rigid-body movement.
4. Geometric isotropy (invariance). The function should not prefer any direction, so the result should not change if the local axes are rotated. This needs complete polynomial terms of a given order.

Elements that satisfy all the conditions converge monotonically; those that satisfy only completeness may still converge (non-conforming elements pass the patch test).

Pascal's triangle

Order 0:           1
Order 1:         x   y
Order 2:       x2  xy  y2
Order 3:     x3  x2y xy2  y3
Order 4:   x4  x3y x2y2 xy3 y4

• Terms are taken row by row so that the polynomial is complete up to a given degree and has geometric isotropy.
• Number of terms must equal the number of nodal degrees of freedom: 3 terms (1, x, y) for the 3-node triangle; 4 terms (1, x, y, xy) for the 4-node rectangle, where $xy$ is chosen symmetrically.
• Terms must be selected symmetrically about the diagonal of the triangle to maintain isotropy.

Displacement function

Two nodal values $u_1$ (at $x = 0$) and $u_2$ (at $x = L$) allow a linear polynomial:

$$u(x) = a_1 + a_2x$$

   u1                u2
   o-----------------o
  x=0               x=L

Apply the nodal conditions:

$$u(0) = a_1 = u_1,\qquad u(L) = a_1 + a_2L = u_2 \Rightarrow a_2 = \frac{u_2 - u_1}{L}$$

Substitute back:

$$u = u_1 + \frac{u_2 - u_1}{L}x = \left(1 - \frac{x}{L}\right)u_1 + \frac{x}{L}\,u_2$$

So $u = N_1u_1 + N_2u_2$ with

$$N_1 = 1 - \frac{x}{L},\qquad N_2 = \frac{x}{L}$$

Properties

• $N_i = 1$ at node $i$ and $0$ at the other node (Kronecker delta property).
• $N_1 + N_2 = 1$ everywhere (partition of unity); this allows rigid-body motion.
• Each $N_i$ is linear, so the displacement is continuous between elements ($C^0$), but the strain is constant per element.
• In natural coordinate $\xi = 2x/L - 1$: $N_1 = \dfrac{1 - \xi}{2}$ and $N_2 = \dfrac{1 + \xi}{2}$.

Strain-displacement matrix

$$\varepsilon = \frac{du}{dx} = \frac{dN_1}{dx}u_1 + \frac{dN_2}{dx}u_2 = \frac{1}{L}\left[-1\ \ 1\right]\begin{Bmatrix} u_1 \\ u_2 \end{Bmatrix}$$ $$[B] = \frac{1}{L}\left[-1\ \ 1\right]$$

Mesh and element areas

Area varies as $A(x) = 600 - x$ (mm$^2$, $x$ in mm from the fixed end).

 |#### 600 -> 400 -> 200
 |###################>---> 15 kN
 1        2          3

Element stiffness

$$\begin{aligned} k_1 &= \frac{A_1E}{L} = \frac{500(200\times10^3)}{200} = 5\times10^5\ \text{N/mm} \\ k_2 &= \frac{300(200\times10^3)}{200} = 3\times10^5\ \text{N/mm} \end{aligned}$$

Global equations with $u_1 = 0$

$$\begin{bmatrix} 8\times10^5 & -3\times10^5 \\ -3\times10^5 & 3\times10^5 \end{bmatrix}\begin{Bmatrix} u_2 \\ u_3 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 15\,000 \end{Bmatrix}$$

Both elements carry the same force, 15 kN, so:

$$u_2 = \frac{15\,000}{5\times10^5} = 0.03\ \text{mm},\qquad u_3 = u_2 + \frac{15\,000}{3\times10^5} = 0.03 + 0.05 = 0.08\ \text{mm}$$

Stresses

$$\sigma_1 = \frac{15\,000}{500} = 30\ \text{MPa},\qquad \sigma_2 = \frac{15\,000}{300} = 50\ \text{MPa}$$

Exact solution

$$u_{tip} = \int_0^{400}\frac{P\,dx}{E(600 - x)} = \frac{15\,000}{200\times10^3}\ln\frac{600}{200} = 0.075(1.0986) = 0.08240\ \text{mm}$$

The FE tip value is 0.0800 mm, which is 2.9 % below the exact value. The stress at the fixed end, $15\,000/600 = 25$ MPa, and at the tip, $75$ MPa, differ from the element values (30 and 50 MPa); the error reduces with more elements.

Answer: $u_2 = 0.030$ mm, $u_3 = 0.080$ mm; $\sigma_1 = 30$ MPa, $\sigma_2 = 50$ MPa; exact tip displacement 0.0824 mm.

Consider a triangle with nodes 1, 2, 3 at $(x_1, y_1)$, $(x_2, y_2)$, $(x_3, y_3)$ numbered anticlockwise. Each node has two displacements $(u_i, v_i)$, so the element has 6 degrees of freedom.

         3 (x3,y3)
        / \
       /   \
      /     \
     1-------2

Displacement field

Choose a complete linear polynomial (from Pascal's triangle):

$$u = \alpha_1 + \alpha_2x + \alpha_3y,\qquad v = \alpha_4 + \alpha_5x + \alpha_6y$$

Apply the nodal conditions for $u$:

$$\begin{Bmatrix} u_1 \\ u_2 \\ u_3 \end{Bmatrix} = \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{bmatrix}\begin{Bmatrix} \alpha_1 \\ \alpha_2 \\ \alpha_3 \end{Bmatrix}$$

Solving by Cramer's rule, where the determinant of the matrix equals $2A$ ($A$ = area of triangle):

$$u = N_1u_1 + N_2u_2 + N_3u_3,\qquad N_i = \frac{1}{2A}\left(a_i + b_ix + c_iy\right)$$

with

$$\begin{aligned} a_1 &= x_2y_3 - x_3y_2, & b_1 &= y_2 - y_3, & c_1 &= x_3 - x_2 \\ a_2 &= x_3y_1 - x_1y_3, & b_2 &= y_3 - y_1, & c_2 &= x_1 - x_3 \\ a_3 &= x_1y_2 - x_2y_1, & b_3 &= y_1 - y_2, & c_3 &= x_2 - x_1 \end{aligned}$$ $$2A = x_1(y_2 - y_3) + x_2(y_3 - y_1) + x_3(y_1 - y_2)$$

The same $N_i$ apply to $v$, so

$$\begin{Bmatrix} u \\ v \end{Bmatrix} = \begin{bmatrix} N_1 & 0 & N_2 & 0 & N_3 & 0 \\ 0 & N_1 & 0 & N_2 & 0 & N_3 \end{bmatrix}\{d\}$$

Properties

• $N_i = 1$ at node $i$ and $0$ at the other two nodes.
• $N_1 + N_2 + N_3 = 1$ at every point (and $\sum N_ix_i = x$, $\sum N_iy_i = y$).
• $N_i$ varies linearly along each edge and depends only on the two nodes of that edge, so displacements are continuous between adjacent elements.
• $N_i$ equal the area coordinates $L_i = A_i/A$.

Why constant strain

The strains are

$$\varepsilon_x = \frac{\partial u}{\partial x} = \frac{1}{2A}\sum b_iu_i,\quad \varepsilon_y = \frac{1}{2A}\sum c_iv_i,\quad \gamma_{xy} = \frac{1}{2A}\sum (c_iu_i + b_iv_i)$$

$b_i$ and $c_i$ depend only on the nodal coordinates, not on $x$ or $y$. So the strains (and stresses) are the same everywhere inside the element.

(a) Area

$$2A = x_1(y_2 - y_3) + x_2(y_3 - y_1) + x_3(y_1 - y_2) = 10(20 - 50) + 60(50 - 10) + 30(10 - 20)$$ $$2A = -300 + 2400 - 300 = 1800 \Rightarrow A = 900\ \text{mm}^2$$

The positive value shows the nodes are numbered anticlockwise.

(b) Shape functions

Constants:

At P (30, 25), $N_i = \dfrac{a_i + b_i(30) + c_i(25)}{1800}$:

$$\begin{aligned} N_1 &= \frac{2400 - 900 - 750}{1800} = \frac{750}{1800} = 0.4167 \\ N_2 &= \frac{-200 + 1200 - 500}{1800} = \frac{500}{1800} = 0.2778 \\ N_3 &= \frac{-400 - 300 + 1250}{1800} = \frac{550}{1800} = 0.3056 \end{aligned}$$

Sum $= 0.4167 + 0.2778 + 0.3056 = 1.0000$, as required.

Check on location: $x = \sum N_ix_i = 0.4167(10) + 0.2778(60) + 0.3056(30) = 30.0$ and $y = 0.4167(10) + 0.2778(20) + 0.3056(50) = 25.0$.

(c) Temperature at P

$$T_P = 0.4167(100) + 0.2778(60) + 0.3056(80) = 41.67 + 16.67 + 24.44 = 82.78\ ^\circ\text{C}$$

Answer: $A = 900\ \text{mm}^2$; $N = (0.417,\ 0.278,\ 0.306)$; $T_P = 82.8\ ^\circ$C.

   4 (0,b)  o-----------o  3 (a,b)
            |           |
            |           |  b
            |           |
   1 (0,0)  o-----------o  2 (a,0)
                  a

Displacement function

Four nodes per scalar field need four terms. Taking terms symmetrically from Pascal's triangle gives the bilinear polynomial

$$u(x, y) = \alpha_1 + \alpha_2x + \alpha_3y + \alpha_4xy$$

(The $xy$ term is chosen; $x^2$ or $y^2$ would make the function unsymmetrical.)

Apply the four nodal conditions $u(0,0) = u_1$, $u(a,0) = u_2$, $u(a,b) = u_3$, $u(0,b) = u_4$:

$$\begin{aligned} \alpha_1 &= u_1 \\ \alpha_2 &= \frac{u_2 - u_1}{a} \\ \alpha_3 &= \frac{u_4 - u_1}{b} \\ \alpha_4 &= \frac{u_1 - u_2 + u_3 - u_4}{ab} \end{aligned}$$

Substituting and collecting the coefficients of $u_i$ gives $u = \sum N_iu_i$ with

$$\begin{aligned} N_1 &= \left(1 - \frac{x}{a}\right)\left(1 - \frac{y}{b}\right), & N_2 &= \frac{x}{a}\left(1 - \frac{y}{b}\right) \\ N_3 &= \frac{xy}{ab}, & N_4 &= \left(1 - \frac{x}{a}\right)\frac{y}{b} \end{aligned}$$

The same functions interpolate $v$.

Properties verified

• At node 3 ($x = a$, $y = b$): $N_3 = 1$, $N_1 = N_2 = N_4 = 0$. Similarly for other nodes.
• Sum: $N_1 + N_2 + N_3 + N_4 = \left(1 - \frac{x}{a}\right)\left(1 - \frac{y}{b}\right) + \frac{x}{a}\left(1 - \frac{y}{b}\right) + \frac{xy}{ab} + \left(1 - \frac{x}{a}\right)\frac{y}{b} = 1$.
• Along an edge, e.g. $y = 0$, only $N_1$ and $N_2$ are non-zero and linear in $x$, so the displacement is continuous between neighbouring elements.

Strain variation

$$\varepsilon_x = \frac{\partial u}{\partial x} = \alpha_2 + \alpha_4y,\qquad \varepsilon_y = \frac{\partial v}{\partial y} = \beta_3 + \beta_4x$$

$\varepsilon_x$ varies linearly with $y$ and $\varepsilon_y$ varies linearly with $x$, so the element is better than CST at representing bending. A drawback is that it can only be used for rectangles with sides parallel to the axes; the quadrilateral element (isoparametric) removes this limit.

Governing equation and functional

Steady 1D conduction with heat generation:

$$kA\frac{d^2T}{dx^2} + QA = 0$$

The equivalent functional (analogous to the total potential energy) is

$$\Pi = \int_0^L \left[\frac{kA}{2}\left(\frac{dT}{dx}\right)^2 - QA\,T\right]dx$$

The solution of the differential equation makes $\Pi$ stationary (this can be shown by taking $\delta\Pi = 0$ and integrating by parts, which gives back the governing equation).

Element approximation

For a linear element with nodal temperatures $T_1$, $T_2$:

$$T = N_1T_1 + N_2T_2,\quad N_1 = 1 - \frac{x}{L},\quad N_2 = \frac{x}{L}$$ $$\frac{dT}{dx} = \frac{1}{L}[-1\ \ 1]\begin{Bmatrix} T_1 \\ T_2 \end{Bmatrix} = [B]\{T\}$$

Minimising the functional

Substitute into $\Pi$:

$$\Pi = \frac{1}{2}\{T\}^T\left[\int_0^L kA\,[B]^T[B]\,dx\right]\{T\} - \{T\}^T\int_0^L QA\,[N]^T dx$$

Setting $\dfrac{\partial\Pi}{\partial T_i} = 0$ gives $[k]\{T\} = \{f\}$, where

$$[k] = \int_0^L kA\,[B]^T[B]\,dx = \frac{kA}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$$ $$\{f\} = \int_0^L QA\begin{Bmatrix} N_1 \\ N_2 \end{Bmatrix}dx = \frac{QAL}{2}\begin{Bmatrix} 1 \\ 1 \end{Bmatrix}$$

Element equation

$$\frac{kA}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\begin{Bmatrix} T_1 \\ T_2 \end{Bmatrix} = \frac{QAL}{2}\begin{Bmatrix} 1 \\ 1 \end{Bmatrix} + \begin{Bmatrix} q_1 \\ q_2 \end{Bmatrix}$$

where $q_1$, $q_2$ are nodal heat flows at the ends (natural boundary conditions). The generated heat is shared equally between the two nodes. Convection at an end with coefficient $h$ and ambient temperature $T_\infty$ adds $hA$ to the stiffness term at that node and $hAT_\infty$ to the load.

The assembly, boundary conditions ($T$ specified at a node) and solution follow the same steps as for the bar problem.