Splines and Curves of Least Energy

[tropostudio]

I am using Excel Solver in the Recalc button. That is used when you change the approximation interval. 'n' in cell C4 can vary from 8-2048 in both the Euler Spiral or Horn Curve worksheets. I'll start with the Euler Spiral worksheet because it was easier to put together and the first one I did.
Orange cells are for input. gey cells with orange type are calculations. Cells below the calcs are arrays that are resized and recalculated when you use the Re-Calc button. Total curve length, S, andMaximum Curvature, kmax can be adjusted, but if you are looking for the curve that ranges through x=1 to x=o, you can leave those alone. I usually leave C7, kmin, =0. That is at x=1, and asymptotic to a vertical line

It is possible to do an ad-hoc Solver routine that sets J10 to 1 by varying C6, but it really works better to 'bit bang' C6 until you get a value at J10 close to 1. The yellow highlighting indicates we are looking at varies from theta = o - pi/2 rad. O10 is the scaling factor that normalizes x max to x=1. R10 is the normalized total energy in the curve. You can look through the cell formulas and VBA macros to figure out the est. The Re-Calc Curve button cycles through updating the values and re-drawing the charts.

I'll get to the Horn Curve later, but it's similar. The geometry is just more work and calculation time than it takes to graph an Euler Spiral.

 

[NirvanaManana]

Hi Chris,
I was wondering if you developed your solver any further? This is really interesting work.

 

[tropostudio]

Hi Chris,
I was wondering if you developed your solver any further? This is really interesting work.

Funny this popped up. I haven't had time to focus on this project until a couple days ago. The solver/spreadhseet stuff hasn't been touched, but I know what needs to be done next. As a visual thinker and builder, I wanted an apparatus to facilate the beam-bending experiments that started me down this rabbit hole. To that end, I'm just about done machining a set of six 'Suction Cup Spline Ducks' that can be relocated around a 4'x 8' table top:



The CAD image shows a bottom view, with one unit stuck to the screen. The pultruded fiberglass rods I've been using as 'uniform beams' slide between the pair of v-wheels, which in turn are free to rotate around a sealed bearing. The rubber-tipped toggle clamp legs can be adjusted to tune the holding force of the suction cup (which is an OEM unit used for automobile glass resin injectors).

I'll grid out the table surface, color-code the acorn nut atop the rotater, and use a contrating color for the table surface (battens are white). A digital camera will be centered overhead. The grid will be calibrated to give accurate point mapping through OpenCV. Spline interpolation will be part of this. This is half the fun of any project: as an exhibit builder by trade, I can't 'just do the maths!'

Regarding 'the maths': From a practical standpoint, interpolating a symmetrical curve to a prescribed precision is pretty trivial. You have a start point, a midpoint, and a mirrored endpoint. Curvature is 0 at each end and equal at the midpoint.

An assymetrical curve is tougher. An iterative/convergence process is probably in order. It ceratinly is for a person who's mediocre at math and finds elliptic integrals (basis of the 'Horne Curve') challenging! My thinking: use the ratio of the length of two chords defined by your three points to start. Curvature at each end will be zero. Scale one curve section against the other based on chord length. Compare curvatures at the common point. Make your next 'scaling' guess based on how far the two curvatures differ and do it again. Lather, rinse, and repeat until you get convergence to a desired precision. I'm sure it can be done analytically. Whether that leads to a faster solution on a computer vs using an iterative, numerical approximation? Dunno...

 

[NirvanaManana]

Ah, I was hoping you would have done the work on the asymmetrical curve and that I wouldn't have to!

I'll be defining asymmetric Curves of least energy based on three control points. Looks like I'll just have to get stuck into it. With these problems in Excel the hardest thing is to get your head around defining the problem and how to work out a way of solving it. You have done that already.

 

[philSweet]

The superposition principle is your friend.

 

[tropostudio]

Nirvana Manana - The race is on! I'm pretty slow out of the gate

philSweet - Did you have something in mind that could help here? I can see that superposition applies to combined loads in simple beam bendingt. I'm not sure how it helps when we are trying to model an accurate 'curve of minmum energy through inflection points.' Cubic splines are based on simple beam theory, and is the basis for 'natural cubic' splines in a lot of software. The fact that it does not model real-world bending for any but small deflections is what started this project.

 

[DCockey]

Superpostion requires a linear system. Simple beam theory is linear. The models/equations considered in this thread are nonlinear so superpostion does not work.

 

[tropostudio]

