Direct Stiffness Method: Discrete Finite Elements

A bar element carries only axial force. Take nodes 1 and 2 with axial displacements $u_1$ and $u_2$ and nodal forces $f_1$ and $f_2$.

  f1 -->  o===================o  --> f2
         node 1     L       node 2
         (u1)               (u2)

Direct method

The bar behaves like a spring of stiffness $k = \dfrac{AE}{L}$. The elongation is $\delta = u_2 - u_1$ and the axial force is $F = k\delta$.

Equilibrium of the element requires $f_1 = -F$ and $f_2 = +F$:

$$\begin{aligned} f_1 &= \frac{AE}{L}(u_1 - u_2) \\ f_2 &= \frac{AE}{L}(-u_1 + u_2) \end{aligned}$$

In matrix form:

$$\begin{Bmatrix} f_1 \\ f_2 \end{Bmatrix} = \frac{AE}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\begin{Bmatrix} u_1 \\ u_2 \end{Bmatrix}$$

Check by energy method

With linear shape functions $N_1 = 1 - x/L$, $N_2 = x/L$, the strain is $\varepsilon = \dfrac{du}{dx} = \dfrac{u_2 - u_1}{L}$, so $[B] = \dfrac{1}{L}[-1\ \ 1]$. Then

$$[k] = \int_0^L [B]^T AE [B]\,dx = AE L\cdot\frac{1}{L^2}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} = \frac{AE}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$$

Both methods give the same result.

Properties

• Symmetric, since $k_{12} = k_{21}$.
• Singular ($\det = 0$): the element can move as a rigid body until supports are applied.
• Each row sums to zero, which expresses equilibrium for rigid-body motion.

Element stiffness matrices

For a spring element, $[k^{(e)}] = k\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$.

 wall |--k1--o--k2--o--k3--o--> 600 N
      1      2      3      4
                    ^ 400 N

(a) Global stiffness matrix (N/mm), degrees of freedom $u_1, \dots, u_4$

$$[K] = \begin{bmatrix} 200 & -200 & 0 & 0 \\ -200 & 500 & -300 & 0 \\ 0 & -300 & 400 & -100 \\ 0 & 0 & -100 & 100 \end{bmatrix}$$

(b) Apply boundary condition $u_1 = 0$

Delete row and column 1:

$$\begin{bmatrix} 500 & -300 & 0 \\ -300 & 400 & -100 \\ 0 & -100 & 100 \end{bmatrix}\begin{Bmatrix} u_2 \\ u_3 \\ u_4 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 400 \\ 600 \end{Bmatrix}$$

From row 3: $-100u_3 + 100u_4 = 600$, so $u_4 = u_3 + 6$.

Adding rows 2 and 3 gives $-300u_2 + 300u_3 = 1000$ and so $u_3 = u_2 + 3.333$.

Row 1: $500u_2 - 300(u_2 + 3.333) = 0$ gives $200u_2 = 1000$, so $u_2 = 5.000$ mm.

Then $u_3 = 8.333$ mm and $u_4 = 14.333$ mm.

(c) Reaction

From row 1 of the full system: $R_1 = 200u_1 - 200u_2 = -200(5) = -1000$ N. The wall pulls back with 1000 N, which equals the total applied load $400 + 600$.

(d) Spring forces

Check at node 3: $1000 - 600 = 400$ N, equal to the applied load.

Answer: $u_2 = 5.00$ mm, $u_3 = 8.33$ mm, $u_4 = 14.33$ mm; reaction at node 1 = 1000 N (towards the wall); spring forces 1000 N, 1000 N and 600 N.

Data and element stiffness

$E = 200\times10^3$ N/mm$^2$; $L = 250$ mm for each element.

 |==[ A1=400 ]==o==[ A2=200 ]==o--> 20 kN
 1      L=250   2     L=250   3
               ^ 10 kN

$$\begin{aligned} k_1 &= \frac{A_1E}{L} = \frac{400(200\times10^3)}{250} = 3.2\times10^5\ \text{N/mm} \\ k_2 &= \frac{A_2E}{L} = \frac{200(200\times10^3)}{250} = 1.6\times10^5\ \text{N/mm} \end{aligned}$$

Global equations with $u_1 = 0$

$$10^5\begin{bmatrix} 3.2 + 1.6 & -1.6 \\ -1.6 & 1.6 \end{bmatrix}\begin{Bmatrix} u_2 \\ u_3 \end{Bmatrix} = \begin{Bmatrix} 10\,000 \\ 20\,000 \end{Bmatrix}$$

Adding the two equations: $3.2\times10^5\,u_2 = 30\,000$, so $u_2 = 0.09375$ mm.

From the second equation: $1.6\times10^5(u_3 - u_2) = 20\,000$, so $u_3 - u_2 = 0.125$ and $u_3 = 0.21875$ mm.

Stresses

$$\begin{aligned} \sigma_1 &= E\,\frac{u_2 - u_1}{L} = 200\times10^3\times\frac{0.09375}{250} = 75\ \text{MPa} \\ \sigma_2 &= E\,\frac{u_3 - u_2}{L} = 200\times10^3\times\frac{0.125}{250} = 100\ \text{MPa} \end{aligned}$$

