Tailored elastic properties of beam-based lattice unit structures

In this paper a structural optimization framework is developed to design three-dimensional periodic lattice unit cells that meets specific mechanical requirements. The work is motivated by the high design freedom of additive manufacturing technologies, which enable complex multiscale lattice structures to be printed. An optimized lattice unit cell delivers desired orthotropic elastic material properties, providing a tailored metamaterial. The design variables are the coordinates of lattice skeleton nodes defined within the three-dimensional lattice cell space, and the connectivities between them resulting a strut-skeleton. Genetic algorithm (GA) is combined with posterior particle swarm optimization (PSO) algorithm to establish an integrated topology and shape optimization tool. For the calculation of the elastic properties of the individual lattice cells, an effective Timoshenko beam-based finite element calculation method was developed. The novelty of the work stems from its free topology optimization nature, excluding the strut diameters from the optimization variables. The method is demonstrated by four lattice cell optimization cases, where extreme orthotropic elastic properties were targeted and achieved. The tailored lattice cells represent a metamaterial, that can be used to build a structural component on the macroscopic scale, by stacking the cells periodically together, to fill the macroscopic 3D design space. This framework is a strong basis that can be extended to meet further nonlinear metamaterial requirements, such as energy absorption.


Introduction
The paper aims to design tailored metamaterials in the form of lattice unit cell structures, which are potential to be produced through 3D-printing.Additive Manufacturing enables complex lattice structures to be easily produced, which would not be feasible with other manufacturing techniques.This work focuses only on the numerical simulation and design optimization of periodic lattice cells.The presented approach is Vol:.( 1234567890) independent from the choice of the printed material, as well as from manufacturing technology.Therefore, a homogeneous and isotropic reference material is assumed, which is independent from any 3D printing process-induced effects on the material quality and properties.
A repeating unit cell of the whole structure, which is commonly referred as representative volume element (RVE), could be represented through its homogenized structural properties.In this work, only beambased lattice unit cell geometries are optimized.The homogenized properties of the microscopic unit cell can then be assigned to a metamaterial, and might be used for analysis on the macroscopic lattice assembly level.The assembly refers to the cellular component, which is built together by stacked repeating unit cells.
The importance of using periodic cellular materials is constantly increasing, especially in the aerospace industry (Pantelakis and Tserpes 2020).Pantelakis and Tserpes (2020) give an overview of some basic unit cells, including their main application area and manufacturing techniques.Ashby studied different type of lattice and foam unit structures (Ashby 2006), identifying bending and stretching dominant lattice topologies, considering the effect of relative density.Similarly, Gibson and Ashby (1997) cluster the mechanical capabilities of some typical lattice and foam structures.
The traditional construction and optimization of periodic lattice structures happen through the commonly-used standard cells.Some exemplary standard cells are shown in Fig. 1.
The common trend of lattice structure customization is rather concentrated to size optimization problems of standard cell topologies, with predefined design variables.However, due to the high design freedom of additive manufacturing techniques, almost any arbitrary cell topology is feasible to produce.The complex microscale unit cells with extreme anisotopic properties may be integrated into macroscopic structural 3D problems to outperform construction through standard cells.Such a design workflow is presented by Schwahofer et al. (2022).The novel approach aims to optimize topologies with constant diameter of the lattice struts.

Overview
Various homogenization techniques are available to calculate the effective structural properties of a beambased lattice representative volume element (RVE).These are analytical approaches (Gibson and Ashby 1997), and numerical approaches such as asymptotic homogenization by Arabnejad and Pasini (2013).This technique was used by Marschall et al. (2020) for multiscale analysis of 3D-printed lattice sandwich structures.Dong et al. summarized the available numerical and experimental lattice homogenization techniques (Dong et al. 2017).
In this work, finite element-based homogenization is used for the elastic characterization of the unit cells.FE methods may vary according to the applied discretization method of the unit cell, which can be meshed as common practice with solid tetrahedral or hexahedral elements, and controlled periodic boundary conditions (PBC).In this work, an elastic FE beam model was established to reduce the computational effort of a homogenization.The virtual characterization can be carried out in the six principal deformation directions of the cell.Such a homogenization tool was implemented in Abaqus Scripting environment by Omairey et al. (2019), called easyPBC.

Homogenization through elastic beam elements
The main motivation behind the beam element-based FE modelling is the extremely low computational  Vol.: (0123456789) cost, compared to a solid element-based model.A slender strut-based lattice unit cell can be discretized through a few degrees of freedom.The main difficulty is the modelling of the lattice strut-joints (junctions), where more material is concentrated.The homogenization of two dimensional lattice structures by elastic elements was first established by Tollenaere and Caillerie (1998).Reis et al. (2010) analyzed analytically some standard 2D unit cells with beam elastic elements.Standard 2D unit cells were analyzed through elastic beam elements by Reis et al. (2010), Caillerie et al. (2006), and Hutchinson and Fleck (2006).Two dimensional beam-based macroscopic lattice structures were analyzed by Vigliotti and Pasini (2012a).This was enhanced to a threedimensional domain and homogenization of micro unit cell (Vigliotti and Pasini 2012b) through simple Euler-Bernoulli beams neglecting the additional stiffness of the junction points.Vigliotti et al. (2014) also researched nonlinear constitutive models for a 2D pin-jointed beam-based lattice.In the study of Park and Rosen (2018), as-fabricated modeling was applied considering the FFF printing technology and semi-rigid joint frame.
This paper targets to achieve a beam-based formulation to approximate reliably the elastic properties of an arbitrary 3D strut-based lattice unit.

