Plant stem inspired, topology optimised infill structure for 3D printing

I am a regular user of 3D printing with plastic Fused Deposition model type printers such as the the Prusa printers, Up Minis. I was thinking about how to optimise the internal infill, for minimum weight while being strong enough. Infill is a necessity for these types of 3D printing because it allows the printers to print on overhangs reliably. So I turned to nature as inspiration.

I made a trip to the Imperial College Advanced Hackspace(ICAH) in white city on 19th August 2019 to use the microscopes there to look closely at the Buddleia plant, an invasive species that is widespread on train tracks in the UK. I intended to get transverse sections of the stems and view them under a microscope. However, I encountered difficulties slicing the stems to very fine thickness(~100 microns) as seen in stock photos because I did not have access to a microtome. However, it turned out to be a blessing in disguise though I could only get ~1000micron thick slices.

These are the pictures I took with the digital camera mounted on the microscope. All images captured to an SD card on the camera. The photos have all been edited for clarity with Photoshop to improve sharpness and contrast. No shape distortion carried out. The microscope was intended for Printed Circuit Board(PCB) soldering and assembly. Therefore, the microscope had a large focal range, allowing me to see depth. This was critical in being able to identify the cellular structure in 3D, that I had not seen in the extremely thin microtome sliced images seen in my previous analysis of stems.

In the thick slices, I found that the tubes are not that simple. They do not simply extend to the top like stright tubes. Instead, they are a cellular structure. They can be seen clearly in the photo below and the biology textbook diagram.

I also observe that that the cellular structure in both small and large diameter branches are of the same size(radius) and wall thickness, only fewer in number the smaller branches. This can be seen in the following picture: On the left is the cross section of the larger diameter branch while on the right is a smaller diameter branch further up the stem.

On the left is a cross section from of a large diameter stem while on the right is a smaller diameter stem.

As seen in all the other stem cross sections that I analysed, we can see that the density of cells increases as we head towards the edge of the stem cross section. This size gradient can be seen in the following photo. It is known that cells form on the outer edge. They start as smaller cells and grow in size, forming the centre cells.

Interesting observation are the radial direction fibers, especially seen in the green outer ring of the image below. I wonder what purpose they serve. I suspect it has something to do with preventing buckling, another failure mode the plant would not want to experience.

There are 4, 5 or 6 sided polygons that make up the shape of the cells. The sizes vary. In each region there is a majority size interspersed with smaller polygon to make everything fit. As far as I can tell, the cell walls are all of the same thickness. They just fit together so nicely. Maybe the original cells were circles but where squashed together to form the polygons. Maybe the polygons minimises some sort of parameter such as surface area. Maybe it uses the least surface area of cell walls while maximising the diameter of the tubes?

These cellular structures are of interest to me because I can potentially generate these structures as infill for 3D printing using an algorithm.

I started by looking at hexogonal structure deformation to generate varying sizes of cells, as seen in the plant cells. I came across the CAD program called Rhino and its parametric design plugin, Rhino. It allow the deformation of parametrically generated hex structures with what is known as “attractor points”. The effect of these attractor points can be seen in the image below. The attractor points are at the centre of distortion. The images are taken from here.

I also saw another image demonstrating the use of attractor points too:

centres of attraction highlighted in red crosses

I thought I can generate those distortion patterns myself and see if it is feasible to optimise the structure, like in a plant stem.

I first started by generating a hex field in Python’s Matplotlib as shown in the image below:

Radial Distortion

Then I used the radial distortion function used for correcting lenses to distort the hex field. The function is a simplied version of the one described on Wikipedia. The distortion function is described here:

x_\mathrm{d} = x_\mathrm{u} + (x_\mathrm{u} - x_\mathrm{c})(K_1r^2 + K_2r^4 )

y_\mathrm{d} = y_\mathrm{u} + (y_\mathrm{u} - y_\mathrm{c})(K_1r^2 + K_2r^4 )

(x_\mathrm{d},\ y_\mathrm{d}) = distorted image point,
(x_\mathrm{u},\ y_\mathrm{u}) = undistorted image point,
(x_\mathrm{c},\ y_\mathrm{c}) = distortion center,
K_n = n^{\mathrm{th}} radial distortion coefficient,
r = \sqrt{(x_\mathrm{d}-x_\mathrm{c})^2 + (y_\mathrm{d}-y_\mathrm{c})^2}

I take each point in the hex mesh and apply the formula for radial distortion to it and generate its new postion. The output is shown here, where K_1 = 0.0001 and K_2 = 0 .

