Understanding what NetGen is doing!

When using FreeCAD 0.17 it is not that difficult to make a simulation of something pressuring a cube that is fixed to the ground. I’t a simple case where we follow the typical simulation procedure… Make drawing, define boundary conditions, generate mesh, solve and get results.

Simulation Process Steps
Simulation Process Steps

I can understand that I need boundary conditions since these are the essence of the problem. But the mesh seems out of place. Or not… It is the building block of the finite element method where we describe a complex geometry into multiple small connected elements. The mesher will transform the solid in building blocks that represent the geometry we wish to study. We can use a mirrad of element types (block shapes) and the mathematical formulations they hold. NetGen should not be any different it can have tetrahedrons, hexahedrons, … and each specific nodes required…

But a question which came to my mind was how would the element shape originator and the element formulation originator work. Are they the same?Why I make this question? In freeCAD I generate a mesh with NetGen or Gmsh, however I solve it with Calculix and not the solver from NGSolve or Gmsh.

Before attempting a solution, a mesh is generated and thus the shape of the elements is defined and from the images they look tetrahedrons.

Paper Clip tetrahedron Mesh
Paper Clip tetrahedron Mesh
Paper Clip Circular section Mesh
Paper Clip Circular section Mesh

I found an answer to my questions in a paper written by one of the Authors of NetGen: Joachim Schöberl: An advancing front 2D/3D-mesh generator based on abstract rules.

The mesh is generated according to a specific sequence from a geometry. First special points are gathered on the corners of the target geometry, next edges have to be detected to then create a surface mesh that follows specific element shape rules; the final step is to create a volume mesh from the surface mesh and another specific rules necessary for element quality.

The output should be a cloud of nodes, elements and their connectivity to each other. From the paper it looks like it is purely connected to the geometry but then what connects to the differential equations we wish to solve? For instance how do I say a line is a beam? Or how do I create a body shell which basically is a shell with thickness.

Here comes the solver part. In case of calculix with NetGen used in this example all tetrahedrons should have a specific element type which will define the physics involved being structural, thermal, electric. In this case structural.

Indirectly the program gives a hint of the element type suggesting structural loads and forces in the boundary conditions or for a thermal case (not this one), the appropriate thermal boundary conditions.

As an example I am picking a Cube made in FreeCAD using NetGen mesher and Calculix solver. Before entering the solving process I generated the .inp file which contains the mesh, the boundary conditions and the material data to perform the simulation.

Somewhere in the .inp file you will find the following:

TYPE=C3D10