Optimization techniques
An overview about cellular material design and optimization was written by several authors (Pan et al. 2020;Maconachie et al. 2019;Plocher and Panesar 2019).Truss-based ground structure optimization was firstly established by Zhang et al. (2016), where the macroelement approach was introduced as population-based strut layout optimization schemes.
In this paper, population-based genetic (Singh and Choudhary 2021;Kramer 2017;Tao et al. 2020) and particle swarm algorithm were selected for the fitness-based variation and optimization of unit cell topologies.For strut-based application such methods were used in several studies.Vaissier et al. (2019) optimized the support structure considering the manufacturability and targeting the least material used.In Su et al. (2009) Su et al. present a sparse matrix encoding scheme for ground structure based strut structure topology and sizing optimization with specially developed crossover and mutation operators.
Giger and Ermanni proposed an evolutionary algorithm for the optimization of strut structures (Giger and Ermanni 2006).
In recent years, a strut-based evolutionary level set topology optimization method was explored in large-scale crash design applications by Bujny et al. (2018a;b) targeting structural crashworthiness.Raponi et al. conducted research in a similar field using a Kriging-assisted level set method (Raponi et al. 2019a, b).Fairclough et al. (2021) established a 2D strut topology design tool that performs a ground structure optimization on a predefined 2D discretized domain.A projection-based ground structure optimization method was performed in the work of Deng and To (2020) for large-scale problems.
In the work of Chu et al. (2010), a lattice structure was optimized with Particle Swarm Optimization (PSO) and Least-squares Minimization (LSM) for derivation of the optimal size parameters of predefined lattice unit cell.Parameter optimization of strut members of a lattice component in the macroscale was carried out by Stanković et al. (2015) and Gorguluarslan et al. (2017).Sigmund (1994) first employed an inverse topology homogenization method to tailor 2D materials with prescribed constitutive parameters.A comprehensive review about multi-scale structural optimization was published by Wu et al. (2021).This review lists various further inverse homogenization methods considering elastic mechanical properties (Andreassen et al. 2014;Wang and Sigmund 2020), as well as other physics, such as thermal expansion (Sigmund and Torquato 1996), electrical conductivity (Torquato et al. 2002), or fluid permeability (Challis et al. 2012).Lattice unit cell size optimization on the microscale was carried out in recent papers (Imediegwu et al. 2019;Nightingale et al. 2021).These methods only tackled predefined beam-based unit cells.With such an analysis, tailored elastic properties are obtained.Other methods of tailoring elastic properties of lattice unit structures were explored by Chen et al. (2018) and Zhu et al. (2017) through high-fidelity optimization considering material property gamut.Xia and Breitkopf (2015a;b) carried out research in a similar field.In these studies, two-dimensional lattice units with extreme elastic properties were optimized.Da et al. (2017) published their hybrid cellular automaton method to design 2D micro-metamaterial with extreme properties.

Beam-based homogenization method
In the framework of this paper, only point symmetric geometries are considered, built of junction nodes and connectivity elements (struts).The unit cell might be cubic or a rectangular cuboid.Additionally, the constructed lattice skeletons always have a node at the eight corner points.The cross section of the skeleton members is circular.In the entire analysis framework, equal diameters of the beams within the lattice cell are considered.The modelling of the lattice cells with elastic beams happens through strut section refinement and reinforcement of junction elements.

Refinement with beam elements
The struts of the lattice cubes are modelled through multiple interconnected beam elements.The length and properties of the beam elements vary from each other.A beam is considered to be a junction element if it is positioned at a lattice junction node.The length of a junction element is determined relative to the smallest closed angle with another junction beam that is connected to the same junction point, and the diameter of the struts.
The varying and individually determined junction beam length is required to have a more realistic stiffness reinforcement, caused by the additional material concentrated at the connecting nodes of the cell.This will be revealed in detail in 3.2.
The method for the junction element length determination is illustrated in Fig. 2.Here a 2D view of a strut junction of three intersecting beams can be seen.With the knowledge of the smallest closed 3D angle and the r radius of the beams, the length of the junction is calculated through an easy trigonometric Eq. ( 1), where l 1 is the length of a junction beam, and 12 is the smallest closed angle.

