Revisiting Castigliano with SciPy
January 6, 2013 at 2:16 PM by Dr. Drang
On Friday, a colleague asked me if I had a quick solution for determining the spring stiffness of a tapered leaf spring. Yup. This may be the first time I’ve been able to use an old blog post directly for work.
But as I read through the solution, I realized I’d done the numerical integration in Octave, whereas now I’d prefer to do it in Python and SciPy. The dilemma: do the problem immediately in Octave, where I had a complete solution, or do it in SciPy, where I’d have to search to find and learn the equivalent functionality? I decided to go the latter route because:
- My colleague need a quick answer, but not an instantaneous one.
- I figured numerical integration would have to be one of the most prominent SciPy features—easy to find and use.
- It’s in doing practical problems like this that I really learn.
I’m not going to rewrite the logic of that older post. In a nutshell, the problem of a tapered leaf spring is tailor-made for Castigliano’s Second Theorem, which states that the deflection of point upon which a load is acting is equal to the derivative of the complementary strain energy of the structure with respect to that load.
The tapered leaf spring can be reduced to a cantilever beam that looks like this:
The complementary strain energy is
where is the bending moment function, is Young’s modulus for the material, and is the moment of inertia for the cross-section,
where is the (constant) width of the beam and is the rate at which the thickness increases.
Putting these expressions together and taking the derivative, we get
and the equivalent spring stiffness is
OK, now it’s time to move onto SciPy to handle that integral. I could, of course do the integration analytically, but that would involve integration by parts and I’d almost certainly end up with a big algebraic mess. And I need a number for an answer, so I might as well go numerical now as later.
SciPy’s integration functions are described here. There are basically two sets: one set for which the functional form of the integrand is given, and another set for which the integrand is known only as a set of pairs. We’re going to use the function quad from the former set. In its simplest form, quad takes three arguments: the function to be integrated, the upper bound, and the lower bound. Here’s a short Python script that calculates the stiffness for a steel leaf spring with , , , and .
python:
1: #!/usr/bin/python
2:
3: from scipy import integrate
4:
5: # Given parameters
6: E = 29e6 # Young's modulus, psi
7: L = 25 # length, in
8: b = 3 # width, in
9: t0 = .5 # tip thickness, in
10: t1 = 1 # end thickness, in
11:
12: # Derived parameter
13: alpha = (t1 - t0)/L
14:
15: # Integrand
16: def integrand(x):
17: return x**2/(t0 + alpha*x)**3
18:
19: # Stiffness
20: integral = integrate.quad(integrand, 0, L)
21: print integral
22: print E*b/12/integral[0]
The quad function returns a tuple with the answer and the error, so I have that printed out on a line by itself so we can check the accuracy. The results are
(8518.397569993163, 9.457321115117394e-11)
851.099040686
If you look back at the older post, you’ll see that this is the same answer we got using Octave, which isn’t surprising—this isn’t an especially tricky integral to perform numerically.
You’ll note that the first argument to quad is itself a function that takes only one argument. What if you wanted to define the moment of inertia parametrically, like this:
def I(x, b=3, t0=.5, t1=1):
alpha = (t1 - t0)/L
return b*(t0 + alpha*x)**3/12
You’d do this if, for example, you had several different leaf spring designs to consider, and you didn’t want to keep rewriting the integrand function over and over again.
The solution is to use the partial function from Python’s standard functools library, which allows you to define functions in terms of other functions. Here’s an example:
python:
1: #!/usr/bin/python
2:
3: from scipy import integrate
4: from functools import partial
5:
6: # Design parameters
7: E = 29e6 # Always steel
8: L = 25 # Always this long
9: params = []
10: params.append(dict(b = 3, t0 = .5, t1 = 1))
11: params.append(dict(b = 3, t0 = .5, t1 = 1.25))
12: params.append(dict(b = 3, t0 = .75, t1 = 1))
13:
14: # Moment of inertia function
15: def I(x, b=3, t0=.5, t1=1):
16: alpha = (t1 - t0)/L
17: return b*(t0 + alpha*x)**3/12
18:
19: # Integral
20: def integral(p):
21: myI = partial(I, **p)
22: def integrand(x):
23: return x**2/myI(x)
24: return integrate.quad(integrand, 0, L)
25:
26: # Stiffnesses
27: for i,p in enumerate(params):
28: print "k%d = %f" % (i, E/integral(p)[0])
In this example, I’ve assumed that the leaf spring will always be made of steel and that its length is set. The design parameters that can be altered are the cross-sectional dimensions.
Lines 9-12 set up a list of dictionaries of design parameters. The integral definition in Lines 20-24 takes a dictionary as its argument, defines a one-argument moment of inertia function based on that dictionary, and performs the integration for that moment of inertia function. The loop in Lines 27-28 then goes through the list of parameters and prints out the equivalent spring stiffness for each design. This is not only more compact than writing a new integrand function for each design, it’s also more flexible when the time comes to investigate new designs. This flexibility is where Python and SciPy have a distinct advantage over Octave.
By the way, the output of this script is
k0 = 851.099041
k1 = 1436.267876
k2 = 1127.163920
which shows that if you’re going to add a little thickness to the beam, you’re better off adding it to the fixed end than to the free end. Which you probably would’ve guessed anyway, but here you get to see how much more effective it is.
One last thing: the quad function gets its name from quadrature, a commonly used synonym for integration, especially numerical integration. The quad in quadrature refers not to the number four (not directly, anyway), but to a square. In classic compass-and-straightedge geometry, quadrature meant squaring, the process of constructing a square of the same area as some other shape. That definition morphed into calculating an area without constructing a square, and from there it was a short jump to integration.
If you follow Gus Mueller’s blog, you might remember a short post he did a few months ago about quadratic and cubic Bézier curves. That’s another case where quad comes from square rather than four. Of course, squares are associated with quad because squares have four sides, so ultimately quadrature and quadratic equations get their names from the number four. But it’s a pretty roundabout trip.
