Simply supported beam—Newmark’s method

This would have been the last post in the series if I hadn’t realized recently that another method deserved a post. So this is the penultimate. It was described in this 1943 paper by Nathan Newmark and is known—or was known—as Newmark’s method.

I say was because this method is dead. How dead? Forty-five years ago, when I was an undergraduate taking structural analysis courses, Newmark’s method wasn’t being taught in the department that he had been the head of for decades and whose main building is named after him. That’s dead.

Which isn’t to say it wasn’t a good method. Some of my professors talked about using it when they were young, and it was a practical technique when engineers used slide rules and desk calculators, but it was clear by the 80s that its time had passed. Still, it’s desncribed in the textbook I used as an undergrad, and I did eventually learn it well after I was out of school, so here it is.

Newmark’s method is basically a way of doing double integration, so you can use it to go from loading to bending moments or from bending moments to deflection. As Newmark says in the introduction to his paper,

The numerical procedure is approximate, but leads to exact moments (or deflections) when the loading diagram (or diagram of “angle changes”) is made up of segments that are bounded by straight lines or by arcs of parabolas.

For our simply supported, uniformly loaded beam, the “angle change” diagram—which we’ve been calling a curvature diagram—is a parabola, so we can get an exact solution.

Let’s say we have a beam with any kind of support conditions and some general distributed loading:

Beam with generic loading

The fundamental trick of Newmark’s method is to imagine a set of simply supported beams inserted between the load and the real beam.

Insertion of simply supported beams

The small imaginary beam segments are all of the same length, λ, and the analysis proceeds by calculating the reaction forces at the ends of these beams and then building a table of loads, shears, and moments based on these reactions (reversed in direction) acting on the large original beam.

Here’s a figure from the paper that explains how to calculate the reaction forces for parabolically distributed loading:

Equivalent concentrated loads for parabolas

There are simpler formulas if the loading is distributed linearly.

You may have noticed the double curvature in the right drawing, which means that the load shown there is not parabolic. Newmark suggests using the formulas in this figure even if the loading curve is of higher order—the parabolic formulas are close enough if you divide the beam into several segments.

I’ve focused here on the process of going from loads to bending moments because I think the mechanics of that is easier to understand. But the same process applies to going from curvature to deflection. We went through that in the conjugate beam post.

In fact, Newmark does our problem in the paper. He uses the constant q as the intensity of the load and splits the beam into four segments:

Newmark four-segment solution

This may look impenetrable, but we can simplify it by using only two segments and explaining each step. Here’s the table, starting with the moment diagram:

Newmark two-segment solution

The row labeled Moment is the value of the bending moment at each point. The Curvature row is just the moment at each point divided by EI and with the sign reversed. Those are the easy ones.

The Angle change row is based on the curvature values and the formulas in Figure 5 above. The value at the left end is

L/224[70+6(wL28EI)0]=wL364EI

where we’ve used λ=L/2. This is also the value at the right end.

The value at the center is

L/212[0+10(wL28EI)+0]=5wL396EI

We start the Average slope row with the knowledge that the slope at the center of the beam is zero, so if the angle change across the center is

5wL396EI

then the average slope must go from

5wL3192EIto5wL3192EI

as we move from the left half of the beam to the right half.

(If you’re having trouble with this, remember that angle change is analogous to shear and think of a beam with an antisymmetric shear diagram. If the change in shear across the centerline is F, then the shear must be F/2 in the left half of the beam and F/2 in the right half.)

For the Deflection row, we know that the deflection is zero at the left and right ends. Since the average slope in the left half of the beam is

5wL2192EI

the deflection at the center must be

(5wL3192EI)(L2)=5wL4384EI

an answer we’ve now seen twelve times.


Simply supported beam—dummy unit load method

We’re in the home stretch of this series, so ANIAT will soon go back to complaining about Apple’s UI choices.1 Next week’s WWDC keynote should provide some inspiration.

But today’s post covers our eleventh method for deriving the center deflection of a uniformly loaded simply supported beam: the dummy unit load method. When I was an undergraduate, this was the method we used to determine the deflections of truss structures, but it’s more general than that. My favorite explanation of why it works is in Nicholas J. Hoff’s The Analysis of Structures. I have the original Wiley hardcover of this book, but you can apparently get both paperback and hardcover reprints.

Hoff uses the principle of virtual work in his explanation, and I’d like to quote him here, but unfortunately his explanation is split between one section on the analysis of trusses and another on the analysis of beams. I couldn’t figure out a nice way to put them together coherently, so you’re just going to have to trust me.