Reinforcement of the junction elements
The refinement of the beam mesh of the cell is followed by the reinforcement of the junction beams.The junction points require slightly different treatment of inner joints and junctions on the boundary (1) l 1 = r sin( 12) .
surface of the RVE.In both cases, the collinearity of two adjacent beams is checked at every lattice node.
For junction points on the boundary surface, the periodic effect on the RVE boundaries must be considered.The junction clustering algorithm considers therefore the beam skeleton of the adjacent cell next to the RVE.Three exemplary beam crossing nodes are shown in Fig. 3, indicating the collinear condition for the identification as a junction node.The related beams are then selected for stiffness reinforcement.
The junction beam reinforcement can happen through Young modulus scaling as in Eq. (2) or through beam diameter scaling according to Eq. (3).In both cases, a scaling factor ( a mat or a sec ) is used to achieve the desired stiffening effect of the joint area.
(2) E reinf =E nom * a mat .In this work, both strategies were implemented and tested.For linear elastic structural analysis, the resulting homogenized properties are similar with both material-and section-based junction reinforcement.From now on, the material-based reinforcement will be used during the analysis.
Similarly to the junction beam length, the reinforcement factors of the junction beams may individually vary.The scaling factor is approximated based on the closed angle with respect to the principal deformation direction.The three-dimensional homogenization of the cells happens with respect to the six principal loading directions.These are the three normal loading directions (xx, yy and zz) and the three shear directions (xy, xz and yz).In every case, all the junction beams are assigned to an angle, that is closed with the normal plane of the principal deformation vector.This is called the orientation angle.
The body centered cube in Fig. 4 demonstrates the orientation angle under xx and zz directional normal loading.The orientation angle for the shear load cases are determined as the resultant vector of the shear components.As an example, the orientation direction corresponds to [1, 1, 0], in case of xy shearing.
The orientation angle dependent reinforcing model is based on the observation of how the differently oriented beams behave elastically during a linear strain increment.In horizontally and vertically oriented struts, compression or tension deformation seems to dominate.Furthermore, in struts with a tilt orientation, bending deformation may dominate.Due to the bending dominant loading, these struts may be ineffective in the joint area, as the junction beams are barely deformed due to the material concentration at the joint.The junction beams of the horizontal and vertical struts tend to contribute more to the elastic stiffness under pure tension or compression.This two step discretization through beam refinement and junction reinforcement is demonstrated in Fig. 5.
These observations were considered during the development of the reinforcement method of the junction beams.In order to find the realistic a mat reinforcement factor depending on the orientation angle of the junction beam, a numerical calibration was carried out.Five standard lattice cells were elastically homogenized with a high-fidelity (3) d reinf =d nom * a sec .tetrahedral volume mesh, as well as through the presented reduced-order beam-based finite element analysis.Furthermore, all the lattice geometries were homogenized with different diameter and cube length dimensions.The calibration of the beam model reinforcement parameter a mat was carried out in multiple homogenization directions.In this regard, the tetrahedral solid mesh-based homogenization is considered as a reference, and serves for the validation of the developed reduced order Timoshenko beam-based formulation.
Next, the strut orientation angle dependent model was established.This calibrated model is shown in Fig. 6.During the beam reinforcement process, the orientation angle is determined for all the junction beams, and their elastic modulus is assigned to the corresponding a mat reinforcement factor.
Fig. 4 The orientation angels of struts with respect to xx and zz directional compression Beyond the reinforcement factor, the validity range of the applied beam-based modelling was also analysed.During the calibration, the slenderness of five standard cells were varied, thus it was possible to determine the strut diameter -strut length ratio, where the Timoshenko beam-based model delivers accurate elastic properties.Accordingly, the shortest strut length of the lattice cell should be at least four times larger than the diameter in use.Based on the calibration of the reinforcement factors, any more slender unit structure always shows good approximating performance, even though very slender lattice would not need a junction reinforcement at all.For such thin structures, stiffness treatment of the junctions only minorly affects the calculated elastic properties.

Boundary conditions on beam-based mesh
During unit cell homogenization, periodic boundary conditions are desired.In the case of a solid-meshed FE model, there are multiple degrees of freedom through the cross section of the struts to clearly prescribe the dependencies of the opposing boundary edges (Omairey et al. 2019).
However, in case of the beam-based modelling, the corner and surface nodes are represented through one FE node, which makes it unclear how to prescribe the periodic boundary conditions.In the frame of this study, an extended boundary condition method for beam elements was implemented and tested.This approach prescribes periodic rotation boundary conditions (Eq. ( 4)) to the boundary nodes to consider the effect of the surrounding adjacent lattice cells.
where j = x, y, z .The +i, −i superscript indicates the opposing boundary surfaces of the unit. (4) Another important effect must be handled during the cross section assignment of the elastic beams that lie on the cell edge or surface.The struts on the cube edges have the form of a quarter circle cross section, while the beams lying on the cell surface are realized as a half circular strut.For this reason, the beam-based FE model has to inherit these cross section properties as well.With this treatment, it is possible to avoid unrealistic additional stiffness modelled on the boundary beams, which would lead to exaggerating the homogenized properties of the lattice cell.Through the implemented reduced-order method, the simulation time of a single directional homogenization can be reduced between 10 and 150 times, depending on the complexity of the unit cell topology.
The FE analysis is executed as elastic deformations in the six principal strain directions.The whole unit cell is meshed, and it is always points symmetric through the triple mirror principle.The detailed boundary conditions and equations are documented in Omairey et al. (2019).The additional constraints regarding the periodic rotations of the boundary junctions are stated in Eq. (4).Vol.: (0123456789)

Optimization method
In this paper, a population-based optimization technique was selected as it allows for potential extensions to design tailored metamaterial including energy absorption, as a response targeting crashworthiness applications.The two-step method firstly uses a genetic algorithm to generally optimize the unit cell topology.In this context, topology is used to only describe the way the separate nodes are connected, without considering the specific position of the nodes.
Nodes can belong to specific domains (volume nodes, surface nodes, corner nodes).The result of the GA ideally is an optimal topology, with suboptimal node position being permissible.Due to the high degree of freedom and the difference in paradigm of the continuous optimization for node positions and the discrete optimization for topology, having both aspects converge simultaneously is hard to achieve.The impact of the possible discrete changes in topology (removing/adding nodes/beams) is assumed to outweigh the impact of the continuous change of node positions in each time step.Consequently, the genetic fitnessbased selection operation will favor good topology change more than accurate node positions.However, simultaneous optimization of the node positions in the GA is necessary to achieve a minimum degree of node position optimization.This is required so that a good or optimal topology with suboptimal node positions still outperforms worse topologies that happen to be closer to their ideal node placement by chance.
If the cell with superior topology is not optimized far enough regarding its node positions the wrong cell would be selected for the next step.
Then, a reduced degree of freedom optimization is done posterior to the main GA loop in the form of a PSO algorithm which will find the optimal node position for the specific input topology.The posterior PSO loop can be treated as a shape optimization of the prior topology optimization in the GA.The approach is implemented in a python script responsible for genetic operations and optimization, evaluating each genotype through the beam-based Abaqus-Python interface (Sect.3).