I think the 1962 book 'Flexible Bars' by R. Fritsch-Fay as being the full-on 'mechanics of materials' background for this. It can be read online at Scribd: Flexible Bars Frisch-Fay 1962 | PDF https://www.scribd.com/document/317704558/Flexible-Bars-Frisch-fay-1962 An online PDf can also be found here: https://bigoni.dicam.unitn.it/varie/flexible_bars_frisch-fay_1962.pdf.

From page 41: The portion of the curve from A-B is exactly the curve described by Horne, and the condition shown in my photos.

From page 81: An example backing up his diagrams and math, for a symmetrical condition. The shallower curves can all be derived from a segment of a steeper curve that undergoes scaling and rotation. 'Superimposing' is useful, but the principle of superposition does not apply when deflections are not a linear function of bending moment. Fritsch Fay points out the issue on page 1 of the book.

 

[tropostudio]

Here is a Masters thesis on compound bow limb design that may be applicable here :

Optimal design of the limb in compound bows

There is a link to the author's thesis on the site. He mentions Fay's work, and how difficult the math gets. He uses a 4'th order Runge-Kutta approximation method to calculate limb deflection. Much simpler, and any decent numerical methods book will explain Runge-Kutta approximation.

 

[NirvanaManana]

RK4 brings me back - when I used it first in college we had to do hand calculations. It leaves you with an appreciation for Excel. Used it for solving the motion of a double pendulum. It's very powerful.

I've done the Excel model for a curve of least energy for a three point constraint (like a peg board). It can be easily adapted to solve for a curve constrained by two points one of which is also has a fixed angle (like a diving board, with the board protruding at any angle). These problems are more straightforward than what Tropo is solving for above - which is solving for three constraints, but only one point and two angles. I didn't have to do the macros that Chris employed. I'm away for a week but can post the file if anyone's interested.

 

[tropostudio]

I am interested. Please post your file when you are back!

 

[tropostudio]

Things take time when one still has to work , and is also easily distracted Neutral grey p-lam surface holds the 'ducks' , until the vacuuum breaks, probably because it is matte laminate and has surface roughness. Anyway, I now see a digicam overhead for image grabs and ore processsing. Or a pico-projector to project grids or calclulated reference curves for comparison. Smaller diameter FG splines stress the ducks less, but also are less sensitive to bearing friction - less precise. Manufacturing tolerances in the 'uniform splines' of pultruded FG rod is apparent, or that could be due to slight variations in bearing friction in my 'ducks.' All-in-all, it works as designed. I'm happy to get back to math that I can verify with a physical model.

Still quite certain traditional M/EI simplifications, or Cornu/Euler 'curve of least energy' based on curvature based linearly on arc length along the curve won't match. Even the longest, leanest hullform will have curvatures that make for gross errors.

 

[tropostudio]

Progress has been made on a method for calculating a curve of least energy that has the following conditions: Curvature at endpoints = 0 and a single intermediate point that is [obviously] G2 continuous. The 'true' curve of least energy will be a blend of 2 'Horn Curve' sections - not of 2 Euler Spiral sections. But, the method for belneding the 2 curve sections works identically for either the Horn curve or an Euler [Cornu] spiral.

The Excel spreadsheet needs some VBA macro work and clean-up befoe sharing, but I figured I'd post some screenshots that might illuminate recent work:

A rehash, but here's a shot of an Euler spiral section vs a comparable 'BKP Horn' curve section. Recall from the previpous posts that the Horn curve has less energy than a corresponding Euler spiral. Experiments on a workshop table with pultruded GRP rods are what led me down this rabbit hole.

The latest spreadsheet uses a slightly different method for finding points on the Horn and Euler curves than I started with, but the end result is the same. Looking at graphs of an Euler Spiral and its 'curvature curve' vs a Horn Curve and its curvature curve gives a sense of why the Horn Curve would indeed have lower energy than an Euler Spiral:

Recall that neither curve can be rpresented as an algebraic function or polynomial. They are transcendental functions and have required the use of numerical approximation. In either case in my case, I'm using a user-selected subdivision ('n') range of 32-4096 to claculate points on a 'normed curve' section from theta = 0 to theta = Pi/2.

The method for blending two curves is illustrated here for two Euler Spiral sections, but will be identical when applied to two Horn Curve sections. First, etsablish 3 control points (nodes) that will insure 2 sections are blended that fit within the bounds of our 'normed' curve section (keeps the problem 'simple enough'). After much stumbling around, I observed a tangent condition that 'flip-flopped' supplemantary interior angles gave a decent first-order approximation to a final tangent condition at the shared node (point P2):

