Introduction

The objective of this article is to present some introductory concepts in Finite Element Methods by solving a two-dimensional truss problem. The inspiration of this post is the first chapter of the book A First Course in Finite Elements by Jacob Fish and Ted Belytschko. The Julia programming language is used to implement the algorithm suggested by the book, and the results are complemented with those obtained by Ftool, a software for structural analysis.

Methods

This topic starts with the construction of the element and truss stiffness matrices and its implementation in Julia. Next, the treatment of boundary conditions and the solution for nodal displacements and forces are presented. Finally, The truss to be evaluated as an example is defined.

Bar element

The following image displays the configuration for a bar element and its two nodes, their displacements and the forces acting on them.

element

The local stiffness matrix $\mathbf{K}^e$ relates nodal forces and nodal displacements for the element $e$ in the global coordinate system $xy$ as

$$\eq{ \mathbf{F}^e = \mathbf{K}^e \mathbf{d}^e, }$$

where

$$\eq{ \mathbf{K}^e = \frac{E^e A^e}{l^e} \begin{bmatrix} c^2 & s c & -c^2 & -s c \\ s c & s^2 & -s c & -s^2 \\ -c^2 & -s c & c^2 & s c \\ -s c & -s^2 & s c & s^2 \\ \end{bmatrix}, }$$

and

$$ s^2 = \sin^2\phi^e, \quad c^2 = \cos^2\phi^e, \quad s c = \sin\phi^e \cos\phi^e. $$

The tensile stress is computed, after solving for $\mathbf{d}^e$, as

$$\eq{ \sigma^e = \frac{E^e}{l^e} \begin{bmatrix} -c & -s & c & s \end{bmatrix} \mathbf{d}^e. }$$

In the expressions above for the bar element, $E^e$ is the Young’s Modulus, $A^e$ is the cross-section area and $l^e$ is the element length.

The concepts above are implemented in the following Julia composite type:

Julia Element type
struct Element
    node1::Vector{Float64}
    node2::Vector{Float64}
    A::Float64
    E::Float64
    length::Float64
    R::Vector{Float64}
    K::Matrix{Float64}

    function Element(node1, node2, A, E)
        Δ = node2 - node1
        length = √sum(Δ.^2)
        sinϕ = Δ[2]/length
        cosϕ = Δ[1]/length
        s2 = sinϕ^2
        c2 = cosϕ^2
        sc = sinϕ*cosϕ

        R = Float64[-cosϕ, -sinϕ, cosϕ, sinϕ]
        
        K = A*E/length * Float64[
             c2  sc -c2 -sc;
             sc  s2 -sc -s2;
            -c2 -sc  c2  sc;
            -sc -s2  sc  s2;
        ]
        
        new(node1, node2, A, E, length, R, K)
    end
end

Truss

A truss is a collection of interconnected bar elements and its stiffness matrix is constructed by combining the contributions of each element. In order to build the truss stiffness matrix $\mathbf{K}$, each element stiffness matrix $\mathbf{K}^e$ is first transformed from its local representation displayed in equation $(2)$ to a global representation, and the elements contributions are added together as

$$\eq{ \mathbf{K} = \sum_{e=1}^{n_{el}} \mathbf{L}^{e \mathrm{T}} \mathbf{K}^e \mathbf{L}^e. }$$

$\mathbf{L}^e$ is called gather matrix, it has 4 rows (2 for each element node) and $2 n_\mathrm{nodes}$ columns, where $n_\mathrm{nodes}$ is the total number of nodes in the truss. This matrix is almost entirely composed of zeros, being ones only the entries that map the element nodes to their global numbering.

The following implementation in Julia shows the construction of the truss stiffness matrix in a composite type:

Julia Truss type
struct Truss
    nodes::Matrix{Float64}
    connectivity::Matrix{Int64}
    A::Float64
    E::Float64
    number_of_elements::Int64
    number_of_nodes::Int64
    elements::Vector{Element}
    K::Matrix{Float64}

    function Truss(nodes, connectivity, A, E)
        nel = size(connectivity, 1)
        nds = size(nodes, 1)
        elements = Array{Element}(undef, nel)
        K = zeros(Float64, 2*nds, 2*nds)
        for e in 1:nel
            i₁ = connectivity[e, 1]
            i₂ = connectivity[e, 2]

            element = Element(nodes[i₁,:], nodes[i₂,:], A, E)
            elements[e] = element

            L = zeros(Int8, 4, 2*nds)
            L[1:2, 2*i₁-1:2*i₁] = Matrix{Int8}(I, 2, 2)
            L[3:4, 2*i₂-1:2*i₂] = Matrix{Int8}(I, 2, 2)
            K += L' * element.K * L
        end
        new(nodes, connectivity, A, E, nel, nds, elements, K)
    end
end

Boundary value problem

Some of the truss nodes are fixed (essential nodes) and some are free to displace. In addition to that, external forces are applied to some of the nodes and these will be the loads that make the truss deform. To aid the solution process, the equation that relates forces and displacements is arranged in the following form:

$$\eq{ \begin{bmatrix} \mathbf{K}^\mathrm{E} & \mathbf{K}^\mathrm{EF} \\ \mathbf{K}^\mathrm{FE} & \mathbf{K}^\mathrm{F} \\ \end{bmatrix} \begin{bmatrix} \mathbf{d}^\mathrm{E} \\ \mathbf{d}^\mathrm{F} \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{f}^\mathrm{F} \end{bmatrix} + \begin{bmatrix} \mathbf{r}^\mathrm{E} \\ \mathbf{0} \end{bmatrix}, }$$