Encoding of cells
Each genotype consists of two lists, with one describing the nodes and one defining the beams.In Fig. 7, the data structure for the topology junctions is demonstrated.Each node is given a unique index, its position in cartesian coordinates, an index representing its domain, the number of connections that are allocated for the node and the actual number of connections of the node.The domain describes what location change to the node is possible.A "0" indicates a "volumenode" with free movement, 1-6 are the surface nodes constrained to movement on the corresponding surfaces of the cuboid unit cell, and 10 are the corner nodes that have no free movement.The connections or beams are then defined in a list through pairs of node indices between which a beam is positioned.
The unit cells are optimized on 1/8th sub-cell, which is triple mirrored across the three center planes to obtained the targeted periodical and point symmetric lattice unit design (Fig. 8).

Connection generation
New connections are generated when the allocated number of connections for a node is higher than the actual number.To find a partner for such a node, a list of all other nodes without full connections  Vol:.( 1234567890) is generated and a potential partner is randomly selected.Then the user selectable minimum distance to other connections (for this paper: 2 × diameter ), and a minimum connection length (for this paper: 2 × diameter ) is checked to verify if the new beam is allowed.These conditions also help to generate structurally more meaningful topologies, furthermore to fit to the validation range of the elastic beam calculation framework.This process is repeated until no new connections can be found.After this the full internal connectivity of the cell and the minimum number of connections per node (2 connections for volume nodes, 1 connection for surface and corner nodes) are checked.When these conditions are not met, they trigger a repair function that removes all connections from the nodes and tries to generate a new set of connections from the ground up.The connection generation step also repairs non-connected graphs and beams that have positioned too close or too short through the mating and mutation step.

Mating
In mating different partners are selected to be recombined to form the next generation of cells.The GA optimization loop proposed has been tested with and without mating where both approaches have been successful.When not using mating, the next generation is created through a fitness proportional selection (Kramer 2017).To fill the population size with selected cells elitism can also be enabled in the GA, which allocates a fraction of the population size to a selection from the previous best cells of each generation (Kramer 2017).Mating is implemented as a nonexclusive two-parent scheme.The parents are picked by fitness proportional selection with a child cell having the same parent twice being a possibility.To include the idea of elitism, an option allows the inclusion of a fraction of the total population of previous best cells into the parent pool which in turn allows their selection as parents.
For each parent pair to form their child cell, features need to undergo crossover.The nodes from both parents are matched with each other through an assignment optimization where the total euclidian distance between point pairs is minimized.This assignment approach does not require an equal number of nodes in the parent cells and can consider the nodes in each domain (volume, surface) separately.To solve the assignment problem the weighted frequency method by Habr et al. is implemented (Gottschlich and Schuhmacher 2014;Schrage 2016).
After all the possible assignments have been made, a uniform random factor is used to choose if the excess nodes from the parent with more nodes are used or discarded.The crossover of the positions is done through the arithmetic mean of both matched parent nodes.For the connections as well as the allocated number of connections a uniform random choice picks one parent from which to inherit.

Mutation
Mutation consists of a step-by-step method where each descriptive trait of the unit cell can be mutated according to a user picked probability.Node positions are mutated with a vector uniformly distributed in direction and gaussian in length.The nodes that are affected are selected by a weighted random choice.The movement with the vector is checked to ensure the node-to-node minimum distance constraints are not broken, and to not leave the external dimensions of the cell.In the case of a breach a new vector is generated and tried.The number of connections allocated to each node is also mutated.Each node has a weighted random choice applied and if selected an integer rounded from a gaussian distribution is either added or subtracted from the allocated connection number.Here a constraint for minimum and maximum connections is applied.Each node has a random chance of being deleted.The connection list is cleaned after deletion where all connections to nonexistent nodes are removed.Then a random chance is applied to create a new node respecting the distance constraints.Last, a chance for complete re-connection of a cell is applied.

Optimization objective
The general implementation in the fitness functions is a weight-based transformation from a multi objective optimization to a single combined objective.This approach works well for the different implemented fitness function.The response of each unit cell is evaluated based on the homogenized properties of a continuous periodic lattice of the specific unit cell.

Cell weight approximation
To approximate the weight specific performance of the unit cells the infill percentage is used as a surrogate for density.To find the infill percentage all nodes are considered as spheres and all connections as cylinders.The volumes of entities on the cell surface and edge are counted accordingly.Additionally, the length of each beam is reduced by its diameter multiplied by an empirical factor compensating for overlap of the beams.

Cell infill control
There is no direct method to specify the infill of the cells throughout an optimization.However, the infill is restricted through the constraints and the connection number penalty factor ( c con ) in combination with the beam diameter and the unit cell dimensions.The complexity constraints, namely minimum connection length, minimum node-to-node distance, minimum connection-to-connection distance and minimum connection-to-node distance, affect how the beams can be packed into the space and can reduce the maximum possible infill ratio.
An important aspect regarding infill is the fact that new nodes and beams can be created or removed in the cells in each generation throughout the GA, thus changing the infill of the randomly generated initial population (which is possible) will have no guarantee on the final infill of the best performing cell.

