Another interim report – I found FreeCAD's FEM Workbench has been updated significantly since I last used it, and took a while to get anywhere. Pretty sure I'm not doing it right! Finite Element Method is the application of a 'Finite Element Analysis', more detail about the theory on Wikipedia.
Anyway, in FreeCAD's Part Design Workbench I sketched a 20mm square and extruded it by 100mm to create a solid 3D body. Then I drew a crude spanner outline on one end, and extruded that. I made a number of cuts in the 'spanner' so that it only makes contact with the square rod in realistic places, thanks Clive. (Cuts are made by sketching rectangles and applying the pocket tool.)

Once the rod and spanner combination existed, I switched to the FEM Workplace and:
- created a Solver
- Defined the Solid Material as Aluminium 6061-T6 (only because it happens to be the first metal listed in the database.) The database defines the properties of many common materials: Density, Young's Modulus, Poisson Ratio, and Thermal Properties, or you can enter your own. These measures are used by the Solver to calculate bending, strength and other structural behaviours.
- Added a 'mesh' to the model. The computer fills the model with small triangles. These are used by the solver to break the analysis into manageable chunks. By solving stress in one triangle, the results can be applied to it's neighbours until the whole object is solved. Note the computer doesn't look up the formula for a bending a beam and apply that, rather it generalises the problem so that any shape can be solved. The finer the mesh the more accurate the result. As there's a huge amount of complicated number crunching, a computer is essential. Too much work for humans, and tedious.

- Added a fixed constraint to the far end of the square rod. This is equivalent to clamping it in a vice.
- Added a sideways force to the end of the spanner handle such that the rod will try to rotate. As the other end is constrained, the rod will twist. How much will the rod twist, what will break first, and which parts are over-strong?
At this point the Solver has enough information about the object to perform an analysis. Here's the whole Model as listed in the 'Combo View':

Double clicking the Solver 'CalculiXccxtTools' line opens a dialogue. Clicking the 'Write .inp File' button saves the specification and a 'Run' button appears. Pressing this adds the CalculiX_static_results entry, and the Results_mesh. Double Clicking 'CalculiX_static_results' opens another dialogue that can be used to select various displays, in this example X-displacement:

The min, avg, and Max numbers aren't correct because I messed with the model but the colour map correctly shows the spanner in tension and compression on opposite sides and also that the rod is twisting, mostly at the top and not at all at the bottom. Hardly surprising, I drew a weedy spanner and a hefty rod.
Not worked out how to tabulate results. There are a lot of other features I haven't explored and I don't properly understand the ones I've got working either! Dynamic forces can be analysed, also heat, fluid flow in a pipe, and something to do with electrostatics. No ideas what the new set of tools for post-processing results do!
I can't answer the original square vs hex head question, but how it might be done with FEM is taking shape.
I've also found that Machinery's Handbook has a couple of pages on shafting of various shapes. There are formula for all the common shapes. A mathematician might 'see' which, if any, are stronger just by looking at the formula. Not me, I'd have to feed them all with numbers and compare the results. Another nice job for a computer!
More if I get the time. Too many diversions.
Dave
Edited By SillyOldDuffer on 10/06/2019 17:34:22