Finite element method

Approach

A cubic triangle has 10 nodes (3 corners, 6 mid-third side nodes, 1 centroid) and a complete cubic polynomial has 10 terms ($1,x,y,x^2,xy,y^2,x^3,x^2y,xy^2,y^3$). It is easiest to use the area (triangular) coordinates $L_1,L_2,L_3$ with $L_1+L_2+L_3=1$, where $L_i=\dfrac{a_i+b_ix+c_iy}{2A}$. The displacement is

$$u=\sum_{i=1}^{10}N_i\,u_i$$

and each $N_i$ is a cubic that equals 1 at node $i$ and 0 at the other nine nodes.

Node coordinates in area coordinates ($L_1,L_2,L_3$): nodes 1, 2, 3 = $(1,0,0),(0,1,0),(0,0,1)$; node 4 = $(\tfrac23,\tfrac13,0)$, 5 = $(\tfrac13,\tfrac23,0)$; 6 = $(0,\tfrac23,\tfrac13)$, 7 = $(0,\tfrac13,\tfrac23)$; 8 = $(\tfrac13,0,\tfrac23)$, 9 = $(\tfrac23,0,\tfrac13)$; 10 = $(\tfrac13,\tfrac13,\tfrac13)$. (Mid-side nodes are numbered anticlockwise, the lower number nearer the lower-numbered corner.)

              3
             / \
          8 /   \ 7
           /     \
        9 /   10  \ 6
         /         \
        1---4---5---2

Corner node 1

$N_1$ must vanish on the lines $L_1=0$ (nodes 2, 3, 6, 7), $L_1=\tfrac13$ (nodes 5, 8, and the centre 10) and $L_1=\tfrac23$ (nodes 4, 9). So

$$N_1=C\,L_1\left(L_1-\tfrac13\right)\left(L_1-\tfrac23\right)$$

At node 1, $L_1=1$: $C\times1\times\tfrac23\times\tfrac13=1\Rightarrow C=\tfrac92$:

$$N_1=\tfrac12L_1(3L_1-1)(3L_1-2)$$

Similarly $N_2=\tfrac12L_2(3L_2-1)(3L_2-2)$ and $N_3=\tfrac12L_3(3L_3-1)(3L_3-2)$.

Side node 4 (on side 1-2, $L_1=\tfrac23$, $L_2=\tfrac13$)

$N_4$ vanishes on $L_1=0$, $L_2=0$ and $L_1=\tfrac13$, so $N_4=C\,L_1L_2(L_1-\tfrac13)$. At node 4: $C\times\tfrac23\times\tfrac13\times\tfrac13=1\Rightarrow C=\tfrac{27}{2}$:

$$N_4=\tfrac92L_1L_2(3L_1-1)$$

By the same reasoning (node 5 is the mirror image, vanishing at $L_2=\tfrac13$):

$$N_5=\tfrac92L_1L_2(3L_2-1)$$

and for the other sides, cyclically:

$$\begin{aligned} N_6&=\tfrac92L_2L_3(3L_2-1),& N_7&=\tfrac92L_2L_3(3L_3-1)\\ N_8&=\tfrac92L_3L_1(3L_3-1),& N_9&=\tfrac92L_3L_1(3L_1-1) \end{aligned}$$

Central node 10

$N_{10}$ must vanish on all three sides: $N_{10}=C\,L_1L_2L_3$; at the centre $C/27=1$:

$$N_{10}=27\,L_1L_2L_3$$

Checks

• $N_i=1$ at node $i$ and $0$ at all other nodes (checked for all ten nodes).
• $\sum_{i=1}^{10}N_i=1$ everywhere, so rigid-body motion is represented.
• Along a side the displacement is a cubic defined by the four nodes on that side, so adjacent elements are compatible.

The element is a complete cubic, so it gives a quadratic variation of strain across the element.

Natural coordinates and nodes

Use $\xi=x/a$, $\eta=y/b$ so the element is the square $-1\le\xi,\eta\le1$ (for half-sides $a,b$; if the figure has unit half-sides, $\xi=x$, $\eta=y$). Nodes:

   4-----7-----3
   |           |
   8     +     6      eta ^
   |           |          |
   1-----5-----2          +--> xi