Maximum connection and cell recursion
Without sizing aspect within the optimization, cell recursion could occur (Fig. 9).The only difference in fitness between the simple full cell and the 1/8th sub-cell of the same topology is the beam diameter relative to its size.To stop the algorithm from finding increasingly complex recursions of the same cell topology a maximum connection number n max is set, which when passed ( n is > n max ) induces a penalty ( c con ) for the fitness ( f = optimization value∕c con ) according to Eq. ( 5).This stops additional computational time caused by redundant complex cells, and also stops a pseudo sizing optimization from occurring through this recursion.
The specific choice of n max = 12 for this paper is based on the cube recursion example.Where the simple cube has 3 beams in the 1/8th cell, the first recursion has 12 leading to the choice that this number would generally allow for enough complexity for most optimization goals, while still penalizing the expected recursion problem for quite simple cells.Because this limit is implemented as a cost function as opposed to a hard constraint the other contributing factors in the fitness functions can also outweigh this factor to find the right cell, with only a bias towards a simpler solution.

Stiffness ratio
For this type of fitness function, the stiffness per infill percentage in the main optimization direction is maximized: To control the stiffnesses in the other principal axes relative to the main axis, a ratio is defined, and the cost function shown in Eq. ( 7) is applied for both ratios ( r 1 ,r 2 ).
Different weights ( w 1 ) can be applied to change the influence of these ratio costs on the weight specific main axis stiffness.Through this, selected cell (5) properties which are more important to the specific application can be encouraged.The final fitness is then calculated by dividing by the costs: The formulation of this fitness function with the combination of a main goal and cost functions for the stiffness ratios as well as the equation itself was largely found experimentally.Variations of this such as a reduced power ( 1 ≤ exponent ≤ 2 instead of 3) also led to successful optimizations but were more prone to stagnate in local minima and take more generations to generate well-fitting solutions if found at all.At higher powers ( 4 ≤ exponent ) the selec- tion of cells with high fitness is dominated by only very few candidates due to the stronger scaling with ratio distance from the target ratio, leading to poor exploration.
The hard geometric constraints such as the minimum node-to-node distance, the minimum beamto-beam distance, the minimum beam-to-node distance, the minimum connection length and the minimum beam-to-beam angle, are all enforced outside of the fitness function within the cell altering operations themselves and will thus never be broken. (8)

Poisson ratio
The poisson ratio is implemented as a cost function to approach a simple target value.The main optimization variable is still weight specific stiffness, so that an optimum exists.The connection factor cost is also used.Each of the 6 poisson ratios from the orthotropic material definition 12 , 13 , 21 , 23 , 31 , 32 can be set to a target and weighted independently by a factor w ij .Equation ( 9) describes the cost c ij that is calculated for each poisson ratio.
This cost function can also be used for negative poisson values.The final cell fitness is calculated in the same way as for the stiffness ratio formulation in Eq. ( 8).

Population size
The minimum viable population size is very problem dependant and cannot be generalized from the amount of data collected.In Fig. 10 different exemplary optimization problems with different population sizes can be seen.While the smallest population size shows signs of no continuous evolution, larger populations seem to work well.The population of 60 was also not sufficient for the given problem, but still shows continuous improvement in the cumulative best.If real Fig. 10 Effects of population size on convergence of the genetic algorithm convergence is desired increasing the population size seems effective in eliminating noise and allowing for node position evolution concurrently with topological evolution.But due to the long runtime large population sizes become less feasible with the current implementation even if they show nicer behaving evolution.

PSO loop
The discrete topology optimization through GA has a high frequency of topology mutation.That makes it difficult to achieve smoothly and continuously optimized structure modifications.To counteract a need for the reduction in mutation rate to a point where the GA can find the optimal node position, a fast-converging population-based approach (Hembecker et al. 2007) is applied to the reduced degree of freedom problem of optimizing node positions for an unchanging topology.Here the particle swarm optimization is chosen.PSO is a robust algorithm that is highly regarded and is applicable to a variety of diverse problems (Okwu and Tartibu 2021).It has also been generalized that this approach is robust and efficient even for difficult problems (Hembecker et al. 2007).The PSO algorithm is initialized by having a population of randomly distributed particles, each at a point ( x i (t = 0) ) in the solution space, with an initial random velocity vector ( v i (t = 0) ).The movement of each particle is defined by its velocity for each time step seen in Eq. ( 10) (Okwu and Tartibu 2021).The swarm behavior is introduced through the calculation of the velocity which includes the positions of the global best ( x global Best ) and local best ( x local Best ).The local best denotes the best position that each specific point has reached across all iterations, with the global best denoting the absolute best position found.
A random factor between 0 and 1 ( r 1 , r 2 ) and a weight ( c 1 , c 2 ) are used to induce exploration of the solution space (Hembecker et al. 2007).An inertial mass factor (w) is also applied to the previous velocity v i (t) to set an efficient trade-off between global exploration and local refinement of the solution (Shi and Eberhart 1998).Equation (11) shows the method to calculate the velocity for a point for the next iteration (Shi and Eberhart 1998). (10) To apply the PSO principles to the lattice unit cell, each node is considered as a separate optimization problem with all of them running concurrently with their local best taken from the best performing configuration of a unit cell across its own history, and the global best position taken from the single best solution, for all nodes.The initial population is generated by importing the cell with the highest fitness from the GA loop.Then, a random vector of Gauss distributed magnitude and uniform distributed direction is added to each node.This process is repeated to form the initial population of cells.

Stiffness optimization
In this chapter, the capabilities and performance of the developed reduced order elastic beam-based calculation tool and structural optimization framework are demonstrated through four unit cell topology problems, where diverse extreme anisotropic stiffness properties were targeted.