Reaction

$R_1 = -k_1(u_2 - u_1) = -3.2\times10^5(0.09375) = -30\,000$ N.

Check: $\sigma_1A_1 = 75\times400 = 30$ kN and $\sigma_2A_2 = 100\times200 = 20$ kN, consistent with equilibrium ($10 + 20 = 30$ kN).

Answer: $u_2 = 0.0938$ mm, $u_3 = 0.2188$ mm; $\sigma_1 = 75$ MPa, $\sigma_2 = 100$ MPa (both tensile); reaction = 30 kN directed opposite to the loads.

Geometry

Both bars have length $L = \sqrt{2^2 + 3^2} = 3.606$ m $= 3606$ mm. Direction cosines (measured from the support towards joint 3):

        3 (2,3)
       /|\  \
      / | \  \  <- 10 kN
     /  |  \ \
    1---+---+-2     load 20 kN down at 3

Element stiffness contribution to joint 3

For a truss element, the part belonging to joint 3 is $\dfrac{AE}{L}\begin{bmatrix} c^2 & cs \\ cs & s^2 \end{bmatrix}$.

$\dfrac{AE}{L} = \dfrac{500(200\times10^3)}{3606} = 27\,732$ N/mm.

$$\begin{aligned} K_{xx} &= 27\,732(0.5547^2 + 0.5547^2) = 27\,732(0.6154) = 17\,068\ \text{N/mm} \\ K_{yy} &= 27\,732(0.8321^2 + 0.8321^2) = 27\,732(1.3846) = 38\,402\ \text{N/mm} \\ K_{xy} &= 27\,732[(0.5547)(0.8321) + (-0.5547)(0.8321)] = 0 \end{aligned}$$

Solve

$$\begin{bmatrix} 17\,068 & 0 \\ 0 & 38\,402 \end{bmatrix}\begin{Bmatrix} u_3 \\ v_3 \end{Bmatrix} = \begin{Bmatrix} 10\,000 \\ -20\,000 \end{Bmatrix}$$ $$u_3 = 0.5859\ \text{mm},\quad v_3 = -0.5208\ \text{mm}$$

Member forces

Axial elongation of a bar: $\delta = c\,u_3 + s\,v_3$ (since the other end is fixed). Force $= \dfrac{AE}{L}\delta$.

Check equilibrium at joint 3 (bars in compression push the joint away from their supports): horizontal components $3.005(0.5547) - 21.033(0.5547) = -10.0$ kN, balancing the applied $+10$ kN; vertical components $3.005(0.8321) + 21.033(0.8321) = 20.0$ kN, balancing the applied 20 kN downward.

Answer: $u_3 = 0.586$ mm, $v_3 = -0.521$ mm; bar 1-3 force $= 3.0$ kN (compression); bar 2-3 force $= 21.0$ kN (compression).

A beam element has 4 degrees of freedom: $\{d\} = \{v_1,\ \theta_1,\ v_2,\ \theta_2\}^T$ with corresponding nodal forces $\{f_1,\ m_1,\ f_2,\ m_2\}^T$.

   v1,f1               v2,f2
    |  (th1,m1)         |  (th2,m2)
    o===================o
    1         L         2

Displacement function

Since there are 4 nodal values, a cubic is chosen: $v(x) = a_1 + a_2x + a_3x^2 + a_4x^3$.

Applying $v(0) = v_1$, $v'(0) = \theta_1$, $v(L) = v_2$, $v'(L) = \theta_2$ and solving gives the shape functions (Hermite cubics):

$$\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}$$

Strain energy and stiffness

Curvature is $\kappa = \dfrac{d^2v}{dx^2} = [B]\{d\}$, where

$$[B] = \left[\frac{-6}{L^2} + \frac{12x}{L^3},\ \ \frac{-4}{L} + \frac{6x}{L^2},\ \ \frac{6}{L^2} - \frac{12x}{L^3},\ \ \frac{-2}{L} + \frac{6x}{L^2}\right]$$

The strain energy is $U = \dfrac{1}{2}\displaystyle\int_0^L EI\,\kappa^2\,dx$, which gives

$$[k] = EI\int_0^L [B]^T[B]\,dx$$

Evaluating the integrals term by term (for example $k_{11} = EI\int_0^L \left(\frac{12x}{L^3} - \frac{6}{L^2}\right)^2dx = \frac{12EI}{L^3}$):

$$[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}$$

Meaning of a column

Column 1 is the set of forces and moments needed to give a unit deflection $v_1 = 1$ at node 1 while $\theta_1 = v_2 = \theta_2 = 0$: shear $12EI/L^3$ and moment $6EI/L^2$ at node 1, and the opposite shear $-12EI/L^3$ and moment $6EI/L^2$ at node 2.

The matrix is symmetric and singular (rigid-body motion is possible until supports are applied).

Data

$EI = 200\times10^9 \times 8\times10^{-6} = 1.6\times10^6\ \text{N m}^2$, $L_e = 2$ m for each element.

