Archive | September 2013

Something New

A sneak peak of things to come…

The Method of Joints: Part 2

See Part 1.

With the theory established, the python script was simple enough to write. Prior to this project I hadn’t used numpy before, but after seeing some examples online involving matrix operations, I figured numpy was the perfect tool.

The script is split into 9 major steps:

  1. Read data from the input files
  2. Construct the members table, filled with the coefficients l_{ij} and m_{ij}
  3. Construct the coefficient matrix C and the vector of applied forces P
  4. Solve the matrix equation!
  5. Find the maximum safe load for the truss
  6. Present the results

All the relevant data are kept in a few separate files, which are imported into some ndarrays:

  • containing the (x,y) coordinates of each numbered joint, in mm
  • a list of members with their two associated joints
  • the two supported joints and whether it is a pin or roller support
  • the forces applied to the truss, and their locations
  • the physical properties of the material we’re using

The first four use numpy.genfromtxt() to form their arrays, but is imported into a dictionary:

material = {}
with open(inpath+'', 'r') as fin:
    for line in fin:
        if (line[0] != '#'):
            spl = line.split()
            if len(spl)==2:
                (key, value) = spl 
                material[key] = float(value)

This way the material properties can be accessed easily, for example the Young’s modulus is material['E']

The members table initially has three columns: the member number, the number of joint i, and the number of joint j. The table is extended with a few empty columns, which are filled in with the relevant geometric properties of the member (taken from the (x,y) coordinates of it’s two joints): dx, dy, L, l_{ij}, m_{ij}.

Constructing C is pretty much as described in Part 1: the l_{ij} and m_{ij} values for each member are inserted in the appropriate place, and then the coefficients for the reaction forces are added in to the final three columns. These coefficients are \pm1 as they’re aligned with the axes. The applied forces are inserted into a vector “P
Solving the matrix equation is the easy part and it’s why I chose numpy in the first place:

Q = numpy.linalg.solve(C,P)

I’ve been setting up the applied forces to sum to 1, so Q now contains the portion of the load which is carried by each member. For example if Q[3] = -0.707, this means member #3 is carrying a compressive force of 0.707*(Load).

The maximum load for the truss as a whole is found by iterating through the members and finding the maximum load for the bridge before that member fails, then taking the lowest of these failure limits. This isn’t a rigorous analysis and only yielding and Euler buckling are accounted for, but that’s good enough for a bridge made of plastic straws:

maxL_t = np.zeros([members.shape[0],3])
maxL = 1e6
for k in range(members.shape[0]):
    if Q[k] != 0:
        # yield: F/A &lt;= (safety factor) * s_y
        maxL_t[k,0] = material['sf'] * material['s_y'] * material['A'] / abs(Q[k])
        # buckling: F &lt;= (safety factor) * pi^2 * E * I / l^2
        if Q[k] &lt; 0:
            maxL_t[k,1] = material['sf'] * numpy.pi**2 * material['E'] * material['I'] \
                            / (abs(Q[k]) * members[k,5]**2)
            # which is the smaller limit?
            maxL_t[k,2] = min(maxL_t[k,0], maxL_t[k,1])
            maxL_t[k,2] = maxL_t[k,0]
        if maxL_t[k,2] &lt; maxL:
            maxL = maxL_t[k,2]

And with that the analysis is complete. I’ve set up two modes of output: a plaintext mode, where the results are sent straight to a text file; and a presentation mode, where matplotlib is used to draw the bridge and the results are written to an .html page like this one!

The Method of Joints: Part 1

In one of my courses we’ve been tasked with building a truss bridge out of plastic tubes, slightly thicker than drinking straws, which must span 1 metre and hold 5kg of load. My team had a few basic designs that we wanted to start off with, which meant a whole bunch of calculation to see which one was best. Of course I’m not one to pass up an opportunity like that, and solving systems of linear equations is something that computers are kind-of alright at, so I wrote a python script that solves a given truss and returns the internal forces in each member. With the script complete, a new bridge design takes only about 2 minutes from data input to final results! The slowest step is sketching the bridge design…



The method of joints involves stepping through each joint in a truss and solving the equations of static equilibrium, i.e.

    \begin{align*} \sum{F_x} &= 0,\\ \sum{F_y} &= 0. \end{align*}

These two equations must be solved for each joint, giving 2n equations, where n is the number of joints present. The procedure is quite straightforward but having more than a few joints makes it get tedious very quickly. So, let’s make it faster.

A simple truss

Consider the frame above (with members 1 & 2 the same length). We have six unknowns to solve for: three reaction forces from the supports, and the three internal forces. We also have six equations describing the system, how convenient. These are:

    \begin{align*} 1:\sum{F_x} &= R_{1x} + F_3 \cdot \frac{1}{\sqrt{2}} = 0 \\ 1:\sum{F_y} &= R_{1y} - F_1 - F_3 \cdot \frac{1}{\sqrt{2}} = 0 \\ 2:\sum{F_x} &= R_{2x} + F_2 = 0 \\ 2:\sum{F_y} &= F_1 = 0 \\ 3:\sum{F_x} &= -F_2 - F_3 \cdot \frac{1}{\sqrt{2}} = 0 \\ 3:\sum{F_y} &= F_3 \cdot \frac{1}{\sqrt{2}} - 100N = 0 \\ &\Rightarrow F_3 \cdot \frac{1}{\sqrt{2}} = 100N \end{align*}

