Optimal design of functionally graded lattice structures using Hencky bar-grid model and topology optimization

Presented herein is a novel design framework for obtaining the optimal design of functionally graded lattice (FGL) structures that involve using a physical discrete structural model called the Hencky bar-grid model (HBM) and topology optimization (TO). The continuous FGL structure is discretized by HBM comprising rigid bars, frictionless hinges, frictionless pulleys, elastic primary and secondary axial springs, and torsional springs. A penalty function is introduced to each of the HBM spring’s stiffnesses to model non-uniform material properties. The gradient-based TO method is applied to find the stiffest structure via minimizing the compliance or elastic strain energy by adjusting the HBM spring stiffnesses subjected to prescribed design constraints. The optimal design of FGL structures is constructed based on the optimal spring stiffnesses of the HBM. The proposed design framework is simple to implement and for obtaining optimal FGL structures as it involves a relatively small number of design variables such as the spring stiffnesses of each grid cell. As illustration of the HBM-TO method, some optimization problems of FGL structures are considered and their optimal solutions obtained. The solutions are shown to converge after a small number of iterations. A Python code is given in the Appendix for interested readers who wish to reproduce the results.


Introduction
Additive manufacturing opens up new opportunities to fabricate as-designed functionally graded lattice (FGL) structures by building up material layer-by-layer (Maconachie et al. 2019). Lattice structures (also referred to as cellular structures) have a clear network feature comprising basic elements like solid struts or plates (Gibson 1989). Properly designed FGL structures can achieve high performance with great control in stiffness and strength, weight, energy absorption, heat exchange and acoustic insulation (Zhang et al. 2015;Ferro et al. 2022). Generally, elastic structures with minimum compliance or minimum elastic strain energy are the stiffest (Hassani and Hinton 1999;Kaveh et al. 2008;Zhang et al. 2014). Topology optimization (TO) method is commonly adopted to minimize the compliance or elastic strain energy of FGL structures by iteratively updating the material distribution subjected to loading, boundary conditions and prescribed deflection or stress constraints (Panesar et al. 2018;Liu et al. 2018;Zhu et al. 2021).
Many design strategies have been proposed for obtaining optimal designs of FGL structures with prescribed mechanical properties. Most of the existing strategies may 1 3 276 Page 2 of 23 be grouped into two categories. One category is the greyscale density mapping (GDM) strategy. This design strategy first uses standard TO algorithms such as the solid isotropic material with penalization (SIMP) method (Bendsøe 1989;Zhou and Rozvany 1991;Mlejnek 1992;Sigmund 2001;Andreassen et al. 2011;Wang et al. 2022), the evolutionary structural optimization method (Querin et al. 1998(Querin et al. , 2000Huang and Xie 2007), or the level set method (Wang et al. 2003;Challis and Guest 2009;Challis 2010;van Dijk et al. 2013;Azari Nejat et al. 2022;Lin et al. 2022) with finite element discretization or Isogeometric Analysis (Guerder et al. 2022) to determine a greyscale density solution of material distribution. Next, FGL structures with a selected type-based unit cell are established according to the obtained material distribution. The other category is the ground structure optimization (GSO) strategy. In this strategy, physical structures such as trusses, beams or frames are used to represent the FGL structures. Gradient-based and derivative free methods have been used to optimize the cross-section dimensions, node positions and the structural topology (or element connectivity) of truss, beam and frame subjected to prescribed design requirements such as volume ratio, deflection and stress (Miguel et al. 2013;Nguyen et al. 2013;Zhang et al. 2015;Han and Lu 2018). By comparing these two main strategies, GDM is more efficient because it needs to optimize a smaller number of design variables associated with each structural element. In contrast, GSO requires the optimization of a much larger number of design variables including material properties, node positions and the structure topology (or element connectivity) (Miguel et al. 2013). However, the GSO strategy provides a clearer physical structural configuration without the need of designing the internal structure and the connectivity for each lattice unit cell (Panesar et al. 2018).
In this paper, we propose a new strategy for the optimal design of FGL structures. This new strategy involves using a recently proposed lattice structure model called the Hencky bar-grid model (HBM) (Zhang et al. 2021a(Zhang et al. , b, 2022. A Hencky bar model is a type of bar-spring structural model that was pioneered by Hencky (1921). Owing to its simple and clear physical representation and great flexibility, it had been recently revisited and expanded to study various static, buckling and dynamic problems of different structural types (Wang et al. 2017(Wang et al. , 2020Zhang et al. 2018cZhang et al. , b, a, d, 2019aZhang et al. , b, 2022. HBM comprises bars and springs whose stiffnesses may be adjusted to allow for different material property distributions. The optimization of the spring stiffnesses is performed by using a gradient-based TO method to minimize the compliance while satisfying required constraints. The proposed design framework combines the advantages of both GDM and GSO strategies. Firstly, HBM is easy to couple with standard TO algorithms such as SIMP to efficiently obtain optimal designs involving a smaller number of design variables such as the spring stiffnesses of each grid cell. Secondly, HBM is a simple bar-spring model which does not need detail design to the internal structure and connectivity for every lattice unit cells.
The layout of this paper is as follows: Sect. 2 presents the optimization problem definition. The new version of the HBM is presented in Sect. 3. Section 4 explains how to select the value of the HBM spring stiffnesses for solving a linear elasticity problem. The implementation of the HBM for TO is described in Sect. 5. Section 6 illustrates examples of an FGL beam and an FGL plate that are optimally designed through the proposed HBM-TO method. Section 7 shows a summary of the proposed HBM-TO method and Sect. 8 gives the concluding remarks. Appendix A contains a simple python code that was used to obtain the optimal solutions presented in Sect. 6.

Problem definition
Consider an elastic FGL structure with length L , width L , and a uniform thickness h , which is illustrated in Fig. 1. A non-uniform Young's modulus E(x, y) and a constant Poisson's ratio are assumed. The FGL structure is simply supported along its longitudinal edges and is subjected to a central line load P . Body forces are not considered.
The problem at hand is to determine the optimal design of the cross-section of FGL structures composed by either large or small lattice units for maximum stiffness (Kaveh et al. 2008), i.e. having minimum compliance (Hassani and Hinton 1999) and then design a structure made by functionally graded material based on the material distribution given by the optimal FGL structures. Figure 1b and c illustrates optimal design examples that differ in terms of lattice cell size. Figure 1d illustrates an optimal design example of a functionally graded solid structure. Hencky (1921) proposed a discrete structural model comprising rigid bars connected by frictionless hinges and elastic rotational springs for solving the elastic buckling problem of columns under a compressive axial load. Since then, there has been much work done on the so-called Hencky bar-chain/net/grid model (HBM) for analysis of all kinds of structural forms from beams to frames to arches to plates (Wang et al., 2020). Recently, Zhang et al. (2021b) extended the HBM for solving plane elasticity problems which has the ability to model 2D structures for the full range of Poisson's ratio. In contrast, some previous lattice models such as the models formulated by Born and Karman (1912) and Hrennikoff (1941) can only model 2D structures for Poisson's ratio less than 1/3 (Zhang et al. 2021a;Challamel et al. 2022).

Hencky bar-grid model
In this study, we develop a slightly different version of HBM. A comparison of the new version of HBM and its predecessor reported in Zhang et al. (2021b) is given in Fig. 2. The main difference between the newly proposed HBM and its predecessor is the construction of the internal (i.e. not at the edge or corner) rigid bars. The new HBM has twice the number of rigid bars and springs. This change is essential, as it allows the determination of local stiffness matrix for each HBM grid cell (see Eq. (27)) which will be required by the topology optimization algorithm developed for the adoption of HBM. Figure 3 shows the structure representation, a unit rigid bar-grid representation and a lattice structure representation of the new HBM. As shown in Fig. 3a, the new version of HBM discretize a continuum plane structure into a rigid bar-grid system with a cell size that are connected by frictionless hinges. Within a grid cell, four primary axial springs connect two rigid bars placed at each side of the grid cell with stiffnesses k xx and k yy in the x-and y-directions, respectively, as shown in Fig. 3b. Four secondary axial springs are employed to model the Poisson effect with stiffnesses k xy and k yx in the y-and x-directions, respectively, as shown in Fig. 3b. Four torsional springs are installed at the four corners with stiffnesses k S as shown in Fig. 3b. While primary and secondary springs model the axial stiffness and the Poisson effect, respectively, torsional springs are employed to resist in-plane shear forces. A lattice structure representation related to HBM is shown in Fig. 3c. It is worth noting that the secondary axial springs of HBM are represented by some curved bars. It is possible to use other structures to represent the secondary axial springs.
For the sake of clarity, we shall demonstrate the topology optimization of HBM by using the lattice structures comprising basic lattice cell as shown in Fig. 3c. For a lattice cell whose centre is located at x = i + 1 2 and y = j + 1 2 , its horizontal and vertical axial bars have stiffnesses of k Therefore, an HBM grid with stiffer springs results in a lattice cell with a greater mass as shown in Fig. 3d. When a bar-grid cell as shown in Fig. 3b is stretched or shortened, the resulting axial forces f x , f y are resisted by both primary axial springs and secondary axial springs, i.e.
(1) can be reformulated in terms of in-plane displacements u, v,    Next, the work done by in-plane axial forces F xx , F yy , inplane shear force F xy and the in-plane body forces b x and b y lumped on the joints of each HBM cell as shown in Fig. 3b are given by The total potential energy function is given by where Ω is a set contains all the indexes of the bar-grid cells for the HBM. Now, the local stiffness matrix of an HBM cell is obtained through a variational formulation, by taking Based on Eq. (14), the discrete equations of the derivative of the elastic strain energy U i+ 1 2 ,j+ 1 2 for an HBM cell ( i + 1 2 , j + 1 2 ) are given by and Likewise, the discrete equations of the derivative of the external work W i+ 1 2 ,j+ 1 2 for an HBM cell ( i + 1 2 , j + 1 2 ) are given by In view of Eqs. (18)-(38) and by assembling the expressions for all nodes in the system, the

Determination of HBM spring stiffness
In this section, we shall determine the spring stiffnesses of HBM. For simplicity, it is assumed that each individual HBM grid cell is homogeneous and isotropic. The inhomogeneous material properties of the structure are captured by varying the material properties of different HBM cells. A general version of the continuum equations of elastodynamics (or Navier's equations of elastodynamics) considering inhomogeneous material properties and plane-strain can be written as (Gurtin 1973) and (39) 1 − e E e V e The difference equation in y-direction is obtained from Eq. (54) by swapping variables u, v and indices i, j. Now, if we assemble four adjacent HBM grid cells like Fig. 3b, the difference equation of the middle node can be obtained based on Eqs. (18)-(25) as follows: The difference equation in y-direction is Eq. (55) by swapping variables u, v and indices i, j. In view of Eqs. (50)-(55), it can be seen that solving the linear algebraic system given by the HBM proposed in this work is equivalent to solving an inhomogeneous two-dimensional elasticity problem using FD method.
Noting that the new HBM coincides with the old one (Zhang et al. 2021b) when modelling homogeneous structures. So, it is possible to use the old HBM to model inhomogeneous structures, but the resulting energy functions could not be readily applied to derive the local stiffness matrix for each grid cell.

Topology optimization
The SIMP method is employed to optimize the HBM formulated herein. The spring stiffnesses of each grid cell is formulated as In view of Eqs. (63) and (16), the total strain energy is The TO problem for HBM can be posed as where V f is the prescribed volume fraction and V(̂) and V 0 are the material volume and design volume, respectively. There are various methods for solving the aforementioned optimization problem such as Optimality Criteria (OC) (Yin and Yang 2001), Sequential Linear Programming (Dunning and Kim 2015), and Method of Moving Asymptotes (Svanberg 1987). The OC method is adopted herein due to its efficiency in solving optimization problems with a single objective function (Sigmund 1997). Following the works of Bendsøe (1989), Sigmund (1997) and Sigmund and Petersson (1998)  where δ is a positive value to limit the change of design variable ̂i + 1 2 ,j+ 1 2 between two successive iterations and γ is a numerical damping exponent. These parameters are introduced to avoid drastic change of density between adjacent cells that are impractical (i.e. porous structures). D In order to avoid potential "chess board patterns", a simple filtering technique formulated by Andreassen et al. (2011) is adopted by replacing the HBM grid cell sensitivity matrices ĉ i+ 1 2 +j+ 1 2 with where ̂(i, j, k, l) is a weight factor defined as The k and l indicate the location of corresponding joints in the model. Equations (71) and (72) indicate that the filter is a local averaging operator whereby the HBM grid cell sensitivity matrices ĉ i+ 1 2 +j+ 1 2 take on an average value of all ĉ weighted by ̂ within a radius r . The filtered HBM grid cell sensitivity matrices obtained through Eqs. (71) and (72)

Results and discussions
In this section, we first validate the proposed HBM through modeling continuous functionally graded structures. We then apply the HBM-TO to determine the optimal design of a beam and a plate FGL structures subjected to central line loading. The considered boundary conditions include loads and prescribed displacements. In the following example problems, we assume the material has Poisson's ratio of 0.3 and set the starting value for the elastic modulus E 0 as 1.

Validation of HBM for functionally graded plane bodies
In this section, plane-stress and plane-strain elasticity problems are considered that are confined in a square design domain under the loading and displacements boundary conditions illustrated in Fig. 1. The functionally graded Young's modulus is assumed to be The accuracy of the proposed HBM is assessed by comparing the predicted in-plane displacements u, v, with results obtained through the direct stiffness matrix method (Rao 2005). A general form of the local stiffness matrix for a unit square FEM plane element and a square HBM cell with unit bar length, can be expressed using Eq. (31). We assume isotropic material properties within individual FEM elements and HBM cells. The coefficients of the local stiffness matrix are for plane-stress: for plane-strain: (74) The abovementioned coefficients of the local FEM stiffness matrix can be derived as detailed in Rao (2005). For the plane-stress case, the coefficients of local HBM stiffness matrix are obtained using the HBM spring stiffnesses expression given by Zhang et al. (2021b). For the planestrain case, the coefficients of local HBM stiffness matrix are computed by using Eqs. (50)-(53). It is interesting to note that FEM and HBM have identical values for the coefficients b and d which indicate that both FEM and HBM have a similar physical interpretation of the Poisson effect.
In investigating the performance of HBM, several simulations were run by adopting grid/mesh size = L∕100 . For

Optimal design of FGL beams subjected to a central line load
The HBM-TO is first applied to the optimal design of an FGL beam with a square cross-section ( = 1 ) under a central line load. The bottom edges are fixed in y-direction. The geometry and boundary conditions are shown in Fig. 6. The HBM-TO is carried out as follows.
Noting that the variation of the elastic modulus changes the stiffness of the springs, thus the cell mass changes accordingly. In actual design of FGL structures, the crosssection area Ã , the elastic modulus Ẽ and the density ̃ of lattice bars are dependent on the applied material. Figure 7 shows the optimal cross-section designs of FGL structures obtained by HBM-TO. It can be seen that the solutions are significantly affected by the size of the HBM grid cell and the penalty exponent m . The cross-section mass is more concentrated when is smaller and m is larger. In contrast, the mass of cross-section becomes more distributed when is larger and m is smaller. It can be seen from Fig. 8 that the elastic strain energy converges after 25 iterations and the solutions converge faster when m is smaller. The converged elastic strain energy tends to be lower when m is smaller.
The optimally designed FGL beams can be classified into two groups: "distributed" and "defined"; see also Fig. 9. The first group suggests the beams could be designed with functionally graded properties such as a graded elastic modulus. The second group suggests that the beams could be designed with a defined geometry. The observed results from Fig. 9 show that the "distributed" FGL beams are stiffer than the "defined" FGL beams. The former ones have a smoother transition of mass which could help in preventing failure initiations (Wang 1983;Mahamood et al. 2012). However, the latter ones could be manufactured more readily. Based on Eqs. (56)-(60), the variable ̂ describes the distribution of the elastic modulus of the FGL beams. By least-squares fitting ̂ from HBM-TO FGL beams using = L 20 , m = 0.5 ; and = L 30 , m = 0.5 , one finds the following elastic modulus distribution function, i.e.
(87) The functionally graded beams subjected to a central line load with elastic modulus varying in x-and y-directions (as shown in Figs. 1d and 6) could be designed following the distribution function given in Eq. (86) as follows

Optimal design of FGL plate subjected to a central line load
For the second example, we apply HBM-TO to design the cross-section of an FGL plate with a rectangular cross-section under a central line load. The bottom edges are fixed in the y-direction. The geometry and boundary conditions are presented in Fig. 11. In this example, the same parameters of the previous example are used except = L 10 or L 20 in this example. Figure 12 presents the optimal cross-section designs of the FGL plates obtained from HBM-TO. Similar to what observed in the previous example, results show that the mass is more concentrated; thus, resulting in a more defined geometry when is smaller and m is larger. It can be seen from Fig. 13 that the elastic strain energy converges after 80 iterations and the solutions with smaller m converge faster. Similar to the HBM-TO designed FGL beams, the converged elastic strain energy is lower when m is smaller. The results of optimally designed FGL plates can be classified into two groups. The first group suggests the FGL plates could be designed with a more distributed elastic modulus and a smoother transition of mass which could help in preventing failure initiations (Wang 1983;Mahamood et al. 2012) as shown in Fig. 14a, while the other group (Fig. 14b) suggests the plates could be designed with a specific shape which could be manufactured more readily. By least-squares fitting ̂ from HBM-TO FGL plates using = L 10 , m = 0.5 ; and = L 20 , m = 0.5 , one finds the following elastic modulus distribution function, i.e.
where A contour plot showing the distribution given in Eq. (89) is presented in Fig. 15.

Summary
Developed herein is a novel strategy to obtain optimal designs of FGL structures. It involves the use of the Hencky bar-grid model and topology optimization. The HBM formulated herein extends a previous model with twice the number of internal rigid bars and springs as shown in Fig. 2. The new HBM enables the determination of the local stiffness matrix for each HBM grid cell by using the expressions of the stiffnesses of elastic primary axial springs, elastic secondary axial springs and torsional springs that are employed to connect the rigid bars. A TO with filtering processes is presented based on the obtained local stiffness matrices for optimizing HBM. The objective is to minimize the strain energy of an HBM subjected to boundary conditions and equilibrium and compatibility constraints by varying spring stiffnesses within each grid cell. The obtained HBM with optimal material property distribution is then used to build FGL structures (see Fig. 3c  . Therefore, HBM grid cells with larger spring stiffnesses corresponds to lattice structures with thicker bars and thus more mass. Disclaimer: The authors do not guarantee that the code is free from errors, and they shall not be liable in any event caused by the use of the code.