Assumed displacement

Eight nodal values need an 8-term polynomial (the serendipity family: complete quadratic plus two cubic terms):

$$u=\alpha_1+\alpha_2\xi+\alpha_3\eta+\alpha_4\xi^2+\alpha_5\xi\eta+\alpha_6\eta^2+\alpha_7\xi^2\eta+\alpha_8\xi\eta^2$$

Writing it at the 8 nodes gives $\{u\}=[C]\{\alpha\}$ and $N=[1\ \xi\ \eta\ \dots][C]^{-1}$. The same result is obtained more quickly by building each $N_i$ as a product of straight-line factors that vanish at all the other nodes, as below.

Mid-side node 5 $(0,-1)$

$N_5$ must vanish on the sides $\xi=\pm1$ and $\eta=1$, and have a parabola along the bottom side: $N_5=C(1-\xi^2)(1-\eta)$. At node 5, $C(1)(2)=1\Rightarrow C=\tfrac12$. Likewise for nodes 6, 7, 8:

$$\begin{aligned} N_5&=\tfrac12(1-\xi^2)(1-\eta)\\ N_6&=\tfrac12(1+\xi)(1-\eta^2)\\ N_7&=\tfrac12(1-\xi^2)(1+\eta)\\ N_8&=\tfrac12(1-\xi)(1-\eta^2) \end{aligned}$$

Corner node 1 $(-1,-1)$

$N_1$ must vanish on the sides $\xi=1$ and $\eta=1$ (nodes 2, 6, 3, 7, 4) and at the mid-side nodes 5 and 8, which lie on the line $1+\xi+\eta=0$: $N_1=C(1-\xi)(1-\eta)(1+\xi+\eta)$. At node 1: $C(2)(2)(-1)=1\Rightarrow C=-\tfrac14$. Hence

$$N_1=-\tfrac14(1-\xi)(1-\eta)(1+\xi+\eta)$$

and for the other corners the general form is

$$\boxed{N_i=\tfrac14(1+\xi\xi_i)(1+\eta\eta_i)(\xi\xi_i+\eta\eta_i-1)},\quad i=1,2,3,4$$

that is

$$\begin{aligned} N_2&=\tfrac14(1+\xi)(1-\eta)(\xi-\eta-1)\\ N_3&=\tfrac14(1+\xi)(1+\eta)(\xi+\eta-1)\\ N_4&=\tfrac14(1-\xi)(1+\eta)(-\xi+\eta-1) \end{aligned}$$

Checks

At the centre $(0,0)$: corners give $4\times\tfrac14(-1)=-1$, mid-sides give $4\times\tfrac12=2$; the sum is $1$. In general $\sum N_i=1$ and each $N_i=1$ at its own node and $0$ at the other seven. Each side carries a quadratic variation through three nodes, so adjacent elements are compatible.

Data: $EI=2\times10^5\ \text{MPa}\times5\times10^6\ \text{mm}^4=10^{12}$ N mm$^2=1000$ kN m$^2$ (same for both spans). Nodes: 1 = fixed end (x = 0); 2 = intermediate support (x = 7 m) with an applied moment of 25 kN m (taken clockwise); 3 = right roller (x = 12 m). Span 1: $L=7$ m with 25 kN at 3 m from the fixed end; span 2: $L=5$ m with 15 kN/m UDL. Units: kN, m.

Element stiffness

$$[k]=\frac{EI}{L^3}\begin{bmatrix}12&6L&-12&6L\\6L&4L^2&-6L&2L^2\\-12&-6L&12&-6L\\6L&2L^2&-6L&4L^2\end{bmatrix},\quad \{d\}=\{v_i,\theta_i,v_j,\theta_j\}^T$$

(Upward forces and displacements, and counter-clockwise moments and rotations, are positive.)

Equivalent nodal loads

