Myrkgard logo
Series logo MECH000 Statically Determinate Trusses
Fri, 24 Nov 2017 18:00:00 +0000

Statically Determinate Trusses

In this episode we'll see how to find the unknown forces for a given statically determinate truss under external load. In contrast to the various manual methods, the method presented is general enough to be implemented in software easily. It consists of two steps: deriving a matrix equation from the problem description and then solving that equation. The episode comes with a code example. The example relies on the linear algebra C++ library Eigen for solving the matrix equation using a QR decomposition.

How it Works

We'll explain the general method with an example. In the following drawing, the left side shows our given problem. It's a truss consisting of 3 two force members, connected by 3 joints, and restricted by a pinned and a roller support. The two force members we'll call elements, and the joints we'll call nodes. The structure is loaded with a single applied force. We want to know the element and reaction forces. plot

On the right side the descretized truss is depicted, along with the names we'll use for the various forces. Here, we'll use roman numerals \(I, II, \cdots\) as indices for the nodes and uppercase latin characters \(A, B, \cdots\) for the elements. Up to 2 reaction forces can be present on each node. They are distinguished by the arabic numerals \(1, 2\). Each force vector has an associated angle, measured against the \(x\) axis, around the \(z\) axis. For the element forces the letter \(s\) is used, the reaction forces are denoted by \(r\), the applied forces by \(f\). For their respective angles greek letters \(\alpha, \beta, \cdots\) are used.

Given the number of nodes \(n_n\), the number of elements \(n_e\), and the number of reaction forces \(n_r\), the system is statically determinate if and only if the condition \[n_e + n_r = 2n_n \] is met. If that condition is not met, the problem is not solvable with this method.

