A script for elastic bending (aka the elastica curve)

Hey all,

I thought I'd share this script I've been working on that solves the problem, "What shape does a long, narrow pole/beam/wire make when it undergoes large elastic (non-permanent) bending deformation?" I know that this can already be accomplished using the brilliant Kangaroo plugin, but I wanted a simpler and faster (yet still accurate) single component that could replicate this unique curve using a variety of inputs: the length of the rod/wire, the width/distance between the endpoints, the height of the bend, and the tangent angle at the start. I also wanted make the unknowns (such as height if only length and width are known) easily accessible for plugging into additional components.

The resulting script, being an all-in-one solution, is somewhat unwieldy, but it could easily be broken down into smaller components (custom .gha's which I don't have the ability to code). If someone wants to tackle this, please do! I'm not an expert coder by any means, and as this was only my second time diving into Grasshopper scripting, if the script seems somewhat strange, that's probably why. I did try to comment the code pretty well though. Here's the full description:

--------------------------------------------------

DESCRIPTION:
This beast creates the so-called 'elastica curve', the shape a long, thin rod or wire makes when it is bent elastically (i.e. not permanently). In this case, force is assumed to only be applied horizontally (which would be in line with the rod at rest) and both ends are assumed to be pinned or hinged meaning they are free to rotate (as opposed to clamped, when the end tangent angle is fixed, usually horizontally). An interesting finding is that it doesn't matter what the material or cross-sectional area is, as long as they're uniform along the entire length. Everything makes the same shape when bent as long as it doesn't cross the threshold from elastic to plastic (permanent) deformation (I don't bother to find that limit here, but can be found if the yield stress for a material is known).

Key to the formulas used in this script are elliptic integrals, specifically K(m), the complete elliptic integral of the first kind, and E(m), the complete elliptic integral of the second kind. There was a lot of confusion over the 'm' and 'k' parameters for these functions, as some people use them interchangeably, but they are not the same. m = k^2 (thus k = Sqrt(m)). I try to use the 'm' parameter exclusively to avoid this confusion. Note that there is a unique 'm' parameter for every configuration/shape of the elastica curve.

This script tries to find that unique 'm' parameter based on the inputs. The algorithm starts with a test version of m, evaluates an expression, say 2*E(m)/K(m)-1, then compares the result to what it should be (in this case, a known width/length ratio). Iterate until the correct m is found. Once we have m, we can then calculate all of the other unknowns, then find points that lie on that curve, then interpolate those points for the actual curve. You can also use Wolfram|Alpha as I did to find the m parameter based on the equations in this script (example here: http://tiny.cc/t4tpbx for when say width=45.2 and length=67.1).

Other notes:

* This script works with negative values for width, which will creat a self-intersecting curve (as it should). The curvature of the elastica starts to break down around m=0.95 (~154°), but this script will continue to work until M_MAX, m=0.993 (~169°). If you wish to ignore self-intersecting curves, set ignoreSelfIntersecting to True

* When the only known values are length and height, it is actually possible for certain ratios of height to length to have two valid m values (thus 2 possible widths and angles). This script will return them both.

* Only the first two valid parameters (of the required ones) will be used, meaning if all four are connected (length, width or a PtB, height, and angle), this script will only use length and width (or a PtB).

* Depending on the magnitude of your inputs (say if they're really small, like if length < 10), you might have to increase the constant ROUNDTO at the bottom

REFERENCES:
{1} "The elastic rod" by M.E. Pacheco Q. & E. Pina, http://www.scielo.org.mx/pdf/rmfe/v53n2/v53n2a8.pdf
{2} "An experiment in nonlinear beam theory" by A. Valiente, http://www.deepdyve.com/lp/doc/I3lwnxdfGz
{3} "Snap buckling, writhing and Loop formation In twisted rods" by V.G.A. GOSS, http://myweb.lsbu.ac.uk/~gossga/thesisFinal.pdf
{4} "Theory of Elastic Stability" by Stephen Timoshenko, http://www.scribd.com/doc/50402462/Timoshenko-Theory-of-Elastic-Sta... (start on p. 76)

INPUT:
PtA - First anchor point (required)
PtB - Second anchor point (optional, though 2 out of the 4--length, width, height, angle--need to be specified) [note that PtB can be the same as PtA (meaning width would be zero)] [also note that if a different width is additionally specified that's not equal to the distance between PtA and PtB, then the end point will not equal PtB anymore]
Pln - Plane of the bent rod/wire, which bends up in the +y direction. The line between PtA and PtB (if specified) must be parallel to the x-axis of this plane

** 2 of the following 4 need to be specified **
Len - Length of the rod/wire, which needs to be > 0
Wid - Width between the endpoints of the curve [note: if PtB is specified in addition, and distance between PtA and PtB <> width, the end point will be relocated
Ht - Height of the bent rod/wire (when negative, curve will bend downward, relative to the input plane, instead)
Ang - Inner departure angle or tangent angle (in radians) at the ends of the bent rod/wire. Set up so as width approaches length (thus height approaches zero), angle approaches zero

* Following variables only needed for optional calculating of bending force, not for shape of curve.
E - Young's modulus (modulus of elasticity) in GPa (=N/m^2) (material-specific. for example, 7075 aluminum is roughly 71.7 GPa)
I - Second moment of area (or area moment of inertia) in m^4 (cross-section-specific. for example, a hollow rod would have I = pi * (outer_diameter^4 - inner_diameter^4) / 32
Note: E*I is also known as flexural rigidity or bending stiffness

OUTPUT:
out - only for debugging messages
Pts - the list of points that approximate the shape of the elastica
Crv - the 3rd-degree curve interpolated from those points (with accurate start & end tangents)
L - the length of the rod/wire
W - the distance (width) between the endpoints of the rod/wire
H - the height of the bent rod/wire
A - the tangent angle at the (start) end of the rod/wire
F - the force needed to hold the rod/wire in a specific shape (based on the material properties & cross-section) **be sure your units for 'I' match your units for the rest of your inputs (length, width, etc.). Also note that the critical buckling load (force) that makes the rod/wire start to bend can be found at height=0

THANKS TO:
Mårten Nettelbladt (thegeometryofbending.blogspot.com)
Daniel Piker (Kangaroo plugin)
David Rutten (Grasshopper guru)
Euler & Bernoulli (the O.G.'s)

--------------------------------------------------

Edit: More on the math behind this here.

Cheers,

Will