Span 1 ($a=3$, $b=4$, $L=7$): $\dfrac{Pb^2(3a+b)}{L^3}=\dfrac{25\times16\times13}{343}=15.16$ kN, $\dfrac{Pab^2}{L^2}=24.49$ kN m, $\dfrac{Pa^2(a+3b)}{L^3}=9.84$ kN, $\dfrac{Pa^2b}{L^2}=18.37$ kN m. Span 2: $wL/2=37.5$ kN, $wL^2/12=31.25$ kN m. Equivalent loads are the negatives of the fixed-end forces. The applied clockwise moment gives $-25$ kN m at $\theta_2$ (counter-clockwise positive).

Assembled equation

Boundary conditions: $v_1=\theta_1=0$, $v_2=0$, $v_3=0$. Free DOFs: $\theta_2,\theta_3$.

$$\begin{bmatrix}1371.429 & 400.000 \\ 400.000 & 800.000\end{bmatrix}\begin{Bmatrix}\theta_2\\\theta_3\end{Bmatrix}=\begin{Bmatrix}-37.883 \\ 31.250\end{Bmatrix}$$ $$\theta_2=-0.04568\ \text{rad},\qquad \theta_3=0.06190\ \text{rad}$$

Support reactions $\{f\}=[k]\{d\}-\{F_{eq}\}$

• Fixed end: $R_1=9.567$ kN (up), $M_1=11.439$ kN m (counter-clockwise)
• Intermediate support: $R_2=56.827$ kN (up)
• Right roller: $R_3=33.606$ kN (up)
• Check: $\Sigma R=100.000$ kN $=25+15\times5=100$ kN.

Deflection at mid span

The displacement within a span is the Hermite interpolation of the nodal values plus the fixed-end beam deflection due to the load. Evaluating (checked by refining the mesh with nodes at the mid-points):

• Mid of span 1 (x = 3.5 m): $v=-2.22$ mm (downward)
• Mid of span 2 (x = 9.5 m): $v=-91.65$ mm (downward)

The value is large because $EI=1000$ kN m$^2$ is small for these spans.

Answer: $R_1=9.57$ kN, $M_1=11.44$ kN m, $R_2=56.83$ kN, $R_3=33.61$ kN; mid-span deflections $2.22$ mm (span 1) and $91.65$ mm (span 2), both downward.

The beam is solved by the direct stiffness method with three beam elements.

Data and idealisation

$EI=(20\times10^{5}\ \text{N/mm}^2)(5\times10^{6}\ \text{mm}^4)=1\times10^{13}\ \text{N·mm}^2=10\ \text{kN·m}^2$.

        30 kN     40 kNm        20 kN/m
          |       (clockwise)  ||||||||||||
 ##=======o=========o==============o
  A       2         B              C
  fixed   |<- 5 m ->|    |<- 6 m ->|
 |<- 3 m->|

• Nodes: 1 = A (x = 0, fixed), 2 (x = 3 m, point load), 3 = B (x = 8 m, support, $M=40$ kNm taken clockwise), 4 = C (x = 14 m, roller).
• Elements: 1 (3 m), 2 (5 m), 3 (6 m, carries the UDL from B to C).
• Dofs per node: (deflection $v$ up, rotation $\theta$ anticlockwise). Known: $v_1=\theta_1=0$, $v_3=0$, $v_4=0$. Unknown: $v_2,\theta_2,\theta_3,\theta_4$.

Element stiffness (each element)

$$[k]=\frac{EI}{L^3}\begin{bmatrix}12&6L&-12&6L\\6L&4L^2&-6L&2L^2\\-12&-6L&12&-6L\\6L&2L^2&-6L&4L^2\end{bmatrix}$$

Equivalent nodal loads

Element 3: $wL/2=60$ kN down at each end, $wL^2/12=60$ kNm (clockwise at B, anticlockwise at C). Loads on the structure for $(v_2,\theta_2,\theta_3,\theta_4)$:

$$\{F\}=\begin{bmatrix}-30\\0\\-40-60\\+60\end{bmatrix}=\begin{bmatrix}-30\\0\\-100\\60\end{bmatrix}$$

Reduced equations (after assembly)