The equilibrium conditions for the 3 nodes are described as a system of vector equations: \[\begin{array}{|lclclclclclclclc} \mathbf{s_{I,A}} & + & \mathbf{s_{I,B}} & + & \mathbf{r_{I,1}} & & & = & \mathbf{0} \\ \mathbf{s_{II,A}} & + & \mathbf{s_{II,C}} & + & \mathbf{f_{II}} & & & = & \mathbf{0} \\ \mathbf{s_{III,B}} & + & \mathbf{s_{III,C}} & + & \mathbf{r_{III,1}} & + & \mathbf{r_{III,2}} & = & \mathbf{0} \\ \end{array} \] We move the applied forces to the right-hand side and split the vector equations into their scalar \(x\) and \(y\) components: \[\begin{array}{|lclclclclclclclc} s_{I,A,x} & + & s_{I,B,x} & + & r_{I,1,x} & & & = & 0 \\ s_{I,A,y} & + & s_{I,B,y} & + & r_{I,1,y} & & & = & 0 \\ s_{II,A,x} & + & s_{II,C,x} & + & f_{II,x} & & & = & 0 \\ s_{II,A,y} & + & s_{II,C,y} & + & f_{II,y} & & & = & 0 \\ s_{III,B,x} & + & s_{III,C,x} & + & r_{III,1,x} & + & r_{III,2,x} & = & 0 \\ s_{III,B,y} & + & s_{III,C,y} & + & r_{III,1,y} & + & r_{III,2,y} & = & 0 \\ \end{array} \] Using the directional angles of the forces we get \[\begin{array}{|lclclclclclclclc} s_A cos(\alpha_{I,A}) + s_B cos(\alpha_{I,B}) + r_{I,1} cos(\beta_{I,1}) & = & 0 \\ s_A sin(\alpha_{I,A}) + s_B sin(\alpha_{I,B}) + r_{I,1} sin(\beta_{I,1}) & = & 0 \\ s_A cos(\alpha_{II,A}) + s_C cos(\alpha_{II,C}) & = & -f_{II} cos(\gamma_{II}) \\ s_A sin(\alpha_{II,A}) + s_C sin(\alpha_{II,C}) & = & -f_{II} sin(\gamma_{II}) \\ s_B cos(\alpha_{III,B}) + s_C cos(\alpha_{III,C}) + r_{III,1} cos(\beta_{III,1}) + r_{III,2} cos(\beta_{III,2}) & = & 0 \\ s_B sin(\alpha_{III,B}) + s_C sin(\alpha_{III,C}) + r_{III,1} sin(\beta_{III,1}) + r_{III,2} sin(\beta_{III,2}) & = & 0 \\ \end{array} \] And finally, we rewrite the system of equations as the single matrix equation \[ \begin{pmatrix} cos(\alpha_{I,A}) & cos(\alpha_{I,B}) & 0 & cos(\beta_{I,1}) & 0 & 0 \\ sin(\alpha_{I,A}) & sin(\alpha_{I,B}) & 0 & sin(\beta_{I,1}) & 0 & 0 \\ cos(\alpha_{II,A}) & 0 & cos(\alpha_{II,C}) & 0 & 0 & 0 \\ sin(\alpha_{II,A}) & 0 & sin(\alpha_{II,C}) & 0 & 0 & 0 \\ 0 & cos(\alpha_{III,B}) & cos(\alpha_{III,C} & 0 & cos(\beta_{III,1}) & cos(\beta_{III,2}) \\ 0 & sin(\alpha_{III,B}) & sin(\alpha_{III,C} & 0 & sin(\beta_{III,1}) & sin(\beta_{III,2}) \\ \end{pmatrix} \begin{pmatrix} s_A \\ s_B \\ s_C \\ r_{I,1} \\ r_{III,1} \\ r_{III,2} \\ \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ -f_{II} cos(\gamma_{II}) \\ -f_{II} sin(\gamma_{II}) \\ 0 \\ 0 \\ \end{pmatrix}\] And its general form \[\mathbf{K} \mathbf{a} = \mathbf{f}\]

For a 2-dimensional truss, \(\mathbf{K}\) is a \(2n_n \times 2n_n\) matrix. Vector \(\mathbf{a}\) contains the (unknown) element forces followed by the (unknown) reaction forces. The (known) applied forces are contained in \(\mathbf{f}\).

We can know use one of the standard methods to solve the equation, e.g. using the QR decomposition. If \(\mathbf{K}\) is not invertible, that means that the system is kinematically indeterminate, and thus not solvable. In numerical terms, the determinant of K must be reasonably far away from zero, e.g. \[\left\lVert det(\mathbf{K}) \right\lVert \overset{!}{>} 0.0001\] Here, \(\mathbf{a}\) evaluates to \[\begin{pmatrix} 707.107 \\ 0.000 \\ 707.107 \\ -707.107 \\ -707.107 \\ 0.000 \\ \end{pmatrix} \]

The result is also shown in the plot below. One can see that the elements \(A\) and \(B\) are loaded with tension, element \(C\) isn't loaded at all (so no elements in this example are loaded with pressure), and the pinned support is only loaded in one direction. plot

Code Example

The following code does the same what we've done manually above, also the truss configuration and the forces are all the same. It uses Eigen's built-in function to compute the QR decomposition of \(\mathbf{K}\) first and then solve the matrix equation with it.

#include <iostream>
#include <vector>
#define _USE_MATH_DEFINES
#include <cmath>
#include <Eigen/Dense>
#include <Eigen/QR>

struct Node
{
    double position_x = 0.; // meters
    double position_y = 0.; // meters
};

struct Element
{
    int node_a_idx = 0;
    int node_b_idx = 0;
};

struct ReactionForce
{
    int node_idx = 0;
    double angle = 0.0; // radians
};

struct AppliedForce
{
    int node_idx = 0;
    double strength = 0.; // Newton
    double angle = 0.;    // radians
};

bool isStaticallyIndeterminate(int const ne, int const nr, int const n2)
{
    return  (n2 != ne + nr);
}

bool isKinematicallyIndeterminate(Eigen::MatrixXd const &K, double const threshold = 0.0001)
{
    return abs(K.determinant()) < threshold;
}

int main()
{
    std::vector<Node> nodes = {
        // position x (meters), position y (meters)
        {0., 1.}, // I
        {1., 1.}, // II
        {1., 0.}, // III
    };

    std::vector<Element> elements = {
        // node a, node b
        {0, 1}, // A
        {0, 2}, // B
        {1, 2}, // C
    };

    std::vector<ReactionForce> reactionForces = {
        // node, angle (radians)
        {0, 0.},      // I, 1
        {2, M_PI_2},  // III, 1
        {2, 0.},      // III, 2
    };

    std::vector<AppliedForce> appliedForces = {
        // node, strength (Newton), angle (radians)
        {1, 1000., M_PI_4},
    };

    int const nn = nodes.size();
    int const ne = elements.size();
    int const nr = reactionForces.size();
    int const n2 = nn * 2;

    if (isStaticallyIndeterminate(ne, nr, n2))
        throw 1;

    Eigen::MatrixXd K = Eigen::MatrixXd::Zero(n2, n2);
    for (int i = 0; i < nn; ++i) {
        auto const &nd = nodes.at(i);
        for (int j = 0; j < ne; ++j) {
            auto const &el = elements.at(j);
            if (el.node_a_idx == i || el.node_b_idx == i) {
                auto const nd_other_idx = (el.node_a_idx == i) ? el.node_b_idx : el.node_a_idx;
                auto const nd_other = nodes.at(nd_other_idx);
                auto const x2 = nd_other.position_x;
                auto const y2 = nd_other.position_y;
                auto const x1 = nd.position_x;
                auto const y1 = nd.position_y;
                auto const angle = atan2(y2 - y1, x2 - x1);
                K(i * 2, j)     = cos(angle);
                K(i * 2 + 1, j) = sin(angle);
            }
        }
        for (int j = 0; j < nr; ++j) {
            auto const &rf = reactionForces.at(j);
            if (rf.node_idx == i) {
                K(i * 2, ne + j)     = cos(rf.angle);
                K(i * 2 + 1, ne + j) = sin(rf.angle);
            }
        }
    }

    if (isKinematicallyIndeterminate(K))
        throw 2;

    Eigen::VectorXd f = Eigen::VectorXd::Zero(n2);
    for (auto const &af: appliedForces) {
        f(af.node_idx * 2)     += - cos(af.angle) * af.strength;
        f(af.node_idx * 2 + 1) += - sin(af.angle) * af.strength;
    }

    std::cout << K.colPivHouseholderQr().solve(f) << std::endl;
    /* output:
        707.107
    3.00515e-13
        707.107
       -707.107
       -707.107
    2.43665e-13
    */
}