For a beam, you can get the deflection at a particular point by following these steps:

  1. Work out an equation for the bending moment, M, in the beam for the given set of applied loads.
  2. Imagine that same beam with those loads removed (this is the dummy structure) and a concentrated unit load placed at the point at which we want the deflection. By definition, the magnitude of a unit load is one. The unit load points in the direction of the deflection we want to calculate.
  3. Work out an equation for the bending moment, m, in the dummy beam with the unit load.
  4. The desired deflection is

    δ=0LmMEIdx

Let’s apply this simple principle to our problem:

Simply supported beam with uniform load

The bending moment is

M=w2(Lxx2)

where the x coordinate starts at the left end of the beam and increases to the right, as usual.

This same structure with just a unit load at the center is

Beam with dummy unit load

Its bending moment is

m=12xfor 0xL2

and m is symmetric about the center of the beam.

Therefore,

δ=y(L2)=0LmMEIdx

and we can take advantage of symmetry to say

y(L2)=2w4EI0L/2(Lx2x3)dx=w2EI[13Lx314x4]0L/2=w2EI[L424L464]=5wL4384EI

You can pretty much do this problem in your head if you’re good at figuring out least common multiples. I’m not, but after writing things out, I did realize that the LCM of 24 and 64 is 192.


  1. Heads up, though. I recently thought of a thirteenth method that’s different enough from the others to merit another post. Sorry about that. 


Simply supported beam—finite element method

The biggest problem I face when writing posts in this series is deciding how much background explanation to put in. If I included only the stuff I do to get the answer, the posts would be very terse: one sketch and a dozen or fewer lines of math. A full explanation, though, could easily turn into several semesters’ worth of structural analysis theory. I’m trying to get close to the former while adding just enough background so readers with engineering experience can follow along. Even if you haven’t used the method of a particular post, you should be able to take what’s in the post and find the whatever further information you need to fill in the inevitable gaps.

This problem is most acute in describing today’s analysis. The finite element method—more precisely called matrix structural analysis or the direct stiffness method when applied to structures made of beam, truss, or frame elements—has many analytical parents: the Rayleigh-Ritz method, the slope-deflection equation, and matrix math. Presenting it without at least some of its antecedents would be a crime. But I will try to keep the introductory material to a minimum.

The fundamental piece of this analysis is the beam bending element. This is a length of beam with no interior loading, and its behavior can be described completely by four parameters, called its degrees of freedom. The DOF are the the displacements and rotations at the two ends (or nodes) of the element, usually numbered and directed as follows:

Beam element with DOF

The DOF are often called generalized displacements, and we’ll refer to them as u1, u2, u3, and u4. The deflected shape of this beam is completely determined by the sum of the products of these generalized displacements and a set of shape functions that look like this:

Finite element beam shape functions

Each shape function is a cubic curve that has a unit displacement corresponding to one of the DOFs and zero displacement corresponding to the other three DOFs. Apart from the sign convention, these are the same shapes we saw in our discussion of the slope-deflection equation.

There is a generalized force associated with each DOF; these are the forces and moments at the two ends of the beam, and we’ll call them f1, f2, f3, and f4. Note that when we multiply a u term by an f term, we get an energy expression. For DOF 1 and 3, that’s a displacement multiplied by a force; for DOF 2 and 4, that’s a rotation multiplied by a moment.

For this element, we can write a matrix equation for the relationship between the generalized displacements and generalized forces,

𝒌𝒖=𝒇

where

𝒖={u1u2u3u4},𝒇={f1f2f3f4}

and

𝒌=EIL3[ 12 6L 12 6L 6L 4L2 6L 2L2 12 6L 12 6L 6L 2L2 6L 4L2 ]

Each column is the set of generalized forces that correspond to that shape function.

You may notice that the terms in the second row of 𝒌 match the coefficients in the slope-deflection equation before the yA and yB terms were combined into the ψ span rotation term.

What should we do when, as in our case of a beam with a distributed load, there is interior loading? We convert that internal loading into what are called consistent loads at the ends. These are equal in magnitude and opposite in direction to the reaction forces and moments that would occur if the interior loads were applied to a beam with both ends fixed: wL/2 for DOF 1 and 3 and wL2/12 for DOF 2 and 4.

To solve a problem by the finite element method, we assemble the structural stiffness matrix, 𝑲, and the structural force vector, 𝑭 from all the individual elements and then solve the structural equation,

𝑲𝑼=𝑭

for the vector of structural displacements, 𝑼. As you can see, I tend to use lower-case letters for the element terms and upper-case for the structural terms. That’s a reasonably common, but by no means universal, notation.

In general, this assembly and solution of a matrix equation can be quite complex, which is why the finite element method is typically thought of as a process that requires a computer. But our problem is simple enough to do by hand.

Simply supported beam with uniform load

Recall that in both the slope-defection and Myosotis solutions, we took advantage of symmetry and considered just half the beam. We’ll do that again here.

Half beam for finite element analysis