$$\begin{bmatrix} 5.4044 & -4.2667 & 2.4 & 0 \\ -4.2667 & 21.333 & 4 & 0 \\ 2.4 & 4 & 14.667 & 3.3333 \\ 0 & 0 & 3.3333 & 6.6667 \end{bmatrix}\begin{bmatrix}v_2\\\theta_2\\\theta_3\\\theta_4\end{bmatrix}=\begin{bmatrix}-30\\0\\-100\\60\end{bmatrix}$$

Solution:

$$v_2=1.066\ \text{m},\quad \theta_2=2.255\ \text{rad},\quad \theta_3=-10.891\ \text{rad},\quad \theta_4=14.445\ \text{rad}$$

(equivalently $v_2=10.657/EI$, etc.; the very small $EI$ gives large linear-theory values).

Support reactions ($\{R\}=[K]\{d\}-\{F\}$)

$$R_A=10.298\ \text{kN}\uparrow,\quad M_A=7.930\ \text{kNm (anticlockwise)},\quad R_B=85.627\ \text{kN}\uparrow,\quad R_C=54.076\ \text{kN}\uparrow$$

Check: $\Sigma F_y$: $10.298+85.627+54.076=150.000=30+20\times6$. $\Sigma M_A=0$ also checks.

Deflection at mid-span

Using the element shape functions (and adding the UDL term $wx^2(L-x)^2/24EI$ inside the loaded element):

• Mid-span of A-B ($x=4$ m): $v=4.141$ m (upward).
• Mid-span of B-C ($x=11$ m): $v=-25.752$ m (downward). Without the UDL term the cubic interpolation alone would give $-19.002$ m.

Bending moment (sagging +), from the reactions

 BMD, sagging up (1 row = 12 kNm, 1 column = 0.35 m)
   84 |                                         
   72 |                              *****      
   60 |                             *     **    
   48 |                            *        *   
   36 |                           *          *  
   24 |        **                *              
   12 |    ****  **             *             * 
    0 |-***--------*---------------------------*
  -12 |*            **         *                
  -24 |               **                        
  -36 |                 *     *                 
  -48 |                  **                     
  -60 |                    **                   
  -72 |                      *                  
  -84 |                                         

The diagram runs from A at the left edge to C at the right edge: hogging at A, sagging under the 30 kN load, deep hogging just left of B, a jump of 40 kNm at B, then sagging in B-C with the maximum at 3.30 m from B, and zero at C.

Answer: $R_A=10.298$, $R_B=85.627$, $R_C=54.076$ kN; $M_A=7.930$ kNm; mid-span deflections $4.141$ m (AB) and $-25.752$ m (BC).

Types of parametric elements

• Isoparametric: the same shape functions (same number of nodes) are used for the geometry and for the displacement, $x=\sum N_ix_i$, $u=\sum N_iu_i$.
• Superparametric: the geometry is described by a higher-order interpolation (more nodes) than the displacement, e.g. quadratic geometry and linear displacement.
• Subparametric: the geometry is described by a lower-order interpolation than the displacement, e.g. straight sides with a quadratic displacement.

Hermite shape functions for a beam element

A beam element of length $L$ has two nodes with two DOFs each: deflection $v$ and slope $\theta=dv/dx$. Both $v$ and $\theta$ must be continuous between elements ($C^1$ continuity), so a cubic (4 constants) is assumed:

$$v(x)=\alpha_1+\alpha_2x+\alpha_3x^2+\alpha_4x^3$$ $$\theta(x)=\frac{dv}{dx}=\alpha_2+2\alpha_3x+3\alpha_4x^2$$

Apply the four nodal conditions: $v(0)=v_1$, $\theta(0)=\theta_1$, $v(L)=v_2$, $\theta(L)=\theta_2$:

$$\begin{Bmatrix}v_1\\\theta_1\\v_2\\\theta_2\end{Bmatrix}=\begin{bmatrix}1&0&0&0\\0&1&0&0\\1&L&L^2&L^3\\0&1&2L&3L^2\end{bmatrix}\begin{Bmatrix}\alpha_1\\\alpha_2\\\alpha_3\\\alpha_4\end{Bmatrix}$$