Now this system is simple enough to solve by inspection. Ignoring this though, we can represent this system as a matrix equation C \cdot F = P, where P is a vector of the applied (or known) forces, F is a vector of our unknowns, and C is the coefficient matrix. Thus:

    \begin{equation*} \begin{bmatrix} 0 & 0 & \frac{1}{\sqrt{2}} & 1 & 0 & 0 \\[0.3em] -1 & 0 & -\frac{1}{\sqrt{2}} & 0 & 1 & 0 \\[0.3em] 0 & 1 & 0 & 0 & 0 & 1 \\[0.3em] 1 & 0 & 0 & 0 & 0 & 0 \\[0.3em] 0 & -1 & -\frac{1}{\sqrt{2}} & 0 & 0 & 0 \\[0.3em] 0 & 0 & \frac{1}{\sqrt{2}} & 0 & 0 & 0 \end{bmatrix} \cdot \begin{bmatrix} F_1 \\[0.3em] F_2 \\[0.3em] F_3 \\[0.3em] R_{1x} \\[0.3em] R_{1y} \\[0.3em] R_{2x} \end{bmatrix} = \begin{bmatrix} 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 100 \end{bmatrix} \end{equation*}

As long as the matrix C is non-singular, we can find it’s inverse and recover the solution F = C^{-1} \cdot P. In this case:

    \begin{equation*} \begin{bmatrix} F_1 \\[0.3em] F_2 \\[0.3em] F_3 \\[0.3em] R_{1x} \\[0.3em] R_{1y} \\[0.3em] R_{2x} \end{bmatrix} &= \begin{bmatrix} 0 & 0 & \frac{1}{\sqrt{2}} & 1 & 0 & 0 \\[0.3em] -1 & 0 & -\frac{1}{\sqrt{2}} & 0 & 1 & 0 \\[0.3em] 0 & 1 & 0 & 0 & 0 & 1 \\[0.3em] 1 & 0 & 0 & 0 & 0 & 0 \\[0.3em] 0 & -1 & -\frac{1}{\sqrt{2}} & 0 & 0 & 0 \\[0.3em] 0 & 0 & \frac{1}{\sqrt{2}} & 0 & 0 & 0 \end{bmatrix}^{-1} \cdot \begin{bmatrix} 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 100 \end{bmatrix} = \begin{bmatrix} 0 \\[0.3em] -100 \\[0.3em] 141.4 \\[0.3em] -100 \\[0.3em] 100 \\[0.3em] 100 \end{bmatrix} \end{equation*}

That was easy! Now comes the question of somehow constructing the matrix for an arbitrary truss, without actually writing the equations by hand.

First, note that for a simple truss (constructed only of triangles) with n joints, there are 2n-3 members connecting them. If the truss is supported by three support forces, for example R_{1x}, R_{1y}, R_{2x} above, then the system is statically determinate: there are 2n equations describing the system and 2n unknowns to solve for.

An arbitrary member

Now consider an arbitrary member k in a truss, connected between two joints i and j. The equations

    \begin{equation*} i:\sum{F_x}, \quad i:\sum{F_y}, \quad j:\sum{F_x}, \quad j:\sum{F_y} \end{equation*}

are the only equations with any contribution from member k (as a member can only connect to two joints), and we can number these equations as 2i-1, 2i, 2j-1, and 2j respectively. The member k will appear in exactly these four equations, although the coefficient may be zero as with F_1 in the example above.

Looking at the x-component of forces at joint i, (equation number 2i-1):

    \begin{equation*} \dots + cos \theta \cdot F_k + \dots + P_{ix} = 0. \end{equation*}

We can simplify this a little by using ratios of the member’s dimensions rather than it’s angle. By defining

    \begin{align*} l_{ij} = cos \theta = \frac{\Delta x}{L_k} \\ m_{ij} = sin \theta = \frac{\Delta y}{L_k} \end{align*}

and by moving the known/applied force P_ix to the right hand side of the equation, we get the following:

    \begin{alignat*}{4} & \dots + l_{ij} \cdot F_k + \dots && = && -P_{ix} \quad && (eqn. \; 2i-1) \\ & \dots + m_{ij} \cdot F_k + \dots && = && -P_{iy} \quad && (eqn. \; 2i) \\ & \dots - l_{ij} \cdot F_k + \dots && = && -P_{jx} \quad && (eqn. \; 2j-1) \\ & \dots - m_{ij} \cdot F_k + \dots && = && -P_{jy} \quad && (eqn. \; 2j) \end{alignat*}

In the matrix C, the entire column k will be zero except for the rows corresponding to these equations. So the contributions to C from member k are:

    \begin{alignat*}{3} & C_{2i-1,k} \quad && = \quad && l_{ij} \\ & C_{2i,k} \quad && = \quad && m_{ij} \\ & C_{2j-1,k} \quad && = \quad && -l_{ij} \\ & C_{2j,k} \quad && = \quad && -m_{ij} \end{alignat*}

With C complete, we just construct the vector F as above, with F_k in order, and the three support reaction forces at the end. The vector of applied forces P is simple enough: for a load L applied at node i we just set

    \begin{alignat*}{3} & P_{2i-1} && = && -L_x \\ & P_{2i} && = && -L_y \end{alignat*}

And now we just find C^{-1}, solve F = C^{-1} \cdot P, and we’re done!