where the superscript $\mathrm{E}$ stands for essential and $\mathrm{F}$ stands for free. In this article, essential nodes are considered fixed, therefore $\mathbf{d}^\mathrm{E} = \mathbf{0}$. So, the displacement of free nodes $\mathbf{d}^\mathrm{F}$ can be computed as

$$\eq{ \mathbf{K}^\mathrm{F} \mathbf{d}^\mathrm{F} = \mathbf{f}^\mathrm{F}. }$$

After solving for $\mathbf{d}^\mathrm{F}$, the reactions at the essential nodes can be computed from

$$\eq{ \mathbf{r}^\mathrm{E} = \mathbf{K}^\mathrm{EF} \mathbf{d}^\mathrm{F}, }$$

and the tensile stress can be computed from $(3)$.

All these steps are implemented in Julia through the following function:

Julia solve function
function solve(truss::Truss, fixed_nodes::Vector{Int64}, forces::Dict{Int64, Vector{Float64}})
    nds = truss.number_of_nodes
    nel = truss.number_of_elements

    # Displacement vector.
    d = zeros(Float64, 2*nds)

    # Force + reaction vector.
    f = zeros(Float64, 2*nds)

    # Stress vector.
    σ = zeros(Float64, nel)

    # Masks for essential and free nodes.
    en = fill(false, 2*nds)  # Essential nodes
    for i in fixed_nodes
        en[2*i-1:2*i] .= true
    end
    fn = .!en  # Free nodes

    # Stiffness submatrices.
    KEE = truss.K[en, en]
    KEF = truss.K[en, fn]
    KFE = truss.K[fn, en]
    KFF = truss.K[fn, fn]

    # Force vector at free nodes.
    for i in keys(forces)
        f[2*i-1:2*i] = forces[i]
    end
    fF = f[fn]  # Forces at free nodes

    # Displacement of free nodes.
    dF = KFF\fF
    d[fn] = dF
    
    # Reaction of fixed nodes.
    rE = KEF * dF
    f[en] = rE

    # Stress vector.
    for (i, e) in enumerate(truss.elements)
        i₁ = truss.connectivity[i, 1]
        i₂ = truss.connectivity[i, 2]
        
        de = zeros(Float64, 4)
        de[1:2] = d[2*i₁-1:2*i₁]
        de[3:4] = d[2*i₂-1:2*i₂]

        σ[i] = e.E/e.length * sum(e.R .* de)
    end

    # Reshaped displacement and force vectors
    dr = reshape(d, 2, nds)'
    fr = reshape(f, 2, nds)'
    
    return dr, fr, σ
end

Application problem

The code implemented in Julia will be tested against the truss presented in following figure. I consists of 11 steel bar elements, all with circular cross-section with 2 cm of diameter. The complete truss has 3 m in length and 1 m in height. The elements are connected by 7 nodes, where nodes 1 and 4 are fixed, but free to rotate. Vertical loads of -0.5 kN and -1.0 kN are concentrated at nodes 2 and 3, respectively.

truss

In Julia, the truss is defined by the following code:

Julia truss definition
nodes = Float64[
    0.0 0.0;
    1.0 0.0;
    2.0 0.0;
    3.0 0.0;
    1.0 1.0;
    2.0 1.0;
    3.0 1.0;
]

connectivity = Int64[
    1 2;
    1 5;
    2 3;
    2 5;
    2 6;
    3 4;
    3 6;
    3 7;
    4 7;
    5 6;
    6 7;
]

# Circular cross section.
D = 0.02     # Diameter (m)
A = π*D^2/4  # Area (m²)

# Young's modulus (Pa).
E = 205e9  # Steel

# Truss definition.
truss = Truss(nodes, connectivity, A, E);

# Boundary conditions.
fixed_nodes = [1, 4]
forces = Dict(2 => [0.0, -500.0], 3 => [0.0, -1000.0]);

# Solution.
d, f, σ = solve(truss, fixed_nodes, forces)
T = σ * A;

In the code above, the nodes are defined by their x and y coordinates in the nodes matrix. Elements are defined by each line of the connectivity matrix, where the two columns represent the first and second nodes, respectively.

Results

The truss problem defined in the last topic is solved and the solution is complemented with results obtained by Ftool, a software for structural analysis. Their results will not be compared quantitatively side by side because they are actually identical. So, the following images, displaying results from the code implemented in Julia and from Ftool, represent the solution obtained throught the Finite Element Method.

The first two images displays the truss in deformed condition, scaled by a factor of 1000. The first one computed with Julia and plotted with Makie, the second evaluated by Ftool.

truss deformation
ftool truss deformation

The next two images shows the actual displacements in x and y, respectively.

ftool dx
ftool dx

Ftool does not compute the tensile stress directly, but rather the elements axial forces. These are displayed by the following two images, the first computed with Julia and plotted with Makie, the second evaluated by Ftool.

element load
ftool element load

Conclusion

This post presents an interesting and fun way of introducing Finite Element concepts. The truss problem is easily implemented in any programming language and showcases why the Finite Element Methods is so dominant in structural analysis. However, this introduction does not tell much what is FEM and what it can be used for. Future articles will hopefully make this clearer.

References

  1. Jacob Fish and Ted Belytschko. 2007. A First Course in Finite Elements. Wiley, Chichester, UK.

Appendices