This half-beam has just two DOF, the rotation at the left end and the deflection at the right end. In the assembly process, element DOF 2 matches structural DOF 1 and element DOF 3 matches structural DOF 2. The other element DOF don’t appear in the structural matrix equation because their U terms are zero. Here’s an illustration of the assembly process for the stiffness matrix, accounting for the fact that the beam element has a length of L/2:

Assembly of stiffness matrix

I took a class in matrix structural analysis (CE 361) from Prof. Jamshid Ghaboussi in the fall of 1981. He used the phrase destination array for the list of structural DOF that the element DOF go to. I’ve always thought that was a particularly good description and have used it ever since, even though I don’t think I’ve ever seen it in any finite element textbook. Note again that DOF 1 and 4 in the element have no destination because those displacements are zero in the structure.

When the assembly of both the stiffness matrix and force vector (which is done in a similar way) are complete, we get this equation:

8EIL3[ L2 3L 3L 12 ] { U1 U2 } = { wL248 wL4 }

There are different ways to solve this, but because a 2×2 matrix is easy to invert, that’s how I did it. Premultiplying each side by the inverse of 𝑲 gives us

{ U1 U2 } = L24EI [ 12 3L 3L L2 ] { wL248 wL4 }

Matrix multiplication gives us

{ U1 U2 } = { wL324EI 5wL4384EI }

The signs are the opposite of what we’ve seen before because of the directions we used in defining the degrees of freedom. Physically, though, these are the same answers we’ve seen several times.


Simply supported beam—Castigliano's method

Continuing our series on the many ways to get the center deflection of a uniformly loaded beam, we come to another energy-based technique: Castigliano’s method.

Castigliano’s second theorem provides a relationship between displacements, forces, and the strain energy in a linearly elastic system. The strain energy, U, is the potential energy of the system due to its deformation. For a generic elastic body, like this,

Elastic potato

Castigliano’s second theorem can be as:

If an elastic system is mounted so that rigid-body displacements of the entire system are impossible and certain external point forces P1, P2, … act on the system, in addition to distributed loads and thermal strains, the displacement component δi of the point of application of force Pi in the direction of force Pi is determined by the equation

δi=UPi

I took this from Henry Langhaar’s Energy Methods in Applied Mechanics, a favorite text of mine. I’ve altered the variable names to avoid conflicting with some variable names we’ve used in previous posts in this series.

Looking this over and thinking about it in terms of our problem,

Simply supported beam with uniform load

two concerns come to mind:

  1. Our problem doesn’t have a point force at the center of the beam, where we want to determine the deflection.
  2. Even if we did have a point force at the center of the beam, how do we write an expression for the strain energy in terms of that force?

Let’s tackle the second concern first. Recall from our two Rayleigh-Ritz posts (Fourier and polynomial) that the potential energy associated with beam bending is

0L12EI(y)2dx

This is the strain energy, U. Unfortunately, it’s written in terms of deflection, not force, but we can fix that. We saw early on that the curvature, y, and the bending moment, M, are proportional:

M=EIyory=MEI

We substitute the second of these into the expression for U to get

U=0LM22EIdx

Because the bending moment can be written in terms of the applied loads, we’ve solved the second concern.

We solve the first concern with a trick. We can pretend there’s a load, P, at the center of the beam, and write out an expression for the bending moment, M, that includes P.

Simply supported beam with uniform and concentrated load

We then do the integration to get U, take its derivative with respect to P, and—here’s the trick—evaluate the derivative for P=0. This may seem like cheating, but it’s perfectly legitimate. Imagine P is very very small—it won’t contribute much to the deflection. So why not just make it zero?

Let’s go ahead and see what happens. The upward support reaction force at end of the beam above is (wL+P)/2, so we can draw a free-body diagram for the left portion of the beam like this:

Free-body diagram for Castigliano

where x<L/2. The bending moment in the left half of the beam is therefore

M=12[(wL+P)xwx2]

and its square is

M2=14[(wL+P)2x22w(wL+P)x3+w2x4]

This formula applies only in the left half of the beam, but because of symmetry we know that

U=0LM22EIdx=20L/2M22EIdx=0L/2M2EIdx

So

U=14EI0L/2[(wL+P)2x22w(wL+P)x3+w2x4]dx=14EI[13(wL+P)2x312w(wL+P)x4+15w2x5]0L/2=14EI[L324(wL+P)2L432w(wL+P)+L5160w2]

Taking the derivative with respect to P gives us

UP=14EI[L312(wL+P)L432w]

Therefore the downward deflection at the center of the beam is

δ=UP|P=0=wL448EIwL4128EI=5wL4384EI

as we’ve now seen many times.

Since this is based on Castigliano’s second theorem, you may be wondering what kinds of problems we use his first theorem to solve. It’s a good question, but one I can’t answer. In the 45 years since I first learned of Castigliano’s theorems, I’ve never used his first one.