Solving: $\alpha_1=v_1$, $\alpha_2=\theta_1$,

$$\alpha_3=\frac{3(v_2-v_1)}{L^2}-\frac{2\theta_1+\theta_2}{L},\qquad \alpha_4=\frac{2(v_1-v_2)}{L^3}+\frac{\theta_1+\theta_2}{L^2}$$

Substituting and collecting the coefficients of $v_1,\theta_1,v_2,\theta_2$ gives $v=[N_1\ N_2\ N_3\ N_4]\{d\}$ with

$$\begin{aligned} N_1&=1-\frac{3x^2}{L^2}+\frac{2x^3}{L^3}\\ N_2&=x-\frac{2x^2}{L}+\frac{x^3}{L^2}\\ N_3&=\frac{3x^2}{L^2}-\frac{2x^3}{L^3}\\ N_4&=-\frac{x^2}{L}+\frac{x^3}{L^2} \end{aligned}$$

In terms of $s=x/L$: $N_1=1-3s^2+2s^3$, $N_2=L(s-2s^2+s^3)$, $N_3=3s^2-2s^3$, $N_4=L(-s^2+s^3)$.

Check: $N_1(0)=1$, $N_1(L)=0$, $N_1'(0)=N_1'(L)=0$; $N_2'(0)=1$, $N_2(0)=N_2(L)=N_2'(L)=0$, and so on for $N_3,N_4$. Also $N_1+N_3=1$ (rigid-body translation is represented).

Assumptions (figure data): A is fixed, C is a roller (simple support), $EI$ is constant. Three beam elements: AB (4 m, 10 kN/m UDL), BD (2 m) and DC (2 m, no load); the 40 kN downward load acts at node D. Positive: $v$ upward, $\theta$ counter-clockwise. Units: kN, m.

Element stiffness

Each element uses

$$[k]=\frac{EI}{L^3}\begin{bmatrix}12&6L&-12&6L\\6L&4L^2&-6L&2L^2\\-12&-6L&12&-6L\\6L&2L^2&-6L&4L^2\end{bmatrix},\quad \{d\}=\{v_i,\theta_i,v_j,\theta_j\}^T$$

(Upward forces and displacements, and counter-clockwise moments and rotations, are positive.) with $L=4$ m (AB), $L=2$ m (BD, DC). Nodes: 1 = A, 2 = B, 3 = D, 4 = C.

Equivalent nodal loads

Element AB: $wL/2=20$ kN, $wL^2/12=13.333$ kN m, so the equivalent loads are $-20$ kN and $-13.333$ kN m at A, and $-20$ kN and $+13.333$ kN m at B. The 40 kN load is a direct nodal load at D.

Assembled equation (free DOFs $v_B,\theta_B,v_D,\theta_D,\theta_C$)

With $v_A=\theta_A=0$ and $v_C=0$:

$$EI\,\begin{bmatrix}1.688 & 1.125 & -1.500 & 1.500 & 0 \\ 1.125 & 3.000 & -1.500 & 1.000 & 0 \\ -1.500 & -1.500 & 3.000 & 0 & 1.500 \\ 1.500 & 1.000 & 0 & 4.000 & 1.000 \\ 0 & 0 & 1.500 & 1.000 & 2.000\end{bmatrix}\begin{Bmatrix}v_B\\\theta_B\\v_D\\\theta_D\\\theta_C\end{Bmatrix}=\begin{Bmatrix}-20.000 \\ 13.333 \\ -40.000 \\ 0 \\ 0\end{Bmatrix}$$

Solving:

$$\begin{aligned} v_B&=-\frac{230.00}{EI}\ \text{m},&\theta_B&=-\frac{34.167}{EI}\ \text{rad}\\ v_D&=-\frac{207.083}{EI}\ \text{m},&\theta_D&=\frac{63.958}{EI}\ \text{rad},\quad \theta_C=\frac{123.333}{EI}\ \text{rad} \end{aligned}$$