Stiffness optimization with large strut diameter
The target of the optimization is a weight specific maximum y-axis stiffness ( E yy ) with cost functions to con- trol the other ratios of x/y stiffness ( E xx ∕E yy ) and z/y stiffness ( E zz ∕E yy ).The x directional stiffness was tar- geted to be 10 percent of the maximized y directional stiffness, while one percent in z direction.The design space of the 1/8th cell can be seen in Eq. ( 12), while the main objective is stated in Eq. ( 13).The two ratios for the cost functions (Eq.( 7)) are indicated in Eq. ( 14), including their weights in Eq. ( 15).The complete fitness is determined according to Eq. ( 8). (11) For this optimization the evolution rate seen in Fig. 12 is as expected, with the increased slope from generations 3 to 10 representing the biggest change in the topology by removing insignificant beams.The topology optimization however functions well even after this point and new topologies are found that are better able to match the optimization criteria, which can be seen in the sample cells (c) and (d) in Fig. 11.
Here, the cell complexity is increased from the low complexity starting point of generation 10 by incorporating two beams in the place of 1 to adjust stiffness by changing the loading direction within the cell, so that it becomes possible to fulfill the optimization goal.In this optimization a population size of 200 was used, however in evolution chain of the best performing cells only a single approach is represented that is being altered, while all other niches are short lived without significant contribution towards (15) w 1 = 0.2 and w 2 = 0.2 the solution.This may suggest insufficient diversity within the candidate solutions, possibly due to a lack of niches forming due to no additional method used to encourage this, as well as aggressive elitism being used (Fig. 12).
The particle swarm optimization (Fig. 13) for the cell topology from Fig. 11 generation 38, is very convergent both for the best solution as well as for the average of all.Thus, it can be inferred that no solutions are trapped in a local optimum and all solutions are converging toward the same node positions.
The optimization may have reached the actual optimum with the weighting that was used, but as can be seen in Table 1 did not fully match the prescribed goal ratio for ratio 2. To show that this is the global optimum multiple reruns of the same algorithm would be required which was not yet tested due to the comparatively high run time of this process.
In an attempt to increase this accuracy, parameter constraints such as a minimum beam length of two times the beam diameter are imposed on the movement of the nodes throughout the optimization, which limit the solution space, but are a reasonable trade-off for higher accuracy for use with this fast simulation method.
Nevertheless, the large strut radius of 0.5 compared to the lattice unit cell size of a 10 × 10 × 10 cube, leads to the relatively high inaccuracy of the beam model results compared to a solid-meshed FEM analysis (Table 1).The homogenized elastic properties through the tetrahedral are considered as reference values.In the next demonstrated optimization problems, in 5.1.2,5.2.1 and 5.2.2, a strut radius of 0.1 is used which bring closer agreement between the beam and solid-based calculation of the optimized lattice skeletons.With reduced strut diameter, it is more likely that the optimizer constructs longer struts, that fit better to the validity range of the elastic beam model.
The cell is optimized to fit the unique elastic behaviour defined through the fitness function.Thus, comparing the cell performance with standard cells on basis of the fitness function will not be meaningful.The 'octet' unit cell (Table 5, Cell (f)) yields a fitness of 0.9907 for this function compared to the 24,120 of the actual result.The strong response of the fitness function to the allocated ratios is seen here, which is part of the reason why cells are found which fit these ratios so closely.This behaviour can be moderated through the weight factor when determining the ratio costs ( c ratio ), thus a more lenient ratio optimi- zation goal that leads to cells with the focus more on weight specific stiffness, instead of the specific ratios, can be defined as well.