The next thing to do was to find portions of the 'normed' curve that corresponded to a scaled chord length, curvature, and tangent angle that enforced G2 continuity at the shared node (point P2). I start with looking up the points immediately above and below the 'final point' and getting a more accuarte location for the final point based on interpolation between the upper and lower lookups. The interpolated point can simply be a linear fit between the upper and lower points, or it can be located more accurately by usng an interpolated chord length and angular deviation of that chord with respect to the x-axis. These two estimates will converge as interval count increases [and interval subdivisions become smaller]:

Here are two Euler Spiral curve sections that will be blended shown at both 'normed' scale and 'final' scale ( final chord lengths and matching P2 curvature):

And a final chart with the pair of blended curves, along with the preceding affine transformations for each section:

I look forward to empirically proving this with Horn Curves for 3 points on my fixture table...

...and then to following this with a method with empirical proof for 4+ nodes:

 

[philSweet]

The problem with multiple nodes can be set up in excel for the plug-in XLAM solver. I use the gradient solver, but have played with setting the problem up for some of the other methods in XLAM as well. I did extensive work on this curve because it also happens to be the same curve that describes ideal midship properties of sailboats based on simple RM maximization. For instance - the ideal heeled section shape that gives maximum RM at a particular heel angle and girth length will have a curvature that varies linearly with the lateral distance from the VCG of the heeled boat. Using the half girth as the reference dimension, it quickly gets complicated when the section centerline is immersed and the section has to be symmetrical about the unheeled vertical axis. This constraint is not unlike adding more constraints to your curve with new ducks.

This isn't much of a challenge computationally for the computer, so you don't need to worry at all about execution efficiency. The trick is entirely about presenting the problem in a way that XLAM can digest it and converge on a solution over the entire design space. Note that as you approach the lowest energy shape, small changes in shape will have almost no change in energy, so the incremental effects that are significant will decrease by five or six orders of magnitude as the solution is approached. Baking that in to the iteration scheme takes planning.

I represent the curve by a set of fixed-length links and operate on the bend angles between links as the primary variable set in the optimization scheme. I would start with a known solution, such as your first post curve, and place your ducks where you want them, plus a phantom duck on the middle of the curve. The initial measure of merit for the solver will consist of the energy integral of the curve plus a penalty for the distance from each duck. The initial penalty, say the sum of the square of the distance from each duck to its nearest link, should be low, except for the phantom duck, which will be high; and then run the optimizer to find the solution. Then increase the distance penalty and the solver will move the curves closer to the ducks. Decrease the distance penalty from the phantom duck at the same time. This is a problem presentation that a gradient solver should like. The key is to keep all the terms that are contributing to the measure of merit about the same order of magnitude during the entire process. You want to set up, keep track of, and update separate increment values for each link angle and distance penalty so that any variable changed by it's dedicated increment will have about the same effect on the measure of merit. And as already noted, it will decrease by many orders of magnitude.

As for an iteration scheme, you could randomly pick 5 link angles. For each one, apply it's increment and compute the measure of merit. If the value is negative, repeat with negative increment. Compare the five, pick the biggest winner, and apply this to the curve. Keep track of how often each of the variables wins this race, whether it was an increment or decrement, and adjust the random generator and increment sizes to maintain balance. I plotted this data during the run.

As a final note, I did a pretty extensive error analysis on how to construct, manipulate, and revise the links in response to the evolving shape. It's not terribly important unless you want a really accurate estimate of the curve's length. Note that the energy curve will be really flat at the solution point, and it is an integral over length, so it is important that the links consistently represent the actual curve length of their segments. This can be fine tuned as you go or at the end.

 

[tropostudio]

@philSweet - thanks for your reply. It'll take me a bit to digest your use of Solver here. I'm familiar with it and used it in early versions of the spreadsheets on this subject. I used Goalseek for one operation on the two-piece curve blend because it was simple, and it 'usually' converges. But there are plenty of times where it, and even Solver for that matter, doesn't. For Goalseek, there's no fixing it. For Solver, I can probably set up the problem better (operator error and unfamiliarity). I often have better luck at reaching convergence by manually adjusting the variable values observing the approach to the goal, working from most to least significant digit. Not sure what the method is - seems too obvious to not have a million examples. It's not bisection, not Newton - Raphson...