Deflection at B = 230.0/EI m downward (with $EI$ in kN m$^2$; e.g. $EI=2\times10^4$ kN m$^2$ gives 11.5 mm).

Element end forces $\{f\}=[k]\{d\}-\{F_{eq}\}$

Reactions: $R_A=50.312$ kN (up), $M_A=82.50$ kN m (counter-clockwise), $R_C=29.687$ kN. Check: $50.3125+29.6875=80$ kN = $10\times4+40$.

Bending moment (sagging +) and shear force

• A: $-82.50$ (hogging); at 2 m from A: $-1.875$; B: $+38.75$; D: $+59.375$; C: $0$.
• Shear: $+50.31$ at A, $+10.31$ at B (just before and after), $-29.69$ right of D.

 A(fixed)      B          D         C(roller)
 |-------------|----------|---------|
 BMD: -82.5   +38.75     +59.375     0   (kN m)
 SFD: +50.31 -> +10.31 (A to D) -> -29.69 (D to C)

Stresses at B

At B, $M_B=38.75$ kN m (sagging) and $V_B=10.31$ kN. The bending stress at a fibre at distance $y$ from the neutral axis is

$$\sigma_B=\frac{M_B\,y}{I}=\frac{38.75\times10^6\ y}{I}\ \text{N/mm}^2$$

tension at the bottom and compression at the top. The section size is not given, so for a rectangular section $b\times h$ use $\sigma_{max}=6M_B/(bh^2)$; the average shear stress is $\tau=V_B/A$ (maximum $1.5V_B/A$ for a rectangle).

Answer: $v_B=-230.00/EI$ m (downward), $\theta_B=-34.167/EI$, $M_B=38.75$ kN m; $\sigma_B=M_By/I$.

A constant strain triangle (CST) is used with plane stress ($E=69900$ MPa, $\nu=0.3$, $t=8$ mm). Units: N, mm. Node 1 is the loaded node; the nodes 2 and 3 on the left edge are fixed.

Step 1: Geometry

• Node 1: $(35,\ 25)$ mm
• Node 2: $(0,\ 25)$ mm
• Node 3: $(0,\ 0)$ mm

Area: $A=\tfrac12[x_1(y_2-y_3)+x_2(y_3-y_1)+x_3(y_1-y_2)]=437.5\ \text{mm}^2$

$b_1=y_2-y_3=25,\ b_2=y_3-y_1=-25,\ b_3=y_1-y_2=0$ mm

$c_1=x_3-x_2=0,\ c_2=x_1-x_3=35,\ c_3=x_2-x_1=-35$ mm

Step 2: [B] and [D] matrices (plane stress)

$$[B]=\frac{1}{2A}\begin{bmatrix}25 & 0 & -25 & 0 & 0 & 0 \\ 0 & 0 & 0 & 35 & 0 & -35 \\ 0 & 25 & 35 & -25 & -35 & 0\end{bmatrix},\quad 2A=875\ \text{mm}^2$$ $$[D]=\frac{E}{1-\nu^2}\begin{bmatrix}1&\nu&0\\ \nu&1&0\\0&0&\frac{1-\nu}{2}\end{bmatrix}=10^{4}\begin{bmatrix}7.681 & 2.304 & 0 \\ 2.304 & 7.681 & 0 \\ 0 & 0 & 2.688\end{bmatrix}\ \text{N/mm}^2$$

Step 3: Element stiffness matrix

$$[k]=t\,A\,[B]^T[D][B]= 10^{5}\begin{bmatrix}2.195 & 0 & -2.195 & 0.922 & 0 & -0.922 \\ 0 & 0.768 & 1.075 & -0.768 & -1.075 & 0 \\ -2.195 & 1.075 & 3.700 & -1.997 & -1.506 & 0.922 \\ 0.922 & -0.768 & -1.997 & 5.070 & 1.075 & -4.302 \\ 0 & -1.075 & -1.506 & 1.075 & 1.506 & 0 \\ -0.922 & 0 & 0.922 & -4.302 & 0 & 4.302\end{bmatrix}\ \text{N/mm}$$