Stiffness optimization with reduced strut diameter
For this optimization problem the weight specific maximum z-axis stiffness was targeted, with cost functions of ratio of x/z stiffness (7.9 %) and ratio of y/z stiffness (8.8 %).The inputs for the fitness function are stated in Eqs. ( 16), ( 17), ( 18), ( 19).
The proposed algorithm was able to find a topology which matched the prescribed stiffness rations better (Table 2) than for the first stiffness optimization.However no claim can be made that either have or have not found the absolute optimum, only that a useful unit cell is generated, that suits the application for which the fitness parameters were set.In the fitness evolution with the GA (Fig. 15) different topologies are evolved simultaneously with versions of the topology from cell (c) in Fig. 14 and cell (d) co-evolving and competing for the highest fitness within each generation.This is a desired behaviour and is working well in this optimization with a population size of 100 (Fig. 15).
( The topology which ends up with the highest fitness is first tried in generation 36 of the GA and is improved until generation 62. Due to the high degree of freedom, when allowing topology change and node position change concurrently, this process only slowly converges on the local or global optimum.When comparing this with the convergence of the PSO (Fig. 16) where within 6 generations (population size 80) the node position converges, the advantage of using this joint method is clear.The high randomness of mating with different topologies as well as mutation of the topology in the GA will effectively only have a few to no cells of the "best" topology being solely optimized for node positions within each generation, which is inefficient, while the PSO can do the positional optimization in a more focused manner.

Positive poisson ratio optimization
For the poisson ratio optimization there are potentially more cost functions, if each of the 6 poisson ratios were being optimized, which would lead to a higher problem complexity, but in this use case only 3 were chosen by weighting the others with 0. The unit cell here is a cuboid with the dimensions 6 × 10 × 20 .As a consequence the poisson ratio goal 31 = 2 and 32 = 2 need to be achieved.
The space for the design variables as well as inputs for the cost functions of poisson ratio-based formulation from Eq. ( 9) are shown below (Eqs.( 20), ( 21), ( 22)):   Vol.: (0123456789) In the fitness evolution of the GA for this optimization (Fig. 18) a high randomness and large sections of fitness degradation can be seen.This is a problem inherent property which is especially for the poisson ratio optimizations as they are much more sensitive to a change in topology.The discrete nature of adding or removing beams, even at low mutation rates where per generation only a single topology element is changed, can lead to a large change in fitness.If this type of fitness decrease due to exploration happens for all the well performing cells within a generation a significant fitness decrease is seen.In the algorithm elitism is used, meaning that the previous single best cell is always included in the next generation, and a random selection of previous well performing cells are added to the mating pool.But due to the sensitivity of the cell performance to topology change generations of low performance, such as from generation 39 to 43 (Fig. 18), this is expected and not harmful as the elitism allows for the previously well performing population to be reconsidered.For this complex goal of a poisson ration of 2 the resulting cell topologies become more complex (Fig. 17).Though the same kind of gradual topology change through exploration of the solution space as was seen with the simple examples can be seen here, when from cell c (generation 38) to cell d (generation 59) (Fig. 17), two beams are removed and one is added in an advantageous way over the course of 21 generations.The convergence of the genetic algorithm is plotted in Fig. 18.The quick convergence of the particle swarm optimization (Fig. 19) can be observed as usual.
In this optimization the effects of the weights of the objectives are evident.The matching of 31 and 32 to the goal is almost perfectly achieved, while the topology that allows this is not able to be matched to the goal for 21 , which is regarded less by the optimization due to the selected weights.The accuracy of the beam-based homogenization on the optimized topology shows acceptable accuracy with a maximal error of 13 %.The results of (20)  3.

Negative poisson ratio optimization
In this problem, two poisson ratios were targeted to reach negative values.The cost function members are defined according to Eq. ( 24) to formulate Eq.The weights the fitness as well as the design space are declared in Eqs. ( 23) and ( 25).
(23) x i , y i , z i ∈ [0, 5.0] (24) 31 = −1 and 32 = −1 The convergence of the GA seems to be poor, as can be seen from Fig. 21, where after the single peak at generation 68 no further improvement occurs.Here a saturation of the generations occurs where a strong local optimum dominates the candidates.The solution from the peak however has a high enough potential that after the optimization with the PSO the goal is perfectly met.
From the sample cells in Fig. 20, the variety of different topologies that are attempted, is seen.When a problem such as these negative poisson ratios require (25) w 12 , w 13 , w 23 , w 21 = 0 and w 31 , w 32 = 10  very specific unit cells to even exhibit the wanted behavior at all, the exploration of the algorithm is shown to be sufficient to find a candidate that can do this (Fig. 21).
The way the fitness function is implemented responds strongly to small changes in the goal parameter, so that the results may converge with high accuracy.This leads to the disparity of the GA reaching a relatively low maximum fitness, while the PSO (Fig. 22) converges at around a forty times higher fitness value than the initial GA.
The accuracy of the beam-based approximation is very good for the three normal stiffness properties ( E 11 , E 22 , E 33 ) according to Table 4, but seems to be less accurate for the negative poisson ratios.Due to the complex nature of the deformation under negative poisson ratio, the reduced order beam-meshed model cannot predict the actual poisson number so accurately.Nevertheless, the tetrahedral FE calculations confirmed the negative poisson ratios, which are considered as very good results, considering the very challenging nature of the demonstrated optimization problem.Due to the extremely complex design, the manufacturability of this lattice (Fig. 20) would be probably not feasible.
All in all, the algorithm performs well for different optimization variables between the stiffness optimization and the poisson ratio optimization, where the latter requires much more specific topologies compared with stiffness.O. Schwahofer cell properties.Main objective and cost ratios r 1 and r 2 are set according to equa- tions below (( 26), ( 27), ( 28), ( 29)).The beam radius in this optimization is set to 0.75.For the cubic optimization n max = 10 maximum connection number is used in the fitness function from Eq. ( 8).
(26) x i , y i , z i ∈ [0, 5.0] In Table 5 the simulation results for standard unit cells as well as their respective fitness for this function can be seen.The cube cell (Table 5, Cell (d)) outperforms all others in weight specific performance, with the intuitive explanation that it has ideal material alignment for loading along the major axes when disregarding shear stiffness.
Based on these samples from the solution space an indication of convergence is for the GA to find the cube cell or a better performing cell.For the two attempts at this fitness function, cubic optimization 1 (Figs. 23,25) and cubic optimization 2 (Figs.24,26), the GA successfully reached the cube cell.Consequently, some confidence of repeatability and The initial population for optimization 2 uses more complex (higher number of nodes and number of beams) cells to start the optimization, as can be seen in cells (a) and (b) in Fig. 24, compared to the relatively simple cells from optimization 1 in Fig. 23.This is reflected in the number of generations, with optimization 1 (Fig. 25) requiring 11 generations in contrast to the 30 generations of optimization 2 (Fig. 26).Here the extra generations are required for enough discrete changes (beam and node deletions and creations) to derive the cube from the complex starting point.
While there is some indication that in a 'sufficiently large, long running' optimization the algorithm could be perfectly convergent, this argument is essentially not important for the proposed method.The definition of population size and generation count to be sufficient is unclear for each problem and is likely to become prohibitively large, thus the general goal in the sample optimizations throughout this paper is to find a good solution for the posed fitness function with no consideration if it is optimal or not.
In contrast to a ground structure optimization nodes and beams can be both created and removed, which should mean that with enough generations any solution should be able to be transformed into any other solution.This further supports the idea that the starting population is not relevant in determining whether the GA finds the optimum only how long it will take as the method of mutation and mating limits the changes possible per generation.When looking at the changes to the cells in cubic optimization 1 (Fig. 23), for the best cells in each generation only small changes like adding and removing 1 or 2 beams occur.This is not indicative of the exploration of the algorithm.In Fig. 27 6 randomly selected cells from generation 11 in cubic optimization 1, where the cube cell was found, are shown.The intra generation variability displayed here is very high showing good exploration of the solution space.

General remarks
Based on other observed cell optimizations, convergence to the same result takes place, when the objective was moderately complex and the same geometric complexity constraints were used.With differing geometric complexity (not initial complexity but constraints such as minimum connection length and beam diameter allowing for more nodes and beam to fit into the cell) the optimal design will differ if there is a topology that fulfills the fitness function more closely.However, simple solutions such as the cube cell from Sect.5.3.1 will always be contained in the solution space even when more complex topologies are also considered.The fitness function is formulated in such a manner that it is insensitive to complexity (apart form the maximum connection number) as an expected increased stiffness due to more beams is compensated by the increased infill ratio.
However, for very complex objectives (such as negative passion ratio from the paper), the final design is probably not exactly optimal and a repeated run could converge to different but similarly performing unit cell solutions, even if the same geometric complexity is applied.
The question of convergence is much simpler for the PSO results due to the good convergence for the best cell as well as a general tendency for all other cells in the population (average line) to move to the optimum (unclear if local or global) being evident (Figs. 13,16,19,22).The PSO is irrelevant to this cube cell test as all nodes in the best cell are corner nodes that cannot be moved.

Periodic continuation of a cell
Figure 28 shows the periodic continuation of the converged cell from the first stiffness optimization problem from 5.1.1.Due to the constraints within the algorithm to have minimum distances to each boundary surface, this continuation can be done for all the cells and give reasonable results, that are usable in a real application.Therefore, a component can be built by filling its design space with the optimized cells, representing the orthotropic elastic metamaterial, that was designed by the optimizer.This lattice assembly was also 3D-printed with a commercial SLA printer using resin (Fig. 29).

Summary
In this work, a beam-based calculation framework was presented that approximates the effective properties of microscopic strut-like lattice unit structures.Within a certain validity range of strut diameter and strut length ratio, the homogenized properties through the beam-based simulation shows good accuracy.In the finite element models of the lattice cells, a perfectly isotropic material was assumed, neglecting process-induced effects on the material properties.The highly reduced order beam-based calculations are used to calculate the design responses in a combined GA-PSO structural optimization loop.The optimization tool can design a lightweight lattice unit cell, providing a desired orthotropic elastic and lightweight metamaterial.
The optimization goals that were set for the examples are quite extreme and the proposed approach can still generate a topology to mostly match them.
Through the introduced optimization method, three-dimensional and 3D-printable orthotropic lattice unit cells can be tailored.Similar research for unit cells was performed only in Chen et al. (2018), Zhu et al. (2017) and Xia and Breitkopf (2015a).As it is summarized in Sect.2.2, meta-heuristic optimization methods were often used in literature mostly for nonlinear structural cases, but was not yet explored for such microscopic unit cell problem.The promising linear elastic cell tailoring results are encouraging to add nonlinear responses to fully benefit from the method to find ideal elastoplastic properties of the tailored metamaterial.

Outlook
The algorithm performs well for different optimization variables between the stiffness optimization and the poisson ratio optimization.Thus it is extendable to other quantifiable structural properties such as unidirectional strength and energy absorption, enabling crash worthiness as an objective on the macroscopic component to be considered.The main motivation behind using a genetic algorithm is the easy expansion to nonlinear structural responses under plastic deformation.The integration of other non-structural responses might also be possible, such as wave propagation or thermal properties.Such an approach would enable a compact multidisciplinary formulation of the metamaterial optimization through strut-like lattice cells.
An application-specific unit cell from this approach could outperform conventional unit cells, when it is integrated into a macroscopic structural problem.The enhancement of this framework for the integration of 3D macroscale problems is revealed in Schwahofer et al. (2022).
The ability to explicitly control and optimize the relative density of a tailored cell is not available at the moment.An adjusted fitness with infill objective could be formulated, but the pseudo-discrete material allocation through constant strut diameter could lead to convergence difficulties.
Further development of the optimization method by neighborhood mating in the GA could increase niching in the candidate solutions providing more complex topologies, while a PSO integrated GA could also bring convergence benefits.Some of the optimized cells in Sect. 5 have converged to an extremely complex topology.Manufacturing such objects would lead to printability issues.However, the manufacturability constraints of the optimization algorithm were in this paper deactivated, as the exploration of capabilities of the numerical method was in focus.More advanced lattice cells could be derived by removing the triple-symmetry constraint and introducing a periodic geometrical hard constraint on the 3D unit domain, allowing free inner topology.Projekt DEAL.

Data availibility
The datasets generated during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.

Fig. 1
Fig. 1 Some standard lattice unit cells

Fig. 2 Fig. 3
Fig. 2 The determination of junction beam length and closed angles

Fig. 7
Fig. 7 Data structure of nodes

Fig. 13
Fig. 11 Cell samples from evolution for stiffness ratios (large beam diameter)

Fig. 15 Fig. 16
Fig. 15Fitness evolution of genetic algorithm for stiffness ratios (small beam diameter)

Fig. 14
Fig.14Cell samples from evolution for stiffness ratios (small beam diameter)

Fig. 19
Fig. 17 Cell samples from evolution for poisson ratios

Fig. 20 Fig. 22
Fig. 20 Cell samples from evolution for negative poisson ratios

Fig. 26
Fig. 23 Cell samples from evolution for cubic goal (simple random starting cells)

Fig. 29
Fig. 29 3D-printed demonstrator of the tailored cell from 5.1.1 Index node x pos y pos z pos Index surf num con goal num con is

Table 2
Performance of unit cell stiffness ratio

Table 4
Performance of unit cell: negative poisson ratio

Table 5
Fitness and properties of standard unit cell for cubic stiffness ratios ( r 1 = 1 , r 2 = 1)