$\dfrac{EI}{L_e^3} = \dfrac{1.6\times10^6}{8} = 2\times10^5$ N/m.

Nodes: 1 (left support), 2 (mid-span), 3 (right support). Degrees of freedom: $(v_1, \theta_1, v_2, \theta_2, v_3, \theta_3)$. Boundary conditions: $v_1 = v_3 = 0$.

   ^              10 kN
   |               |
   o======= (1) ===o======= (2) ===o
  1                2                3
  |<------ 2 m ---->|<------ 2 m -->|

Element matrix (each element)

With $L_e = 2$:

$$[k] = 2\times10^5\begin{bmatrix} 12 & 12 & -12 & 12 \\ 12 & 16 & -12 & 8 \\ -12 & -12 & 12 & -12 \\ 12 & 8 & -12 & 16 \end{bmatrix}$$

Reduced system

Free DOFs: $\theta_1, v_2, \theta_2, \theta_3$. By symmetry about mid-span, $\theta_2 = 0$ and $\theta_3 = -\theta_1$.

Row for $\theta_1$ (element 1, second row), with $v_1 = 0$, $\theta_2 = 0$:

$$2\times10^5\,(16\theta_1 - 12v_2) = 0$$

Row for $v_2$ (sum of element 1 third row and element 2 first row), with $v_1 = v_3 = 0$, $\theta_2 = 0$, $\theta_3 = -\theta_1$:

$$2\times10^5\,(-12\theta_1 + 24v_2 + 12\theta_3) = 2\times10^5\,(-24\theta_1 + 24v_2) = -10\,000$$

So $\theta_1 = 0.75v_2$, and $-18v_2 + 24v_2 = -0.05$ gives

$$v_2 = -0.008333\ \text{m},\qquad \theta_1 = -0.00625\ \text{rad}$$

Reactions

$R_1 = 2\times10^5(12v_1 + 12\theta_1 - 12v_2 + 12\theta_2) = 2.4\times10^6(-0.00625 + 0.008333) = 5000$ N, and by symmetry $R_3 = 5000$ N.

Comparison

Exact: $v_{max} = \dfrac{PL^3}{48EI} = \dfrac{10\,000(4)^3}{48(1.6\times10^6)} = 8.333$ mm; exact support rotation $\dfrac{PL^2}{16EI} = 6.25\times10^{-3}$ rad. The FE results match exactly because the cubic element is exact for a load applied at a node.

Answer: mid-span deflection = 8.33 mm downward; support rotations = 6.25 $\times10^{-3}$ rad (left clockwise, right anticlockwise); reactions = 5 kN each; the FE result equals the exact value.

Comparison (2D)

The local frame element matrix is the bar matrix for $u$ combined with the beam matrix for $v, \theta$:

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

Need for transformation

Members of a truss or frame lie at different angles, but the global equations must be written in one common (global) $x$-$y$ system. The element matrix derived along the member axis (local) is converted by

$$[k] = [T]^T[k'][T]$$

For a truss element, with $c = \cos\theta$ and $s = \sin\theta$:

$$[T] = \begin{bmatrix} c & s & 0 & 0 \\ 0 & 0 & c & s \end{bmatrix}$$

For a frame element, $[T]$ is a $6\times6$ matrix with the $\theta$ rotation unchanged. Beam elements all lie along one axis, so no transformation is needed.

Analogy

Element equation: $\dfrac{kA}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\begin{Bmatrix} T_i \\ T_j \end{Bmatrix} = \begin{Bmatrix} Q_i \\ Q_j \end{Bmatrix}$.

Element conductances ($A = 1\ \text{m}^2$)

$$c_1 = \frac{1.2}{0.25} = 4.8\ \text{W/K},\qquad c_2 = \frac{0.04}{0.10} = 0.4\ \text{W/K}$$

 800 C o-----brick-----o-----insulation-----o 30 C
       1       (1)     2         (2)         3

Assembly and solution

$$\begin{bmatrix} 4.8 & -4.8 & 0 \\ -4.8 & 5.2 & -0.4 \\ 0 & -0.4 & 0.4 \end{bmatrix}\begin{Bmatrix} T_1 \\ T_2 \\ T_3 \end{Bmatrix} = \begin{Bmatrix} Q_1 \\ 0 \\ Q_3 \end{Bmatrix}$$

With $T_1 = 800$ and $T_3 = 30$, the middle row (no heat source at node 2) gives

$$-4.8(800) + 5.2\,T_2 - 0.4(30) = 0 \Rightarrow T_2 = \frac{3840 + 12}{5.2} = 740.77\ ^\circ\text{C}$$

Heat flow through element 1: $q = 4.8(800 - 740.77) = 284.3$ W/m$^2$.

Check with element 2: $0.4(740.77 - 30) = 284.3$ W/m$^2$. Check with thermal resistances: $q = \dfrac{800 - 30}{0.25/1.2 + 0.10/0.04} = \dfrac{770}{2.7083} = 284.3$ W/m$^2$.

Answer: interface temperature $= 740.8\ ^\circ$C; heat flow $= 284.3$ W per m$^2$ of wall.