(DOF order: $u_1,v_1,u_2,v_2,u_3,v_3$.)

Step 4: Load vector

Loads act at node 1 only. Assumed directions: 45 N horizontal toward $+x$ and 90 N vertical downward ($-y$). Self-weight is neglected (not given).

$$\{F\}=\begin{Bmatrix}45 \\ -90 \\ 0 \\ 0 \\ 0 \\ 0\end{Bmatrix}\ \text{N}$$

Boundary conditions and solution

The left edge is fixed, so nodes 2 and 3 have $u=v=0$. Only $u_1$ and $v_1$ are unknown, so use rows and columns 1 and 2 of $[k]$.

$$10^{5}\begin{bmatrix}2.195 & 0 \\ 0 & 0.768\end{bmatrix}\begin{Bmatrix}u_1 \\ v_1\end{Bmatrix}=\begin{Bmatrix}45 \\ -90\end{Bmatrix}$$ $$\begin{aligned} u_1&=2.05\times10^{-4}\ \text{mm} \\ v_1&=-0.001172\ \text{mm} \end{aligned}$$

Strains and stresses (constant over the element, so equal at the centroid)

Centroid: $(11.67,\ 16.67)$ mm.

$$\{\varepsilon\}=[B]\{d\}=\begin{Bmatrix}5.858\times10^{-6} \\ 0 \\ -3.348\times10^{-5}\end{Bmatrix}$$ $$\{\sigma\}=[D]\{\varepsilon\}=\begin{Bmatrix}0.45 \\ 0.135 \\ -0.9\end{Bmatrix}\ \text{MPa}$$

Answer: deflection at the load point (node 1): $u_1=2.05\times10^{-4}$ mm (horizontal), $v_1=-0.00117$ mm (vertical, downward). Stresses at the centroid: $\sigma_x=0.450$, $\sigma_y=0.135$, $\tau_{xy}=-0.900$ N/mm$^2$.

Element and coordinates

A four-node quadrilateral has nodes 1 to 4 anticlockwise and 8 DOFs $(u_i,v_i)$. Shape functions are derived for the parent square in natural coordinates $(\xi,\eta)$, $-1\le\xi,\eta\le1$, with nodes

$$1:(-1,-1),\quad 2:(1,-1),\quad 3:(1,1),\quad 4:(-1,1)$$

   4-----------3          eta
   |           |           ^
   |     +     |           |
   |           |           +--> xi
   1-----------2

Assumed displacement

A bilinear polynomial (4 terms, one for each node) is used for each displacement component:

$$u=\alpha_1+\alpha_2\xi+\alpha_3\eta+\alpha_4\xi\eta$$

Applying the nodal values:

$$\begin{Bmatrix}u_1\\u_2\\u_3\\u_4\end{Bmatrix}=\begin{bmatrix}1&-1&-1&1\\1&1&-1&-1\\1&1&1&1\\1&-1&1&-1\end{bmatrix}\begin{Bmatrix}\alpha_1\\\alpha_2\\\alpha_3\\\alpha_4\end{Bmatrix}$$

Inverting:

$$\begin{Bmatrix}\alpha_1\\\alpha_2\\\alpha_3\\\alpha_4\end{Bmatrix}=\frac14\begin{bmatrix}1&1&1&1\\-1&1&1&-1\\-1&-1&1&1\\1&-1&1&-1\end{bmatrix}\begin{Bmatrix}u_1\\u_2\\u_3\\u_4\end{Bmatrix}$$

Substituting into $u=[1\ \xi\ \eta\ \xi\eta]\{\alpha\}$ and collecting coefficients of the nodal values:

$$\begin{aligned} N_1&=\tfrac14(1-\xi-\eta+\xi\eta)=\tfrac14(1-\xi)(1-\eta)\\ N_2&=\tfrac14(1+\xi-\eta-\xi\eta)=\tfrac14(1+\xi)(1-\eta)\\ N_3&=\tfrac14(1+\xi+\eta+\xi\eta)=\tfrac14(1+\xi)(1+\eta)\\ N_4&=\tfrac14(1-\xi+\eta-\xi\eta)=\tfrac14(1-\xi)(1+\eta) \end{aligned}$$