First renders of radial distortion by me

The distortion is similar to radial distortion in lenses. If you know the k values of the lens, you will be able to undistort the image. But (automatically) determining the k values to remove radial distortion is a whole different research topic.

Stock photo showing severe radial distortion.

I think the closest I got to generating the structure of the stem was as shown below: On the left, where k_1 = -3e-5 we see a somewhat similar distribution of cell sizes as seen in the plant stem. The cells are largest in the center and gets smaller towards the edges.

I use a global distortion function so that the distortion is continuous, and the cells do not separate. However, you will see that the cells become malformed, where the angle between the sides of the polygon is no longer 120 degrees. This from a structural point of view may actually weaken the structure because I suspect only the perfect hexagon is strongest. I wonder if the weakening by deforming the hexagon cells actually defeats the purpose of using variable sized cells. This is an area to investigate.

Now I try the same thing, with 2 off-center distortion points:

2 off-center distortions

Electric Distortion

Impact of electric potential as the Q value varies. Q is substituted by a in my equation.
By Lookang – Own work, CC BY-SA 3.0,

I chose to try using another distortion formula, the one I think is used in the Rhino plugin. The function is used to describe the Electric potential due to a point charge as you move further an further away from the point charge.

x_\mathrm{d} = x_\mathrm{u} + (x_\mathrm{u} - x_\mathrm{c})( \frac{a}{r} )

y_\mathrm{d} = y_\mathrm{u} + (y_\mathrm{u} - y_\mathrm{c})( \frac{a}{r} )

(x_\mathrm{d},\ y_\mathrm{d}) = distorted image point,
(x_\mathrm{u},\ y_\mathrm{u}) = undistorted image point,
(x_\mathrm{c},\ y_\mathrm{c}) = distortion center,
a =  electric distortion coefficient,
r = \sqrt{(x_\mathrm{d}-x_\mathrm{c})^2 + (y_\mathrm{d}-y_\mathrm{c})^2}

I tried with multiple distortion centers, as shown below. We can see that the influence of the distortion points falls away very rapidly as we move away from the center. This formula I think has quite a localised impact, unlike the radial distortion function, which has an impact all the way to the edge of the hex field.

2 distortion centers(blue triangles) using electric distort. a = 30

I also then tried with 6 centres of distortion, arranged in a grid pattern. I used a much smaller distortion parameter a = 3(I did not log down the exact parameter value for this plot)

grid of distortion points, in a 3 by 3 pattern, a = 3

I then went on to randomly vary the a parameter values of each distortion point. a is a random value in {3,4,5}.