The element type name used is C3D10, as expected expected since normally the solver has a ton of element types, and there should also be in the Calculix solver a similar variety  like for instance Ansys with Solid187 (https://www.sharcnet.ca/Software/Ansys/17.0/en-us/help/ans_elem/Hlp_E_SOLID187.html), Solid186 or Solid226; for a list see this reference:

https://www.sharcnet.ca/Software/Ansys/17.0/en-us/help/ans_elem/

More details about C3D10 element can be seen in:

http://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node33.html

And a list off all Calculix elements can be seen here:

http://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node25.html

 

** written by FreeCAD inp file writer for CalculiX,Abaqus meshes
** Nodes
*Node, NSET=Nall
1, 10, 10, 0
2, 10, 10, 10
3, 10, 0, 0
4, 10, 0, 10
5, 0, 10, 0
6, 0, 10, 10
7, 0, 0, 0
8, 0, 0, 10
9, 0, 5, 10
10, 5, 10, 10
11, 10, 5, 10
12, 5, 0, 10
13, 0, 5, 0
14, 5, 10, 0
15, 10, 5, 0
16, 5, 0, 0
17, 10, 10, 5
18, 0, 10, 5
19, 10, 0, 5
20, 0, 0, 5
21, 5, 5, 10
22, 5, 5, 0
23, 5, 10, 5
24, 5, 0, 5
25, 10, 5, 5
26, 0, 5, 5
** Volume elements
<span style="color: #0000ff;"><strong>*Element, TYPE=C3D10, ELSET=Evolumes</strong></span>
25, 7, 4, 6, 8, 24, 21, 26, 20, 12, 9
26, 7, 6, 4, 1, 26, 21, 24, 22, 23, 25
27, 5, 7, 1, 6, 13, 22, 14, 18, 26, 23
28, 4, 7, 1, 3, 24, 22, 25, 19, 16, 15
29, 4, 2, 1, 6, 11, 17, 25, 21, 10, 23
** Define element set Eall
*ELSET, ELSET=Eall
Evolumes
***********************************************************
** Element sets for materials and FEM element type (solid, shell, beam, fluid)
** written by write_element_sets_material_and_femelement_type function
*ELSET,ELSET=SolidMaterialSolid
Evolumes
***********************************************************
** Node sets for fixed constraint
** written by write_node_sets_constraints_fixed function
** FemConstraintFixed
*NSET,NSET=FemConstraintFixed
3,
4,
7,
8,
12,
16,
19,
20,
24,
***********************************************************
** Materials
** written by write_materials function
** Young's modulus unit is MPa = N/mm2
** FreeCAD material name: 1C60
** SolidMaterial
*MATERIAL, NAME=SolidMaterial
*ELASTIC
210000, 0.300
***********************************************************
** Sections
** written by write_femelementsets function
*SOLID SECTION, ELSET=SolidMaterialSolid, MATERIAL=SolidMaterial
***********************************************************
** At least one step is needed to run an CalculiX analysis of FreeCAD
** written by write_step_begin function
*STEP
*STATIC
***********************************************************
** Fixed Constraints
** written by write_constraints_fixed function
** FemConstraintFixed
*BOUNDARY
FemConstraintFixed,1
FemConstraintFixed,2
FemConstraintFixed,3
***********************************************************
** Node loads Constraints
** written by write_constraints_force function
*CLOAD
** FemConstraintForce
** node loads on shape: Box:Face4
1,2,-0.0000000000000E+00
2,2,-0.0000000000000E+00
5,2,-0.0000000000000E+00
6,2,-0.0000000000000E+00
10,2,-1.6666666666667E+02
14,2,-1.6666666666667E+02
17,2,-1.6666666666667E+02
18,2,-1.6666666666667E+02
23,2,-3.3333333333333E+02
***********************************************************
** Outputs --&amp;gt; frd file
** written by write_outputs_types function
*NODE FILE
U
*EL FILE
S, E
** outputs --&amp;gt; dat file
*NODE PRINT , NSET=Nall
U
*EL PRINT , ELSET=Eall
S
***********************************************************
** written by write_step_end function
*END STEP
***********************************************************
** CalculiX Input file
** written by write_footer function
** written by --&amp;gt; FreeCAD 0.17.12321 (Git)
** written on --&amp;gt; Sun Oct 8 21:43:43 2017
** file name --&amp;gt;
** analysis name --&amp;gt; Analysis
**
**
**
** Units
**
** Geometry (mesh data) --&amp;gt; mm
** Materials (Young's modulus) --&amp;gt; N/mm2 = MPa
** Loads (nodal loads) --&amp;gt; N
**

 

The Paper Clip

A paper clip as simple as it is. It always is handy when we need to join paper sheets together. A soft metal wire is bended in such a way that when we insert a few paper sheets in-between the bended profile, it will act as a spring load pressing the paper sheets together.

Paper Clip
Paper Clip

The mesher used in this case was NetGen with standard options:

Max Size: 1000
Growth Rate: 0.3
Nbr. Segments per Edge: 1
Nbr. Segments per Radius: 2

Resulting Mesh:

Node count: 80181
Triangle count: 25166
Tetrahedron count: 39417

Paper Clip tetrahedron Mesh
Paper Clip tetrahedron Mesh

But what did NetGen actually do? It did a mesh but many questions should now start to appear. Some questions came to my mind, like what elements did it used, what is the general quality of the elements like for instance:

  • Element quality criteria;
  • Aspect ratio
  • Jacobean
  • Warping factor
  • Maximum corner angle
  • Skewness
  • Orthogonal quality

Also I wonder what element types are used and where can I find documentation. I will try to understand this better in my next posts…

The red candle cup.

Day and night is present in our lives since the moment of our birth. Place a lit candle in this cup and a world of reflection opens, giving a fantastic ambience. Light is all about sensation, math, physics and some mysticism at the mixture. Light has been studied by many geniuses from Albert Einstein to Richard Feynman giving birth to fantastic ideias about the universe around us, but a side from the speed of light, it’s diffraction  the common human knows little about it.

Many activities explore light, modelled it and play with it. Photography explores light to an outstanding level. Texture, highlights, color, black and white, grey toned images are just a few examples.

Unfortunately I don’t have much answers regarding light simulation and what can be done with it. Interestingly games play with light from rendering shadows, light glare, and more effects. Photo editing software also create lens flare effects, maybe these can be considered to be some sort of light modelling.

Either by exploring light on a nice summer day or by taking that photo at sun down, light remains the source of inspirations that holds still many secrets to tell.

Red_and_silver_Squared_Glued_Glass_Cup
Red and silver Squared Glued Glass Cup

Fibonacci Spiral in Python

I always like to play around with my models and see what happens next. Don’t you? Problem is that most of the time we don’t know how to do it and with what resources… Let me show you how to model a Fibonacci spiral easy.

First we will use Python. Don’t worry. Install Python or just install Anaconda and open Spider.

Once in spider paste the following code to the editor and hit play.

The result should be something like this:

Fibonacci Python Generated Spiral
Generating a Fibonacci spiral.

There is a nice page that explains how the python script works and gives a deeper insight into the Fibonacci series:

Week 1: Fibonacci Sequence and the Golden Ratio

# chromoSpirals.py
# ----------------
# Code written by Peter Derlien, University of Sheffield, March 2013
# Draws spiralling patterns of circles using the Golden Angle.
# ----------------
# Import from the numpy and matplotlib packages.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.collections import PatchCollection
import matplotlib.patches as mpatches
ox=0.5; oy=0.4 # centre of plot
ndiscs=300
ndiscs=input('No. of discs (e.g. 300)? ')
ndiscs=int(ndiscs)
ncols=input('no. of colours (1 to 34)? ')
ncols=int(ncols)
offset=0.0
offset=input('offset (in radians) from golden angle? ')
offset = float(offset)
tau=(1+5**0.5)/2.0 # golden ratio approx = 1.618033989
#(2-tau)*2*np.pi is golden angle = c. 2.39996323 radians, or c. 137.5 degrees
inc = (2-tau)*2*np.pi + offset
theta=0
k=0.1 # scale factor
drad=k*(1+5**0.5)/4.0 # radius of each disc
minv=maxv=0 # minv and maxv will be used later to display inputs chosen
# now collect in list 'patches' the locations of all the discs
patches = []
for j in range(1,ndiscs+1):
r = k*j**0.5
theta += inc
x = ox + r*np.cos(theta)
y = oy + r*np.sin(theta)
if y &gt; maxv:
maxv=y
elif y &lt; minv:
minv=y
disc = mpatches.Circle((x,y),drad)
patches.append(disc)
# start building the plot
fig = plt.figure()
ax = plt.axes([0,0,1,1])
# create text to show which inputs the user has chosen
font = "sans-serif"
maxv=maxv*0.95
nd = 'ndiscs: '+ str(ndiscs)
plt.text(minv, maxv, nd, ha="center",family=font, size=14)
setting = 'angle offset: '+ str(offset)
plt.text(minv, minv, setting, ha="center",family=font, size=14)
nc = 'ncols: '+ str(ncols)
plt.text(maxv, maxv, nc, ha="left",family=font, size=14)
# build colour cycle, using a number between 0 and 100 for each colour
colcycle=[]
s=100/ncols
for j in range(ndiscs):
colcycle.append((j%ncols)*s)
# bring together the information for locations and colours of discs
collection = PatchCollection(patches, cmap=matplotlib.cm.jet, alpha=1.0)
collection.set_array(np.array(colcycle))
ax.add_collection(collection)
ax.set_xticks([]); ax.set_yticks([]) # suppress display of axes
plt.axis('equal')
plt.show() # display the plot we have built 

Natures numbers – Fibonacci

It is really beautiful when we can see math in the world around us! How can a flower find the most optimised structure to balance growth, sunlight and moisture in each seasonal energy burst? I guess millions and millions of year and we get the golden ratio in the fibonacci series.

Fibonacci Spiralled Flower
A beautiful image portraying Fibonacci Structures

 

Mathematical Spirals

Today we use super computers and special algorithms for complex problems, however in the past mathematicians could model problems and solve them with very little resources. It is always good to get to the basics and see art in the form of math.

Leonardo of Pisa nicknamed “Fibonacci”, brought state of the art Arab mathematics to medieval Europe. One of his contribution was his book Liber Abacci where it described a series that we commonly see in nature. It’s all about spirals but they follow a rule…

The Fibonacci series can be seen as a numeric sequence starting with two ones where each subsequent number is equal to the sum of the preceding two numbers: 1, 1, 2, 3, 5, 8, 13, and so on.

FibonacciSpiral

Although this number series was known in India previously, Leonardo helped spread this knowledge worldwide.

The most interesting thing is that  he stumbled upon it trying to solve a rabbit problem! I think the book proposed the following:

 ‘A man put a pair of rabbits in a place surrounded on all sides by a wall. How many pairs of rabbits are produced from that pair in a year, if it is supposed that every month each pair produces a new pair, which from the second month onwards becomes productive?’

 

 

Playing a bit with the results in the Brompton bike handlebar case

On my previous post I got to discover that the new design for the 2017 Brompton bike model is much more robust. Didn’t get to understand if it was designed on purpose or was a technical difficulty. My hunch was that they enjoyed the design and it made construction simple. However with experience they concluded that it was a problematic design and went for a new design. They could have just reduced the central bar length but instead almost removed it to a more challenging design at least for manufacturing.

Bent and unbent version
Bent and unbent version

FreeCAD produces nice images and this one is an undeformed vs deformed shape. Here I think it is visible that near the stem is where most of the stresses will be generated. So if you have a handlebar from a bellow 2017 model take special care to not pressure in an exaggerated manner the handlebar down.

Special thanks: I would like to thank the FreeCAD forum team and in particular Bernd for providing troubleshooting for FreeCAD  and showing how to do the deformed and undeformed image.

The Brompton Handlebar Case II

While browsing the Internet I found some past references to Brompton Handlebar cracks near the bending elbows. I noticed my colleague Brompton had a different design and thus concluded that probably the Brompton bike team, decided to turn the  handlebar more robust.

Brompton handlebar
Brompton handlebar

Check out my previous post to know what I am talking about. Click here

I notice that I make some force to the ground direction when pedalling and curving, special when the road has bumps.

HandleBar model constraints
HandleBar model constraints

In wich elbow will the stress higher and the cracks be formed? I will be continuing on this subject, so please stay tuned.

The bridge builder game effect

Going back in time, I just remembered about a game that gave me quite a good understanding of stress distribution. Bridge builder. Using simple truss combinations to allow passage of vehicles through it. Some times the most chalenging designs required the use of circular shapes to transfer stress from one point to another, distributing it through as many bars as possible…

Bridge Builder
Bridge Builder

This sun lounger supporting legs are made of a single component and has one main function, to support the weight of the occupier and a bonus function to support the arms. Stresses are evenly distributed throughout the structure and the design feels more organic.

Stress Distribution
Stress Distribution