General form:

$$\boxed{N_i=\tfrac14(1+\xi\xi_i)(1+\eta\eta_i)},\quad i=1,\dots,4$$

Properties

• $N_i=1$ at node $i$ and $0$ at the other three nodes.
• $\sum N_i=1$ (rigid-body motion), and $\sum N_ix_i=x$, so a linear field is reproduced exactly.
• Along each side the variation is linear, so neighbouring elements are compatible.

Use for a general quadrilateral

For an arbitrary quadrilateral the same $N_i$ are used as isoparametric functions: $x=\sum N_ix_i$, $y=\sum N_iy_i$, $u=\sum N_iu_i$. The derivatives with respect to $x,y$ are obtained with the Jacobian:

$$[J]=\begin{bmatrix}\partial x/\partial\xi&\partial y/\partial\xi\\\partial x/\partial\eta&\partial y/\partial\eta\end{bmatrix},\qquad \begin{Bmatrix}\partial N/\partial x\\\partial N/\partial y\end{Bmatrix}=[J]^{-1}\begin{Bmatrix}\partial N/\partial\xi\\\partial N/\partial\eta\end{Bmatrix}$$

For a rectangle of sides $2a\times2b$ centred at the origin, $\xi=x/a$ and $\eta=y/b$.

Data: $E=210$ GPa $=210\times10^6$ kN/m$^2$, $I=6\times10^6$ mm$^4=6\times10^{-6}$ m$^4$, so $EI=1260$ kN m$^2$. Nodes: 1 = fixed end (x = 0), 2 = point load (x = 3 m), 3 = first roller (x = 5 m), 4 = last roller (x = 10 m). Elements: 1-2 ($L=3$ m, $EI$), 2-3 ($L=2$ m, $EI$), 3-4 ($L=5$ m, $2EI$, UDL 10 kN/m). Units: kN, m.

Element stiffness

$$[k]=\frac{EI}{L^3}\begin{bmatrix}12&6L&-12&6L\\6L&4L^2&-6L&2L^2\\-12&-6L&12&-6L\\6L&2L^2&-6L&4L^2\end{bmatrix},\quad \{d\}=\{v_i,\theta_i,v_j,\theta_j\}^T$$

(Upward forces and displacements, and counter-clockwise moments and rotations, are positive.) Equivalent loads for the UDL element: $-25$ kN and $-20.833$ kN m at node 3; $-25$ kN and $+20.833$ kN m at node 4. The 20 kN load is a nodal load at node 2.

Assembled equation

Boundary conditions: $v_1=\theta_1=0$, $v_3=0$, $v_4=0$. Free DOFs: $v_2,\theta_2,\theta_3,\theta_4$.

$$\begin{bmatrix}2450 & 1050 & 1890 & 0 \\ 1050 & 4200 & 1260 & 0 \\ 1890 & 1260 & 4536 & 1008 \\ 0 & 0 & 1008 & 2016\end{bmatrix}\begin{Bmatrix}v_2\\\theta_2\\\theta_3\\\theta_4\end{Bmatrix}=\begin{Bmatrix}-20.000 \\ 0 \\ -20.833 \\ 20.833\end{Bmatrix}$$

(kN, m units.) Solving:

$$v_2=-4.329\ \text{mm},\quad \theta_2=0.00309\ \text{rad},\quad \theta_3=-0.00669\ \text{rad},\quad \theta_4=0.01368\ \text{rad}$$

Deflection under the 20 kN load = 4.33 mm downward; slope = 0.00309 rad (3.09 x 10$^{-3}$, counter-clockwise).

End forces and bending moments

Reactions: fixed end 5.018 kN and moment 6.230 kN m; first roller 44.210 kN; last roller 20.772 kN. Sum $=70$ kN = $20+10\times5$ (check).

Bending moment diagram (sagging +)

 x (m):   0        3        5                 7.92      10
          |--------|--------|-----------------|---------|
 BMD:    -6.23    +8.82    -21.14            +21.57      0   (kN m)