We need to limit the distortion to within the design space, not going out of bounds of the original mesh. To do that, I use a penalisation multiplier \frac{r^{'}}{r} ) to prevent distorted points from leaving the design space. Essentially, it neutralises the effect of the distortion function as it nears the edge of the design space.

Diagram showing the distortion penalisation variables.

x_\mathrm{d} = x_\mathrm{u} + (x_\mathrm{u} - x_\mathrm{c})( \frac{a}{r} ) (\frac{r^{'}}{r} )

y_\mathrm{d} = y_\mathrm{u} + (y_\mathrm{u} - y_\mathrm{c})( \frac{a}{r} ) (\frac{r^{'}}{r} )

r = distance from center of distortion (x_c,y_c) to edge of design space, collinear with point to be distorted (x,y) ,
r^{'} = distance from point to be distorted (x,y) to edge of design space, collinear with center of distortion (x_c,y_c) .


You may guess where I am heading toward now. Optimisation. We can vary the distortion parameters and locations of the centres of distortion to maximise stiffness of the structure. The designs considered so far have been 2D hex fields. In order to measure their stiffness when extruded into a column like a plant stem, we can take a short cut by considering the second moment of area as a measure of stiffness. The Euler–Bernoulli beam theory for slender beams gives the maximum bending stress of a beam as:

\sigma_{b} = - { My \over I_c }

where M is the bending moment at the location of interest along the beam’s length, I_c is the second moment of inertia of the beam’s cross section and y is the distance from the beam’s neutral axis to the point of interest along the height of the cross section. In order to minimse compliance, we need to minimise stress. When designing the beam(artificial plant stem), we can modify I_c ,the second moment of inertia. In fact, by maximising I_c, we can minimise \sigma_{b} and therefore compliance.

Calculating second moment of Inertia, I_c

The second moment of area is calculated by(ref):

Diagram showing y, dA and x axis

I_x = \int{y^2 ~ dA}

where y are the coordinates of element dA with respect to the x axis. In order to work with a bitmap image which I generate with Matplotlib, I consider every pixel as an element dA and calculate I_x by:

\sum_{all~dark~pixels} y^2 ~ dA

Now comes the optimisation part

We want to maximise I_x by varying the distortion parameter a for each distortion centre. The location of these centres are predefined(for now).

max~ I_x(A)

A is a vector of all distortion parameters,
I_x is a the second moment of area, the function of all the distortion parameters

Below is a simple case using the Constrained Optimization BY Linear Approximation COBYLA optimisation algorithm from the Scipy.Optimise Library with a single distortion point(at the centre). I like this algorithm because it usually converges very fast, taking big steps fast.

See how it converges to become like a hollow tube. However, we see that it does not become a complete tube because as the lines get squished together and overlap, their second moment of area I_x actually decreases. The optimiser sees no further improvements to I_x and converges. The algorithm does not allow overlapping regions, perfect for generating cellular structure.

The analytical optimal structure is a hollow cylinderical tube and this program achieves that. It is interesting to see the slight oscillation the algorithm goes through at the end before converging.

Irregular design space cross section

Since there is an analytical optimal solution for a cylinderical design space, I have decided to explore a more irregular U shape:

U shape design space

With a single distortion point:

Single distortion centres (blue triangle)

When we have 3 distortion centres where each parameter is controlled by the optimiser:

3 distortion centers(blue triangle)

However, there are some problems that pop up. For example, there are lines that overlap when in some cases where there are more distortion centers. For example, in the following animation:

overlapping lines for a 3 centre of distortion

At termination point, as show in the red circle below, there is an overlapping region that would be problematic if we need to make it into 3D cells.

The final convergence point yields very distorted polygonal cells. Infact, they do not resemble hexogons anymore. In the future, it may be worth considering adding a penalisation factor to force the hexogons to maintain their angles between sides.

The algorithm also teminates only when there is little change in the objective function. Perhaps it would be worth terminating early once a weigh objective is reached, such as 40% of original mass. This avoid creating these terribly distorted hex structures.

Next steps

I idea is to have these distortion centres in a 3D hexogonal lattice structure where there are centres of distortions embeded in the 3D space. The distortion parameters would be controlled by optimiser. The resultant structure would form the infill for 3D printed components, increasing strength to weight ratio even further, possibly more than traditional SIMP based topology optimisation algorithms. The ideal software to do this would be Rhino+Grasshopper plugin to generate 3D structures. Below is a collection of infills commonly used, with a very regular structure. We can do better than this, concentrating infill where it is most needed and reducing where not needed.,w=4000,gravity=0.5x0.5/wp-content/uploads/2018/06/26144450/nine-common-styles-of-infill-kronr-pinshape-180618.jpg?w=1140&ssl=1
A sample of 3D printed infills. Source: all3dp

It will also be critical to consider merging cells to save weight. In regions where there is no need for cells(walls), the cells can be merged. I think one source of natural inspiration would be bubble merging:

Bubble merging.

Lei Zhu suggests looking at Active set optimisation algorithm, an algorithm that allows selective constraints to save on computational resources. He also suggests looking at Level Set algorithm to look up boundary merging, like bubble merging. It is useful for nucleation, and can be used to merge cells. This is really promising: I can reduce the number of cells altogether in empty areas, and keep them only in areas where they are needed.

I think it is worth looking into other distortion functions such as Non-uniform rational basis spline (NURBS) :

NURBS surface, modified by moving those points. The shape follows the position of those points. Source: Rhino3d

Real Life Application

Testing in real life conditions is a must for all new ideas. I am privileged to be involved in Project Ardvark in the Imperial College Drone Society, led by Morgan Nightingale to build a long range fixed wing drone for National parks monitoring in Kenya. These are early days and I am looking forward to possibly using these ideas for 3D printed components on the Drone. It will be interesting to see how easy it will be to 3D print these unique infills and how they perform under load. So far, I have not considered buckling in the simulations and I do not know how it will perform.


I would like to thank my Supervisor, Dr Nan Li for guiding me on this exploration and for here very detailed critque of my work during our numerous face to face meeting.

I would also like to thank Lei Zhu for getting involved in the discussion and suggesting several algorithms for use.

I thank God, the maker of all things for his inspiration.

1 Comment

Leave a Reply

Your email address will not be published. Required fields are marked *