Yuan Liu
email: yuan.liu3@mail.mcgill.ca

Guolei Zheng

Nikita Letov
email: nikita.letov@mail.mcgill.ca

Yaoyao Fiona Zhao
email: yaoyao.zhao@mcgill.ca

Md-20-1169 Zheng

A Survey of Modeling and Optimization Methods for Multi-Scale Heterogeneous Lattice Structures

Keywords: modeling, optimization, multi-scale, heterogeneous lattice structures

This paper aims to provide a comprehensive review of the state-of-the-art modeling and optimization methods for multi-scale heterogeneous lattice structures (MSHLS) to further facilitate the more design freedom. In this survey, a design process including optimization and modeling for MSHLS is proposed. Material composition and multi-scale geometric modeling methods for representation of material and geometry information are separately discussed. Moreover, the optimization methods including multi-scale and multi-material optimization design methods, as well as the simulation methods suitable for MSHLS are respectively reviewed. Finally, the relationship, advantages and disadvantages of MSHLS modeling and optimization methods are summarized with discussion and comparison, which provides a guidance to further take advantage of MSHLS to improve the performance and multifunctional purpose of production for software developers and researchers concerning the design approaches and strategies currently available.

Introduction

Recently, Additive Manufacturing (AM) technologies, which fabricate parts by joining materials usually layer by layer, have attracted great interests from both industry and academia. Compared to other manufacturing methods, such as machining or casting, AM processes have the following unique capabilities. Firstly, parts with extremely complex shape can be built by AM process without increasing fabrication cost.

Secondly, AM technologies are suitable for processing multiple materials either simultaneously or sequentially; therefore, parts with complex material compositions can be fabricated by this manufacturing method. Thirdly, manufacturing preparation time can be substantially reduced, since the part is directly fabricated from its 3D model. These unique capabilities of AM technologies have brought great application potentials in several major industries such as aerospace [START_REF] Angrish | A critical analysis of additive manufacturing technologies for aerospace applications[END_REF] and medical implants manufacturers [START_REF] Jardini | Cranial reconstruction: 3D biomodel and custom-built implant created using additive manufacturing[END_REF]. For example, in the aerospace industry, lightweight, strong and sometimes electrically conductive parts are more desired. AM process can produce lightweight components by replacing solid material with lattice structures. Gradient electrical conductivity can also be achieved by changing the composition of materials at each fabrication point or each fabrication layer. Major airplane manufacturers such as Boeing, Airbus, and Northrop Grumman have all identified AM technologies to be the emerging and revolutionary manufacturing method [START_REF] Bourell | The roadmap for additive manufacturing and its impact[END_REF].

To take advantages of the design freedom enabled by AM processes, several Design for Additive Manufacturing (DfAM) methods have been developed [START_REF] Rosen | Computer-Aided Design for Additive Manufacturing of Cellular Structures[END_REF]. These methods aim to optimize the performance and/or other key product life-cycle considerations such as manufacturability, reliability, and cost with respect to the capabilities of AM. A recent published literature review [START_REF] Tang | A survey of the design methods for additive manufacturing to improve functional performance[END_REF] shows most existing DfAM methods such as topology optimization method or inverse homogenization. However, it only focusses on a single component design with a single type of material. This review also indicates that in order to further improve the performance of AM fabricated parts, structures with multiscale heterogeneous material distributions are the good choice.

As shown in Fig. 1, the heterogeneity referred in this paper has two aspects: heterogeneous structure (lattice cell and strut size and dimension vary) but with single material composition, and heterogeneous materials (multi material constituent in the lattice structure). A heterogeneous material is composed of different constituent materials and exhibits continuously varying composition and/ or microstructure. In this paper, we refer to lattice structures which are made of heterogeneous structures or heterogeneous materials as heterogeneous lattice structure. The heterogeneous information of such lattice structure can be described and categorized into composition and microstructure [START_REF] Kumar | An approach to modeling and representation of heterogeneous objects[END_REF]. The composition of the graded composite material can be described by the volume fractions of individual constituents that compose the material [START_REF] Li | A Review on Functionally Graded Materials and Structures via Additive Manufacturing: From Multi -Scale Design to Versatile Functional Properties[END_REF]. Heterogeneity are employed in design, often when both aspects of the resulting lattice structure would be across scale. The length-scale breakdown of multiscale structures with hierarchical architectures is shown in Fig. 2. On the macroscale ( ∼ 5 cm) (Fig. 2c), the bulk metallic structures is comprised of a network of hierarchical stretch-dominated octet unit cells which are designed to carry load via axial stress (Fig. 2d). Each stretch-dominated unit cell (Fig. 2e,f) from the hierarchical lattice network contains ∼ 200 μm hierarchical strut members (Fig. 2g) that are comprised of a network of stretch-dominated unit cells (Fig. 2h). These first-order unit cells (Fig. 2i) are comprised of microscale thin-walled hollow tube nickel-phosphorus struts with the thickness ranging from 50 to 700 nm at the lowest hierarchy (Fig. 2j).

Therefore, the multiscale lattice structure usually has the following features: the overall structure can be further divided into multiple design scales. The structures defined in the upper scale is made of the material whose microstructure is defined at the lower scale.

Comparing to those structures made of homogeneous materials or even Functional Graded Materials (FGM) [START_REF] Zhang | Additive manufacturing of functionally graded objects: a review[END_REF], multi-scale heterogeneous lattice structures (MSHLS) provides more design freedom for designers to further control the distribution of effective properties of materials on macroscale by controlling both material compositions and their microstructures. Thus, the performance of this type of structures can be further improved comparing to those structures made of conventional materials.

The results of some recent primary research [START_REF] Sivapuram | Simultaneous material and structural optimization by multiscale topology optimization[END_REF][START_REF] Zhu | Two-scale topology optimization with microstructures[END_REF] have already proven that the multiscale heterogeneous structures can exhibit better performance. This type of structures is extremely useful especially for multifunctional purposes [START_REF] Yan | Two-scale optimal design of structures with thermal insulation materials[END_REF].

Even though, the merits of MSHLS have already been unveiled, it is still difficult for designers to take full advantages of these type of structures in their product design.

There are two major barriers-modelling and optimization. Among these two major barriers, modelling is the most important, as it is the key to link the fabrication, simulation and optimization. However, most existing commercial CAD software cannot be used to model the lattice structure with heterogeneous material or structural variance.

It is mainly due to the underlying geometric modeling methods used by the software which cannot efficiently handle and manipulate the heterogeneous material or property information. To deal with this issue, several heterogeneous object modelling methods [START_REF] Sigmund | Materials with prescribed constitutive parameters: An inverse homogenization problem[END_REF] have been developed. However, most of these methods only focus on a single scale, where only the distribution of material compositions is considered. These methods cannot be used to describe the microstructure distribution of materials on multiple design scales. Moreover, the size of geometric model represented by these methods are usually large.

On one hand, MSHLS provides a great design freedom for designers. On the other hand, it also brings challenges for optimizations. There are two main challenges to optimize MSHLS. Firstly, there are a large number of design parameters that needs to be considered. These parameters are from different design scales. If designers want to solve the defined optimization directly, it usually requests high computational resource.

Secondly, the design parameters from different design scales may be coupled together, which make the optimization problem even harder to be solved.

The two major issues summarized above must be solved in order to allow designers to take further advantage of MSHLS. The aim of this review is to reflect on the MSHLS design process detailing the various phases. Actually, the current literature on the design of MSHLS are mainly focused on discussing a specific phase of the process lacking a comprehensive view. To create such a view, this work has been conceived to build a link among the main research findings available in the literature and related to the optimization design and modeling of MSHLS.

In general, this review can be beneficial for various experts working in the AM field and provides a comprehensive review of the state-of-the-art approaches to further take advantage of MSHLS to improve the performance of designed product. For example, designers could retrieve indications for better organizing their decisional process when looking for the most suitable MSHLS considering functional and manufacturing constraints. IE [14]. Within such a modeling space, every point in the base space is attached with a vector for the description of its material compositions and property, which can be symbolically described as:

( )

IE

). Without considering the differences in the definition scales of heterogeneous objects, heterogeneous objects regarded as components are analyzed, recognized and defined uniformly from three levels: integration, element and information (Fig. 3 Multimaterial-property configuration [START_REF] Liu | Material-Unit Network for Multi-Material-Property and Multiscale Components[END_REF].). Among which, the element level is subdivided into two levels, namely, the component element and material unit and the information level consists of geometric, material and property information [START_REF] Liu | Material-Unit Network for Multi-Material-Property and Multiscale Components[END_REF]. In general, heterogeneous object modeling embraces two fundamental processes: material composition modeling and geometric modeling [START_REF] Zhang | Additive Manufacturing of Functionally Graded Material Objects: A Review[END_REF]. Geometric modeling is concerned with the geometric representations of the objects and material modeling is targeted at modeling material composition distributions defined over the geometric domain. Therefore, material composition and multi-scale geometric modeling methods for representation of material and geometry information are separately discussed in Subsections 2.1 and 2.2.

Element

Component

Component Element

Material

Geometric Property

Specialized Element

Unit

Integration

Information

Material Unit

Fig. 3 Multi-material-property configuration [START_REF] Liu | Material-Unit Network for Multi-Material-Property and Multiscale Components[END_REF].

Material composition modeling

Decomposition-based models

Based on the idea of space segmentation, the space body is divided into grid elements in the decomposition-based models, such as voxel element [START_REF] Hiller | Design and analysis of digital materials for physical 3D voxel printing[END_REF][START_REF] Huang | A digital material design framework for 3D-printed heterogeneous objects[END_REF][START_REF] Leung | Digital Material Design Using Tensor-Based Error Diffusion for Additive Manufacturing[END_REF] and volume mesh element [START_REF] Wang | Finite element-based approach to modeling heterogeneous objects[END_REF][START_REF] You | Adaptive meshing for finite element analysis of heterogeneous materials[END_REF], and the material composition is defined within each element. The material information of the voxel element is stored in the voxel center point, which is suitable for the irregular distribution of the material. As shown in Fig. 4, 2D [START_REF] Andreassen | How to determine composite material properties using numerical homogenization[END_REF] and 3D [START_REF] Dong | A 149 line homogenization code for three-dimensional cellular materials written in matlab[END_REF] heterogeneous lattice structures with multi-material composition are represented by voxel-based model. It not only supports the efficient query of the material composition, but also facilitates the realization of visualization. However, the pitfall of the voxel model is that the accuracy of the method is directly related to the voxel's resolution. Moreover, voxel methods are inexact in terms of the geometric and material accuracies. On the contrary, the material information of the volume mesh-based model is stored in the grid node, and the material information of the non-grid node is saved by interpolation function. Therefore, the representation is a more compact data structure, which solves the huge storage problem in voxel unit to a certain extent. Representation of material composition using these decomposition-based depends on the resolution of the mesh or voxels that may not conform to the material distribution, leading to discretization error. In addition, any change in material function would lead to rediscretization of the whole geometry to re approximate the new material distribution.

Set-based models

According to the geometric representation, set-based models are divided into two types. The first one is rm-set, which is an extension of solid modeling to handle heterogeneous objects by using r-sets (a form of B-rep modeling) with information about material variation. This geometry model is further decomposed into atlases over which material variations are mapped. In order to define FGM models, a set of Boolean operators [START_REF] Sun | Reasoning Boolean operation based modeling for heterogeneous objects[END_REF] for modifying the geometry and composition stored in the data structure, therefore, the heterogeneous objects with more complex geometry and material distribution are obtained. Another set-based model is based on implicit function such as level-set [START_REF] Sethian | Structural Boundary Design via Level Set and Immersed Interface Methods[END_REF], R function [START_REF] Biswas | Heterogeneous material modeling with distance fields[END_REF] and convolution surfaces [START_REF] Gupta | Heterogeneous object modeling with material convolution surfaces[END_REF] to generate exact geometric data representations. These models use the functional representation (F-Rep) as the basic model for both the point set geometry and the material distribution.

Essentially, the implicit function is used to model spatial partitions and the function determines whether a given point is "inside" or "outside" a subset. When a point is within the subset, the corresponding material distribution is then assigned to that point/point set. Tsukanov and Shapiro [START_REF] Tsukanov | Meshfree modeling and analysis of physical fields in heterogeneous media[END_REF] proposed the construction of field modeling using sampling distance field and interpolated physics field by way of extended Rfunctions, however, the limitations of R-functions results in the difficulties in complex shaped heterogeneous object modeling. A level set based scheme was proposed by Wang et al. [START_REF] Wang | A level-set based variational method for design and optimization of heterogeneous objects[END_REF], which adapted a variational model as the objective function to locate any point in the material region of a well-defined gradient or on the boundary edges and surfaces of discontinuities. The set of discontinuities was represented implicitly, using a multi-phase level set model. The compact and exact in data representation capacity makes level set model suitable for multi-material structural representation. In addition, parameterized level set method can be used for CSG (Constructive Solid Geometry) type feature-based CAD modeling. The rectangular frame and the interior hole are the two zero-value level set contours, i.e. the material source profiles. All results demonstrated in Fig. 5 have realized the smooth material transitions between the source profiles [START_REF] Liu | Level set-based heterogeneous object modeling and optimization[END_REF].

Control point-based models

The material information is stored in the control points of Bezier, B spline, NURBS, parameter curve, surface and volume, and the representation method of parameter control points is obtained [START_REF] Yang | A B-spline-based approach to heterogeneous objects design and analysis[END_REF][START_REF] Samanta | Optimized normal and distance matching for heterogeneous object modeling[END_REF][START_REF] Sasaki | Adaptive direct slicing of volumetric attribute data represented by trivariate B-spline functions[END_REF][START_REF] Yoo | Heterogeneous object modeling using the radial basis functions[END_REF]. Control point-based models have some appealing properties. They are compact in both geometry and material representations.

Given the parametric coordinate ( ) , , u v w , the material composition of a point can be also efficiently interrogated. In addition, this representation method can effectively represent complex (2D or 3D dependent) material distributions. Local modifications on both geometries and material definitions are also straightforward. The drawback of this method is that it relies heavily on spatial parameterizations and for arbitrary 3D objects such parameterizations remain a rather non-trivial task. Kou et al. [START_REF] Kou | A hierarchical representation for heterogeneous object modeling[END_REF] proposed a more general definition for 1D heterogeneous features (Fig. 8) described in Eq. (2). (2

)
where P is an arbitrary point on the 1D feature, i P is constructive point (e.g., starting point, ending point of the line and the control points of the B-Spline curve), i W

is the blending weight of the ith constructive point in material gradation. For different types of 1D features, the weight generation methods might be different. Similarly, these approaches could be extended to 2D and 3D features. To more clearly compare the methods that investigate the representation abilities of material heterogeneity, both the advantages and disadvantages of existing models for heterogeneous materials are summarized in Table 1. It serves as a reference tool for selecting appropriate models in the downstream applications. pervasive, it is foreseen that some materials modeling modules will be embedded in future CAD systems so that the geometry and materials of a product can be designed concurrently [START_REF] Huang | A multiscale materials modeling method with seamless zooming capability based on surfacelets[END_REF]. Rosen et al. [START_REF] Huang | Material feature representation and identification with composite surfacelets[END_REF][START_REF] Jeong | Microstructure Feature Recognition for Materials Using Surfacelet-Based Methods for Computer-Aided Design-Material Integration[END_REF][START_REF] Wang | Multiscale Heterogeneous Modeling with Surfacelets[END_REF] firstly proposed a new multi-scale geometric and materials modeling methods that uses a surfacelet-based implicit representation to efficiently capture internal and boundary information of materials. It serves as the foundation for modeling structure-property relationships for materials design.

Therefore, multi-scale geometric modeling methods are reviewed in the following section.

Multi-scale geometric modeling 2.2.1 Multi-scale problem

Before approaching the topic of multi-scale geometric modeling, it is important to cover issues emerging from multi-scale mathematical models describing reality. These issues come from the fact that on different scales different laws of physics act and different mathematical models are used to describe what happens on these scales and how different entities react to each other on these scales. The main challenge of multi- scale modeling is the lack of a unified theory that could correctly describe geometry both down from its structure (~10 -9 m) up to boundary conditions and processes present in it (~10 0 m) [START_REF] Raabe | Multi-scale modeling in materials science and engineering[END_REF]. A common example of a multi-scale problem comes from theoretical physics and lies in finding a theory of everything -a hypothetical ultimate framework that would allow describe physical relations on every possible scale.

Nowadays, in physics there are two major models describing our world: quantum field theory (comes from merging quantum mechanics with special relativity) and general relativity (comes from merging Newtonian mechanics with special relativity); there is nothing similar between these models as the first one works well for nano-and lower scales while the second one works for macro-and higher scales [START_REF] Schombert | Quantum gravity[END_REF].

The ability of modeling in meso-and macro-scale is crucial in many research fields involving computer graphics and geometric modeling. For example, Prada et al. [START_REF] Prada | Superficial 3D mesh generation process using multimedia software for multiscale bone analysis[END_REF] apply multi-scale modeling concepts for 3D mesh generation process for multiscale bone analysis. Engquist et al. [START_REF] Engquist | Multiscale methods in science and engineering[END_REF] identified the need for an appropriate multi-scale geometric modeling tool for porous materials. Predictive nonlinear theories are being applied for multiscale modeling of heterogeneous materials [START_REF] Matouš | A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials[END_REF].

Multi-scale geometric modeling methods are also critical in medical research. For example, Youssef [START_REF] Youssef | Parallelization of a bio-inspired computational model for the simulation of 3-D multicellular tissue growth[END_REF] introduced a bio-inspired framework for the simulation of 3D multicellular tissue growth. However, there is no geometric modeling tool to implement this framework due to high geometric complexity. The need for such a geometric modeling approach for modeling a human heart both as a whole and in details (as even small defects are crucial) is identified [START_REF] Sacks | On the need for multi-scale geometric modelling of the mitral heart valve[END_REF]. Another common application for multi-scale geometric modeling is visualization of online geographic maps: when a user zooms in to a map, e.g. Google Maps, more and more details appear to the extent that even 3D shapes of buildings render [START_REF]Developer Guide: Maps Static API[END_REF].

Geometric modeling of heterogeneous lattice structures

Geometric modeling in engineering is mainly used for visualizing engineering systems such as parts, assemblies, etc., and represented by computer-aided design (CAD) software. While conventional CAD software is normally sufficient for engineering design in industry, there is no geometric modeling tool that would be able to appropriately represent complex heterogeneous lattice structures due to inability of existing geometric modeling tools to represent a model in multi-scale, which forms a research gap that is yet to be filled. Another issue is the computational cost of such geometric modeling tools because heterogeneous lattice structures have high geometrical complexity resulting in high computational needs using conventional methods [START_REF] Dong | A survey of modeling of lattice structures fabricated by additive manufacturing[END_REF].

Multi-scale modeling becomes a crucial add-on to conventional geometric modeling tools when dealing with lattice structures. It is required to consider not only the whole structure in its macro-scale, but also each joint and strut of the lattice as they form at meso-scale [START_REF] Zhang | Extended multiscale finite element method for mechanical analysis of heterogeneous materials[END_REF].

Level of detail

In computer graphics, a far-away object in a virtual environment appears obscured (similar to the effect people with myopia experience), which means that there is no need for representing all details of the object on a screen [START_REF] Marschner | Fundamentals of computer graphics[END_REF]. However, the details appear to be more concrete when approaching the object or zooming into it. Such effect of manipulating the level of detail (LOD) can be achieved, for example, by reducing the polygon count for distant objects by vertex and edge removal [START_REF] Luebke | Level of detail for 3D graphics[END_REF]. However, decreasing LOD of a lattice structure down to a certain level can result in complete homogenization of lattice, i.e. removing from view all geometrical features corresponding to lattice structures [START_REF] Ptc | Material homogenization for lattice simulation in additive manufacturing[END_REF].

Normally, LOD is manually or automatically associated with CAD-features of solid bodies as illustrated in Fig. 9, e.g. extrusions, revolutions, etc. [START_REF] Borrmann | Multi-scale geometric-semantic modeling of shield tunnels for GIS and BIM applications[END_REF]. However, lattice structures are designed not with CAD-features but with nodes connected by struts which are not so well-defined from the CAD-perspective [START_REF] Vayre | Designing for additive manufacturing[END_REF]. Moreover, in heterogeneous lattice even periodicity of the lattice is not obvious. As mentioned, simplifying a lattice structure by reducing the LOD of its model eventually results in complete homogenization of the structure. Such model does not provide the user with sufficient information about structure and topology.

Another way to simplify a lattice structure is to represent it as a beam-element model. This approach, however, avoids representing detailed geometrical information that can be important for an appropriate analysis of a lattice structure, especially in node regions as beam models simply do not take high material concertation at nodes into the account [START_REF] Tang | A survey of the design methods for additive manufacturing to improve functional performance[END_REF].

Setting up a limit for the maximum distance for an object to be shown on screen is a common practice in computer graphics. This limit is called a draw distance [START_REF] Connell | 3D for Graphic Designers[END_REF]. It is mainly used for optimizing the overall performance of a 3D model [START_REF] Limper | The pop buffer: Rapid progressive clustering by geometry quantization[END_REF].

The described model complexity reduction methods have proven their effectiveness in the video game industry when developing game engines, such as

AnvilNext developed by Ubisoft Montreal, in particular for rendering 2D or 3D graphics and providing a physics engine [START_REF] Gregory | Game engine architecture[END_REF]. Some other geometric modelling approaches also can be adapted from game engines that are being used for videogames design.

Considering applying the LOD concept and multi-scale modeling to lattice structures, Liu and Shapiro [START_REF] Liu | Multiscale shape-material modeling by composition[END_REF] propose a framework for multi-scale modeling of parts with lattice structures. The framework's first scale describes the shape-material model;

the second scale describes identity mapping for the solid regions of the part, as well as it includes downscaling and neighborhood functions for lattice structures; the third scale includes downscaling and neighborhood functions for lattice structures in case beams of the lattice structures are made of the fine scale structures. Such division into scales and phases essentially corresponds to the LOD concept, as the complexity of the solid model iteratively rises with downscaling.

Voxel modeling

One more topic worth considering when discussing multi-scale geometric modelling is application of multi-scale voxel modeling. The word 'voxel' derives from combining words 'volumetric' and 'pixel' and it is defined as 'a unit of volume containing the value of the corresponding raster element in a 3D space' [START_REF] Kaufman | Volume graphics[END_REF]. Voxelized objects allow the same set of operations with polygonal objects [START_REF] Cohen-Or | Fundamentals of surface voxelization[END_REF]. Voxel-based object simplification has been used for eliminating high-frequency details of the object ever since the introduction of voxels [START_REF] He | Voxel based object simplification[END_REF]. Similar voxel-based approaches are used for simplification and repair of polygonal models [START_REF] Nooruddin | Simplification and repair of polygonal models using volumetric techniques[END_REF]. Voxels have an advantage in terms of downsampling and acquisition of real-world data, i.e. any geometrical complexity is feasible which, however, comes at a cost due to rise of computational resources required for rendering [START_REF] Laine | Efficient sparse voxel octrees[END_REF]. Note that as lattice structures cannot be produced without AM techniques, there is no need in voxel size higher than 3D-printer resolution [START_REF] Telea | Voxel-based assessment of printability of 3D shapes[END_REF]. Voxel models do not handle zooming efficiently [START_REF] Laine | Efficient sparse voxel octrees[END_REF].

One of the most popular voxel-based simplification methods involves using sparse voxel octrees which are based on generating multi-scale voxels which could be visible or invisible depending on the resolution and size of the screen. Moreover, recent researches show that octree-based neural networks can be applied for 3D shape analysis and learning and thus can be adapted for feature recognition. For example, Wang et al.

introduced an octree-based convolutional neural networks (CNN) for 3D shape analysis [START_REF] Wang | O-cnn: Octree-based convolutional neural networks for 3d shape analysis[END_REF]. Liu et al. [START_REF] Liu | Learning a hierarchical latent-variable model of 3d shapes[END_REF] proposed an octree-based variational shape learning model for the same purpose. Wu et al. [START_REF] Wu | 3d shapenets: A deep representation for volumetric shapes[END_REF] developed a Bayesian volumetric shape recognition based on voxel octrees.

Recent research on voxel-based surface approximation investigates the possibility of applying the previously discussed LOD concept to voxel models. The results have proven themselves encouraging, as a 12 GB voxel model can be approximated to a 2 GB while still accurately visualizing the details [START_REF] Marcus | Level-of-Detail Independent Voxel-Based Surface Approximations[END_REF]. However, in this method, the depth of rendering (which is required for fully adaptive multi-scale voxel modeling [START_REF] Tian | Adaptive voxels: interactive rendering of massive 3D models[END_REF]) is assumed to be given and the voxelization algorithm is not adaptive, which does not allow performing multi-scale modeling. Moreover, current approaches are unable to represent crucial features of a part on a larger scale with voxels [START_REF] Seemann | Simplification of Multi-Scale Geometry using Adaptive Curvature Fields[END_REF] and become significantly slow on a smaller scale [START_REF] Kauker | VoxLink-Combining sparse volumetric data and geometry for efficient rendering[END_REF].

Lattice structure is essentially a structure of some size that is defined by a topology of a significantly smaller size. Thus, voxelization of lattices introduces visible discreteness of its model as illustrated in Fig. 10, which affects user experience. Voxels normally have cubic shape which is explained by render simplicity: one cubic voxel can be defined only by its position in 3D space (assuming that every voxel has the same size). Cubic shape also fits well for simulations, e.g. a cube is a common unit element for stress analysis within structures as noted in Fig. 11a since it easily allows modeling of normal stresses 𝜎𝜎 𝑥𝑥 , 𝜎𝜎 𝑦𝑦 and 𝜎𝜎 𝑧𝑧 , as well as shear stresses 𝜏𝜏 𝑥𝑥𝑦𝑦 , 𝜏𝜏 𝑥𝑥𝑧𝑧 , 𝜏𝜏 𝑦𝑦𝑧𝑧 , 𝜏𝜏 𝑦𝑦𝑥𝑥 , 𝜏𝜏 𝑦𝑦𝑧𝑧 and 𝜏𝜏 𝑧𝑧𝑦𝑦 . However, representing cylindrical shapes illustrated in Fig. 11b as a combination of cubes introduces an obvious distortion which gets lower as the number of cubes approaches infinity, which is not feasible from the computational point of view. Similarly, not cubic but cylindrical sector unit elements are used for simulations of cylindrical shapes, which allows modeling of normal stresses 𝜎𝜎 𝑟𝑟𝑟𝑟 , 𝜎𝜎 𝜃𝜃𝜃𝜃 and 𝜎𝜎 𝑧𝑧𝑧𝑧 , as well as shear stresses 𝜎𝜎 𝜃𝜃𝑧𝑧 and 𝜎𝜎 𝑧𝑧𝜃𝜃 [START_REF] Xu | A three-dimensional soil-water coupled FE analysis of hollow cylinder test concerning non-uniform deformation[END_REF]. Strand [START_REF] Strand | Surface skeletons in grids with non-cubic voxels[END_REF] conveyed research in the area of using non-cubic voxels and tested various voxel shapes such as body-centered cubic (BCC) and face-centered cubic (FCC)

(a) (b)
Note that BCC is similar to a truncated octahedron and FCC is similar to a rhombic dodecahedron. Non-cubic grids appear more sparse than cubic grids which provide the most information about the structure. However, there is no evidence of research that had been conveyed on combining different shapes of voxels within one model. Such research would contribute to not only AM of heterogeneous lattice structures but also to 4D printing technique, i.e. AM with smart materials [START_REF] Sossou | Design for 4D printing: A voxel-based modeling and simulation of smart materials[END_REF].

There is evidence that volumetric meshes can support multi-scale modeling.

Volumetric meshes are similar to polygon meshes that are used in the majority of CADsoftware with the difference in discretizing the whole solid body and not only its surface, i.e. the whole body is considered to be subdivided with polyhedrons instead of polygons [START_REF] Rom | Volume Mesh Generation for Numerical Flow Simulations Using Catmull-Clark and Surface Approximation Methods[END_REF]. For instance, volumetric mesh models of macromolecules are able to provide sufficient visual information in chemistry [START_REF] Feng | Multiscale geometric modeling of macromolecules II: Lagrangian representation[END_REF]. However, the computation of curvature values of polygon and polyhedron meshes is non-trivial due to their discrete nature.

Moreover, a volumetric finite element is always convex since it is bounded by a finite number of planes. Therefore, FEV introduces even more distortion when applied to non-convex geometric shapes such as lattice structures and surfaces with genus 1 and higher [START_REF] Tewari | Meshing genus-1 point clouds using discrete one-forms[END_REF]. This implies the need in a proper meshing algorithm that takes convexity and curvature into the account and affects quality of meshed models, especially ones requiring multi-scale modeling. Note that non-convex volume selection is a topic of interest in geometric modeling [START_REF] Fuchs | Non-convex polyhedral volume of interest selection[END_REF][START_REF] Gortler | Discrete one-forms on meshes and applications to 3D mesh parameterization[END_REF].

Table 2 compares polygon mesh methods with voxel-based method when applied to MSHLS. Note that the table does not include F-Rep method as it is incompatible with other modeling formats and cannot store topology information [START_REF] Pasko | Procedural function-based modelling of volumetric microstructures[END_REF]. The hybrid modeling approach introduced by Tang et al. [START_REF] Tang | A hybrid geometric modeling method for lattice structures fabricated by additive manufacturing[END_REF] has certain features (e.g. using boundary representation for lattice struts) that can assuage the disadvantages of both geometric modeling methods described in Table 2.

Note that as the size of a lattice structure approaches nano-scale, new properties arise, some of which are of quantum nature [START_REF] Bauer | Nanolattices: an emerging class of mechanical metamaterials[END_REF]. As there is no tool that can efficiently produce a lattice structure with tolerance down to its nanoscale, this work does not support geometric modeling of lattice structures at nanoscale.

Optimization methods of heterogeneous lattice structures

Lattice structures are often produced to solve engineering optimization problems and AM itself prompts research in computational optimization [START_REF] Tang | A hybrid geometric modeling method for lattice structures fabricated by additive manufacturing[END_REF]. In heterogeneous lattice structures, optimization of structural topology and thickness of lattice struts results in satisfaction of given functional requirements [START_REF] Tang | Bidirectional Evolutionary Structural Optimization (BESO) based design method for lattice structure to be fabricated by additive manufacturing[END_REF]. Heterogeneity of a structure provides unique properties which are required to be optimal and which cannot be provided by homogeneous structures, e.g. gradual elasticity of the structure or structural stiffness [START_REF] Martínez | Procedural voronoi foams for additive manufacturing[END_REF]. Research on optimizing the internal structure to provide an optimal strength-to-weight ratio is identified as crucial, for example, in the aerospace sector [START_REF] Vasiliev | Anisogrid composite lattice structures-Development and aerospace applications[END_REF]. Topology optimization methods are often being applied for lattice structures [START_REF] Liu | Current and future trends in topology optimization for additive manufacturing[END_REF]. For example, Zhu et al. [START_REF] Zhu | Two-scale topology optimization with microstructures[END_REF] proposed a two-scale topology optimization framework of microstructures. One of the ways to optimize a lattice structure is by adding porosity to the internal structure. Moreover, the approaches that are being used for the porosity generation (such as the ones that are based on Voronoi tessellation [START_REF] Lu | Build-to-last: strength to weight 3D printed objects[END_REF]) produce a tessellation that is visually similar to the tessellation of cells in a living organism. Based on this, it is crucial to analyze multi-scale and multi-material optimization methods of heterogeneous lattice structures covered in Subsection 3.2 and Applying the LOD concept is not straightforward and does not achieve as good results as for polygon mesh modeling 3.3, respectively, as well as the ability to perform computational simulation on them as covered in Subsection 3.1.

Simulation methods of multi-scale heterogeneous lattice structures

This section of the design workflow is the evaluation of the MSHLS performance through simulation methods to support optimization design.

The simplest way is to directly apply the rules of mixtures, where the property of heterogeneous composite can be considered as a linear or non-linear combination of its base materials [START_REF] Kim | On the rule of mixtures for the hardness of particle reinforced composites[END_REF][START_REF] Sola | Functionally graded materials for orthopedic applications -an update on design and manufacturing[END_REF]. Comparing to other ways to evaluate the effective properties of the heterogeneous objects, the rule of mixtures is the simplest and most straightforward method. The effective properties X of a point with material composition values ( , ,..., ) n X X X could be calculated by Eq. (3).

0 1 1 2 2 ... n n X X X X X ρ ρ ρ ρ = = + + + (3)
However, rules of mixtures are not accurate and flexible enough for MSHLS.

Some numerical analysis methods suitable for MSHLS including physics-based method, homogenization method, finite element analysis (FEA) method, and isogeometric analysis (IGA) are reviewed as follows.

Physics-based method

There are many established methods and implementations for simulating deformable soft bodies. Considering the large deformations and relatively low stiffness of the materials involved, the physics-based dynamics are often significant. When our proposed MSHLS consist of multiple materials, it's important to consider deformation simulation for further application such as feature films, video games, and virtual surgery, among others. Much of the development in simulating soft bodies has been driven by the computer graphics community. In this method each voxel is modeled as a lattice point with mass and rotational inertia. Voxels are connected (from their centers) by three-dimensional beam elements with appropriate translational and rotational stiffness leading to realistic deformation under applied forces and moments. These beams form a 3D lattice frame, which acts as a backbone (or a control structure) of the whole shape.

The frame is what holds the multi-material objects together and that governs the deformation of the whole geometry. Based on this method Sossou G et al. [START_REF] Sossou | Design for 4D printing: A voxel-based modeling and simulation of smart materials[END_REF][START_REF] Sossou | Design for 4D printing: Modeling and computation of smart materials distributions[END_REF] developed a computational design tool, which has been implemented in the Rhinoceros® add-on Grasshopper®(GH).

Homogenization method

The homogenization usually refers to a way to replace the composite with a kind of equivalent material model, which can overcome the difficulty in the analysis of the boundary value problem with high heterogeneities. It is used to obtain the effective properties of homogenized material for periodic heterogeneous continuous media in many physical and engineering applications [START_REF] Dong | A survey of modeling of lattice structures fabricated by additive manufacturing[END_REF]. Homogenization methods can also be applied to calculate the effective properties of the heterogeneous objects at any points [START_REF] Allaire | Homogenization of a conductive, convective, and radiative heat transfer problem in a heterogeneous domain[END_REF][START_REF] Hassani | A review of homogenization and topology optimization I -Homogenization theory for media with periodic structure[END_REF]. It can be applied to the cells with any arbitrary material compositions and can also be applied for the calculation of several different types of material effective properties including stiffness, heat transfer coefficients and Poisson's ratio. This method is applied for material effective properties simulation of lattice cell with geometric heterogeneity [START_REF] Xu | Design of lattice structures with controlled anisotropy[END_REF] but also material composition heterogeneity [START_REF] Andreassen | How to determine composite material properties using numerical homogenization[END_REF][START_REF] Dong | A 149 line homogenization code for three-dimensional cellular materials written in matlab[END_REF].

Based on the calculated effective properties of the lattice unit, the overall performance of designed objects can be easily evaluated. A typical procedure of simulation of solidlattice hybrid structure with homogenization method has been introduced by Dong [START_REF] Dong | Design and optimization of solid lattice hybrid structures fabricated by additive manufacturing[END_REF].

However, this method has its own limitations. It cannot be applied to the structure based on heterogeneous density gradients, and it is not easy to manage in case of complex geometries.

FEA-based method

FEA-based simulation method is usually used in interactive inverse design method to create continuous heterogeneous material distributions on 1D beam [START_REF] Yu | A Computational Method for the Design of an Additively Manufactured Personalized Artificial Spinal Disc With Physiological Stiffness Under Multiple Loading Conditions[END_REF], truss element with truss-like lattice structure [START_REF] Lumpe | Computational Design of 4D Printed Shape Morphing Multi-State Lattice Structures[END_REF] or 3D volumetric meshes with arbitrary geometry [START_REF] Xu | Interactive material design using model reduction[END_REF][START_REF] Bickel | Design and fabrication of materials with desired deformation behavior[END_REF][START_REF] Bickel | Capture and modeling of non-linear heterogeneous soft tissue[END_REF]. The material distributions are designed to conform to prescribed displacements and internal elastic forces at the selected vertices of the 1D element or 3D mesh [START_REF] Fang | Geometry-based direct simulation for multi-material soft robots[END_REF]. Such a capability has numerous applications as many mechanical components, structures, and mechanisms have to produce predictable displacements under known forces or pressure distributions. With smooth varying heterogeneous material distributions, multiple properties (functions) can be obtained in the same object, and high-stress regions at the material boundaries can be avoided. The method is applicable not only to small deformations, linear mesh elements but also to nonlinear materials under large deformations. Each mesh element is assumed to be homogeneous and made of an isotropic material parameterized by Young's modulus and Poisson's ratio. Two applications of FEA-based simulation aided design with 1D truss element and 3D tetrahedral mesh, respectively are shown in Fig. 12. FEA method with 1D element has a moderate computational cost and can be used for MSHLS with heterogeneous density gradients. Its limitations are the following: when the strut is not thin enough, it is difficult to provide accurate results and the ability to predict the stress distribution at struts joints is poor. While, FEA method with 3D mesh is of high computation cost and time consuming.

IGA-based method

A common approach to the analysis and optimization of metamaterials and lattices is multiscale simulation through homogenization of unit cell behavior. However, due to the functional grading and curved topology, which lead to non-periodicity of the microstructure, as well as the nonlinearity of a soft lattice, common homogenization approaches cannot be applied, and the lattice needs to be simulated at full scale. In addition, nonlinear simulation based on continuum mechanics and 3D finite elements is not feasible due to the huge computational effort involved and linear modelling with 3D truss elements will not capture the nonlinear behavior of a soft lattice [START_REF] Weeger | Isogeometric shape optimization of nonlinear, curved 3D beams and beam structures[END_REF]. While, IGA uses B-spline function and NURBS basis function, which are commonly used in CAD environment, to accurately describe any complex geometry and approximate unknown solution fields in the analysis [START_REF] Tavakkoli | Isogeometric topology optimization by using optiMality criteria and implicit function[END_REF]. It provides an accurate and efficient numerical discretization of the model and enables a seamless integration of the simulation method with the design approach through the concept of IGA. The rod formulation and IGA method also allow to accurately represent graded design parameters of the cross-sections, such as a spatially variable material distribution, e.g.

in terms of Young's modulus. These properties are highly desirable for many types of applications, such as reusable energy absorbing devices based on nonlinear elastic rather than plastic deformations, vibration mitigation through curved ligaments, tailored energy absorption response through functional grading of material properties or hierarchical microstructures, tissue-like medical implants, soft robotic actuators and devices, 4D printing, or structures with complex constitutive behavior such as negative Poisson's ratio effects. Weeger O et al. [START_REF] Weeger | Digital design and nonlinear simulation for additive manufacturing of soft lattice structures[END_REF] proposed a digital design and manufacturing framework for soft lattice structures (Fig. 13). The design approach was implemented in the 3D modeling software Rhinoceros® (also Rhino® or Rhino3D®) using the algorithmic design environment Grasshopper™ and the open-source plugin IntraLattice. Based on the simulation methods introduced above, a comparison between different simulation approaches is made and summarized in Table 3, and it helps designers select a suitable simulation method for MSHLS.

Multi-scale optimization methods of heterogeneous lattice structures

As to the optimization of the heterogeneous lattice structures, the existing methods can be roughly divided into three groups according to their design scales. It mainly includes the optimization design of microstructure, optimal distribution of material compositions or microstructures in macroscale and the two-scale concurrent optimization.

Optimization design of microstructures

The inverse homogenization technique can be applied to determine the optimal distribution of material in micro-cell at a certain point of the microstructures to achieve their desired properties [START_REF] Gao | Topology optimization of micro-structured materials featured with the specific mechanical properties[END_REF]. In this process, a microstructure periodically arranged in

is smooth and continuous

The stiffness matrix of higher order basis function is dense and does not satisfy the orthogonal condition macro design domain actually behaves like a material, whose length scale is much smaller than that of the macroscopic structure. This method has originally been developed for the structural related properties, such as extreme elastic modulus, Poisson's ratio, and extreme thermal expansion coefficients. Also, it can be easily extended to many other fields for multifunctional purposes. Li et al. [START_REF] Li | Integrated design of cellular composites using a level-set topology optimization method[END_REF] produced periodic cellular composite for each layer, by integrating the numerical homogenization into a level set approach (Fig. 14a). Kazemi et al. [START_REF] Kazemi | Multi-material topology optimization of lattice structures using geometry projection[END_REF] proposed a computational method for the design of architected truss lattice materials where each strut can be made of one of a set of available materials (Fig. 14b). Schumacher et al. [START_REF] Schumacher | Microstructures to control elasticity in 3D printing[END_REF] solved the inverse problem to the numerical coarsening method, obtaining a microstructure that coarsens to a given stiffness tensor (Fig. 14c). The difference among the three methodology is that the representations of microstructures are continuous method (Fig. 14a) and discrete method, in which the discrete method uses 1D truss element (Fig. 14b) and voxel element (Fig. 14c), respectively.

Optimal distribution of material compositions or microstructures in macroscale

Optimal distribution of material compositions. Instead of only focusing on microscale, the second group of optimization methods is developed to determine the optimal distribution of material compositions or microstructures. Due to the limitations of traditional manufacturing technologies, most of the existing optimization methods mainly focus on the optimization of material compositions distributions. During the optimization of material compositions, the material distribution functions are usually preselected. Then, the optimization formulation is generated to find the optimal coefficients which are defined in the preselected functions. For example, in the optimization of Functional Graded Materials (FGM) which is a typical heterogeneous object, power law and exponential functions are widely used [START_REF] Sola | Functionally graded materials for orthopedic applications -an update on design and manufacturing[END_REF][START_REF] Elishakoff | Three-dimensional analysis of an all-round clamped plate made of functionally graded materials[END_REF][START_REF] Na | Volume fraction optimization of functionally graded composite panels for stress reduction and critical temperature[END_REF]. Instead of using functions, SIMP-based topology optimization results with relative density distribution have also been used for the design of the heterogeneous lattice structure.

For example, Brackett et al. [START_REF] Brackett | Topology optimization for additive manufacturing[END_REF] suggested mapping volume fractions of the lattice unit cells onto the intermediate densities of an un-penalized SIMP solution, making the greyscale density solution of topology optimization possible to manufacture with AM.

This work utilized the tessellation of the unit cell for the lattice generation, i.e. where a selected unit cell template was tessellated across the design domain in a regular fashion.

To generate lattice with varying cell sizes, Brackett et al. [START_REF] Brackett | An error diffusion based method to generate functionally graded cellular structures[END_REF] offered an error diffusion-based method which enabled the mapping of irregular unit cell lattice onto a grey scale input. This was achieved by generating dithered points from the greyscale image followed by connecting them, using either Delaunay triangulation or Voronoi tessellation. Tang et al. [START_REF] Tang | Multifunctional design of heterogeneous cellular structures[END_REF] developed a function-performance-property-design parameter model (F-P-P-D) (Fig. 15) for heterogeneous lattice structures where the result of density-based topology optimization has been mapped to the relative density of lattice structures. Its result shows the optimized heterogeneous lattice structures can significantly achieve both low thermal conductivity and high stiffness without increasing its weight. Alzahrani et al. [START_REF] Alzahrani | Design of truss-like cellular structures using relative density mapping method[END_REF] proposed method utilizes the relative density information (Fig. 16) obtained from a solid topology optimization to automatically determine the diameter of each individual strut in the structure, which collectively represents the set of design variables. This allows the method to produce lattice structures that can perform reliably under multiple loading conditions and also reduce the computational cost associated with the design of these structures. Li et al. [START_REF] Li | Optimal design and modeling of gyroid-based functionally graded cellular structures for additive manufacturing[END_REF] presented a novel optimization strategy for designing functionally graded cellular structures with desired mechanical properties. In the strategy a complex functionally graded gyroid lattice structure is generated by implicit surface reconstruction algorithm in 3D space, which maps the corresponding material properties (Fig. 17a) to geometric parameter (Fig. 17b) perfectly. Panesar et al. [START_REF] Panesar | Strategies for functionally graded lattice structures derived using topology optimisation for additive manufacturing[END_REF] presented a number of strategies that enable lattice structures to be derived from topology optimization results suitable for AM. As shown in Fig. 18, the structures with different strategies such as solid (Fig. 18a), intersected lattice of D-P (Fig. 18b), intersected lattice of BCC (Fig. 18c), graded lattice of D-P (Fig. 18d), scaled lattice of D-P (Fig. 18e), and uniform lattice of D-P (Fig. 18f) are used for FEA and made a comparison for several objectives. Instead of using relative density distribution results, stress distribution results of the structures could also be used for optimization. Teufelhart and Reinhart [START_REF] Teufelhart | Optimization of strut diameters in lattice structures[END_REF] optimized the strut diameters of an irregular lattice structure by capitalizing on the flux of force within a solid structure. They showed that enhancement in performance can be achieved when compared to a regular (uniformly tessellated) counterpart, mainly due to a much smoother distribution of stresses owing to the stress conformal nature of the lattice. Also Tang et al. [START_REF] Tang | Bidirectional Evolutionary Structural Optimization (BESO) based design method for lattice structure to be fabricated by additive manufacturing[END_REF] proposed a BESO-based design method for truss-like lattice structure optimization. In this method, Functional Volumes (FVs) and Functional Surfaces (FSs) are first determined based on an analysis of the functional requirements.

FVs can be further decomposed into several sub-FVs. These sub-FVs can be divided into two types: FV with solid and FV with lattice.

Optimal distribution of microstructures. Besides material compositions, the microstructure distribution is another type of design parameter which are considered during the optimization process. For instance, to achieve the given elasticity gradients, an inverse-homogenization technique-based design optimization is proposed for functional graded cellular structures [120]. In this method, a precomputed data-based of tiled cell structures is built to cover a wide range of elastic properties. Then, a global optimization algorithm is proposed to synthesize the cells from the pre-computed database to achieve desired elasticity distribution with the consideration of connections between neighborhood cells. To simultaneously update the topology of unit cell as well as the morphology of cells in the different region of design domain, Liu et al. proposed a new approach to generate functional graded cellular structures based on moving morphable components/voids (MMV/MMC) topology optimization [START_REF] Guo | Doing topology optimization explicitly and geometrically-a new moving morphable components based framework[END_REF][START_REF] Liu | Additive Manufacturing-Oriented Design of Graded Lattice Structures Through Explicit Topology Optimization[END_REF]. In this method, the topology of the unit cell is represented by a set of explicit functions, while the basis perturbation functions have been applied to control the morphology of cell distribution. This innovative design parameterization method enables designers to control both cell topology and its morphology distribution simultaneously without a large number of design parameters. Similarly, Schumacher et al. [START_REF] Schumacher | Microstructures to control elasticity in 3D printing[END_REF] proposed a method for fabricating deformable objects with spatially varying elasticity using 3D

printing. A database of microarchitectures is assembled prior to large-scale optimization, consisting of different families derived through topology optimization targeting defined material properties. A tiling algorithm then maps the discrete microarchitectures to the large-scale domain in order to fulfil functional objectives. Liu et al. [START_REF] Liu | Rapid modeling and design optimization of multi-topology lattice structure based on unit-cell library[END_REF] considered the connectivity between lattice cells, and divided the lattice cells of constructed unit-cell library into seven series to ensure that lattice cells of each series could be connected. Panetta et al. [START_REF] Panetta | Elastic textures for additive fabrication[END_REF] explored a wide space of truss-like, symmetric 3D patterns to obtain a small family. This pattern family can be printed without internal support structure on a single-material 3D printer and can be used to fabricate objects with prescribed mechanical behavior. As shown in Fig. 19, the microstructures from same family are mapped to each voxel according to the property distribution of the designed objects.

Two-scale concurrent optimization

However, above works are restricted to a single scale (either macro or microscale), while the concurrent optimization at both scales is limited. The concurrent topology optimization refers to simultaneously devise the best structural topology with the most compatible material microstructure, with the aim of increasing the design freedom to find the optimal performance in the macrostructure [START_REF] Li | A new multiscale topology optimization method for multiphase composite structures of frequency response with level sets[END_REF]. The development of the AM can provide a feasible solution to fabricate the concurrent multiscale designs with either single material [START_REF] Wang | Multiscale reliability-based topology optimization methodology for truss-like microstructures with unknown-but-bounded uncertainties[END_REF] or multiphase materials [START_REF] Vogiatzis | Topology optimization of multi-material negative Poisson's ratio metamaterials using a reconciled level set method[END_REF][START_REF] Vicente | Concurrent topology optimization for minimizing frequency responses of twolevel hierarchical structures[END_REF]. In the existing concurrent topology optimization framework, macro and micro densities are introduced as independent design variables for macro structure and material microstructure, respectively. Optimizations at two scales are integrated into one system with the homogenization theory. Penalization approaches are adopted in both scales to ensure clear topologies, i.e., SIMP in the micro scale and PAMP (Porous Anisotropic Material Penalization) in the macro scale. Two classes of design variables (density fields) are independently defined [START_REF] Deng | Concurrent topology optimization of multiscale structures with multiple porous materials under random field loading uncertainty[END_REF].

As illustrated in Fig. 20, the two-scale concurrent topology optimization framework, where the material microstructure is assumed to be the same throughout the macroscopic structure, is the focus of this work. The macrostructure is assumed to be composed of uniform materials, whose material microstructures could be regarded as composed of repetitive unit cells in which a voxel-based representation specifies whether a finite element contains base material (voxel=1) or not (voxel=0). and Breitkopf [START_REF] Xia | Recent Advances on Topology Optimization of Multiscale Nonlinear Structures[END_REF], where several different cellular materials are designed for a layered structure following a two-step design procedure. Another concurrent topology optimization model with multiple porous materials is proposed by Deng and Chen [START_REF] Deng | Concurrent topology optimization of multiscale structures with multiple porous materials under random field loading uncertainty[END_REF].

By combining the concurrent topology optimization framework and the discrete material optimization interpolation model, the distribution of subdomains on the macroscale and the topology of the material microstructure in each subdomain can be automatically determined. L-Beam optimized design by a generic approach [START_REF] Sivapuram | Simultaneous material and structural optimization by multiscale topology optimization[END_REF] allowing any number of unique microstructures indicates that increasing the number of allowable microstructures leads to a solution with an improved objective. Li et al. [START_REF] Li | Design of Architected Materials for Thermoelastic Macrostructures Using Level Set Method[END_REF] introduced a larger number of material microstructures in architected thermoelastic materials optimization, the overall performance was improved due to the expanded design space. However, when multiple architected materials with spatial variations in a structure are considered, a challenge arises in topological solutions, which may not be connected between adjacent material architecture. Du et al. [START_REF] Du | Connecting microstructures for multiscale topology optimization with connectivity index constraints[END_REF] introduced connectivity index (CI) to quantify the topological connectivity, and added it as constraint in multiscale topology optimization to achieve connected architected materials.

However, the unit cell of microstructure is composed of single base material with single base material property. In order to further explore design freedom provided by multi-material AM (e.g., Polyjet), Zhu et al. [START_REF] Zhu | Two-scale topology optimization with microstructures[END_REF] proposed a two-scale topology optimization with microstructures consisting of two base materials (soft material and rigid material). As shown in Fig. 21, a continuous gamut of material properties is derived from a precomputed database of discrete microarchitectures prior to large-scale optimization. Topology optimization is then used to optimize the material properties of the large-scale domain within the gamut to satisfy functional objectives. The continuous solution from the topology optimization is then mapped to discrete microarchitectures from within the database. This effective approach requires a large computational resource to assemble the databases prior to large scale optimization.

Multi-material optimization methods of heterogeneous lattice structures

From the multi-material point of view, the material space of the multi-material objects is no longer confined to a single material, it has been extended to multi-material space. Therefore, compared with the topology optimization of a single material object, multi-material topology optimization (MMTO) refers to finding the optimal layout form of variety of materials in the design space under a given design domain, a given constraint and a load. Then, the best way to pass force is obtained to achieve a predetermined function or movement, and to achieve a certain design goal. For MMTO, the problem could be generalized to

( ) 1,2, , Min , 0 Subject
to k M D i j i i k g u Ku f ϕ = ⊂ Ω ≠ Ω ∩ Ω = ∅ Ω ≤ =  ( 4
)
where ϕ is the objective, such as compliance, stress, or volume; k Ω is the topology to be computed for material " k " to be computed; M is the number of materials; D is the domain within which the topology must lie; u is the finite-element displacement vector K is the finite-element stiffness matrix; f is the external force vector, and i g is the constraints on volume, stress, buckling, etc.

The MMTO problem assumes that every point within the design space has a distinct material associated with it (or is void). This differs from functionally graded material optimization, where a mixture of base materials is allowed [START_REF] Mirzendehdel | A pareto-optimal approach to multimaterial topology optimization[END_REF]. It mainly includes Solid Isotropic Microstructure with Penalty (SIMP), Level Set Method (LSM) and Discrete Element Method (DEM).

SIMP for MMTO

SIMP method needs different material interpolation schemes to meet the requirements of different base materials and material properties. Thomsen et al. [START_REF] Thomsen | Topology optimization of structures composed of one or two materials[END_REF] firstly extended multiple materials to SIMP in 1992. Sigmund proposed [START_REF] Sigmund | Design of multiphysics actuators using topology optimization-Part II: Two-material structures[END_REF] a twomaterial interpolation scheme for designing electrically and thermally driven microactuators. Kruijf et al. [START_REF] De Kruijf | Topological design of structures and composite materials with multiobjectives[END_REF] find the optimal structures with the maximum stiffness and minimum resistance to heat dissipation and tailor composite materials with effective thermal conductivity and bulk modulus. Two-phase ill-ordered base materials (i.e. one has a higher Young's modulus, but lower thermal conductivity while another has a lower Young's modulus but higher conductivity) are assumed in order to observe competition in the phase distribution defined by stiffness and conduction. Hvejsel et al. [START_REF] Hvejsel | Material interpolation schemes for unified topology and multi-material optimization[END_REF] proposed two multi-material interpolation schemes as direct generalizations of the well-known SIMP and RAMP material interpolation schemes originally developed for isotropic mixtures of two isotropic material phases. Jeong et al. [START_REF] Jeong | Separable stress interpolation scheme for stress-based topology optimization with multiple homogenous materials[END_REF] presented a new regional constraint method based on the sorting algorithm and the applicability and limitations of the newly developed framework were discussed in the context of its application to several stress-based topology optimizations with multiple materials.

Gaynor et al. [START_REF] Gaynor | Multiple-material topology optimization of compliant mechanisms created via PolyJet three-dimensional printing[END_REF] manufactured compliant mechanism designs based on three-phase (2 solid phases plus void) topology optimization using PolyJet AM technology, which can print bulk materials covering a wide range of elastic moduli (Fig. 22). The advantages of SIMP are as followings: the design variables (cell densities) in the SIMP model are directly linked to the optimization problem, that is, the display of the topology is dependent on the design variables. The model optimization algorithm converges well and can directly calculate the discrete design sensitivity based on finite element method. Moreover, it is suitable for any shape of the design domain and is easier to integrate with existing commercial finite element software. However, SIMP method extended to m-phase materials requires (m-1) design variables for each finite element, which results in complex interpolation express and large computational cost [START_REF] Zuo | Multi-material topology optimization using ordered SIMP interpolation[END_REF]. In order to enforce the selection of, at most or exactly, one material at each design subdomain, a large number of sparse linear constraints are needed.

LSM for MMTO

The level-set function is originally developed by Osher and Sethian [START_REF] Osher | Fronts propagating with curvaturedependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] with the fundamental goal of tracking the motion of curves and surfaces. In this work, the structural boundary is implicitly represented by the zero level sets. For the multi-material structures, multiple level set functions ( )

, 1, 2 k k M φ = 
are employed to denote different phases. These level set functions are utilized to define the following subdomains [START_REF] Cui | A level-set based multi-material topology optimization method using a reaction diffusion equation[END_REF]: Wang et al. [START_REF] Wang | Color" level sets: a multi-phase method for structural topology optimization with multiple materials[END_REF] also introduced a novel level-set approach referred to as a "color" level-set representation which has its unique benefits: it is flexible to handle complex topologies. The multi-level set functions ( )

( ) ( ) ( ) Φ x 0 Φ x 0 , Φ x 0 1, 2 k k k k k k k k M x x x D  > ∀ ∈ Ω Γ  = ∀ ∈ = Γ   < ∀ ∈ Ω     (5)
, 1, 2 k k M φ =  are used to represent 2 M n =
distinct material phases and substantially reduces the number of model functions. Liu et al. [START_REF] Liu | Level set-based heterogeneous object modeling and optimization[END_REF] employed multiple level set functions to build the geometry, utilizes zero-value level set contours as material source profiles, and realized FGM blending with a signed distance-based blending function. More importantly, this model supports the concurrent structure and material optimization because of the unified level set framework for both structure and material composition representation. The initial problem setups are demonstrated in Fig. 24a and Fig. 24c. The heterogeneous source profile in Fig. 24a employs the Young's modulus combination of 1.3 and 0.33, while the heterogeneous source profile in Fig. 24c employs the Young's modulus combination of 1.3 and 0.98. The optimization results are demonstrated in Fig. 24b and Fig. 24d, respectively.

DEM for MMTO

Unlike continuous topology optimization method, i.e. SIMP. Discrete elements such as 1D truss element (ground structure method) [START_REF] Zhang | Material nonlinear topology optimization using the ground structure method with a discrete filtering scheme[END_REF] or 3D bar element (geometry projection method) [START_REF] Norato | A geometry projection method for continuum-based topology optimization with discrete elements[END_REF] can be directly applied to MMTO for MSHLS. Recently, several papers about this topic have been published. Zhang et al. [START_REF] Zhang | Multi-material topology optimization with multiple volume constraints: a general approach applied to ground structures with material nonlinearity[END_REF] proposed an efficient MMTO formulation considering material nonlinearity. The proposed formulation handles an arbitrary number of candidate materials with flexible material properties, features freely specified material layers, and includes a generalized volume constraint setting. An application to MMTO of truss networks considering multiple load cases and nonlinear constitutive behavior is shown in Fig. 25. Kazemi et al. [START_REF] Kazemi | Topology optimization of structures made of discrete geometric components with different materials[END_REF] presented a geometry projection method for the simultaneous topology optimization and material selection of structures made by the union of discrete geometric components, where each component is made of one of multiple available materials.

Michell cantilever example shows that the proposed formulation can be readily extended to any number of materials (Fig. 26). The advantages of DEM are as followings: Optimal design of truss like structures with discrete design variables share many mathematical properties with multi-material topology optimization. For ground structure method, the final optimal truss structure is realized by selecting an optimal substructure from a pre-defined ground structure.

Typically strut cross-section dimensions (such as diameter) are the design variables.

Local strain or stress values from the FEA results are used to compute the strut sizes. If a strut size approaches the specified threshold, it can be removed, effectively changing the topology of the design. This conveniently changes a topology optimization problem into a sizing problem. Unlike ground structure methods, the discrete elements in geometry projection method need not be connected during the optimization, hence giving more freedom to the optimizer to optimally place and size the bars. By having an independent explicit geometry, this method is capable of producing designs with members of predefined cross-sections. Moreover, it can accommodate the case where the thicknesses of the bars are added at their intersections (overlapping bars), and it can compute the volume fraction and its sensitivities directly from the bar's geometry and not from the composite density. However, the discrete nature of the problem makes it difficult to be solved. There exist a large number of both mathematically well-founded as well as heuristic approaches like genetic algorithms, swarms and differential evolution techniques that are not viable alternatives for the vast majority of topology optimization problems. However, none of them is well-suited for the problems with thousands and up to millions of variables that are encountered in topology optimization approaches.

Based on the optimization design methods introduced above, a comparison between different optimization approaches is made and summarized in Table 4. Based on this comparison, it helps designers select a suitable optimization method for different cases.

Discussion

A general design process of MSHLS including input of design and functional requirements, initial design, optimized design and geometric&material composition models is proposed in Fig. 27. This process is based on the design and fabrication methods for AM proposed by Rosen [START_REF] Rosen | Computer-Aided Design for Additive Manufacturing of Cellular Structures[END_REF], Tang [START_REF] Tang | A survey of the design methods for additive manufacturing to improve functional performance[END_REF] and Yang [149]. In this paper, we only focus on the design aspect. With the input of design and functional requirements, initial design including either Functional Volumes (FVs) and Functional Surfaces (FSs) with solid&skin or with lattice can be obtained. Then the initial design result is the input to the simulation step. The simulation result then goes back to modeling stage for the optimization of the initial design to get optimized design. During this design process, optimization methods (multi-scale optimization and multi-material optimization), simulation support (mixtures rule and numerical analysis), and modeling methods (material composition modeling and multi-scale geometric modeling) as main design methods support optimized design phase to obtain final geometric & material composition models of MSHLS.

Optimized Design

Initial Design

Design and Functional Requirements

Geometric&Material Composition Models

Multi-scale Geometric Modeling

Material Composition Modeling

Modeling Methods

Simulation Support

Multi-material Optimization

Multi-scale Optimization

Optimization Methods

Mixtures Rule Numerical Analysis

Fig. 27 The design process of multi-scale heterogeneous lattice structure.

To facilitate the use of the above-mentioned modeling, simulation and optimization methods for MSHLS design, it is necessary to describe the scope of application of the method. According to the minimum basic elements describing the is not only applicable to Boolean operation but also to MMTO (Section3.3.2).

Summarising the overview of modelling methods, there is a gap between geometrical complexity of heterogeneous lattices and inability of current geometric modelling methods to handle such complexity. Moreover, volumetric modelling allows direct transition to structural simulations but is even more complicated. Thus, there is a need in a novel volumetric geometric modelling approach to support MSHLS. Such approach should be extensively tested on numerous use-case of lattice structures and other complex geometries, measuring its performance and speed comparing to existing approaches. Moreover, the concept of level of detail should be adapted for lattice structures, as they do not have CAD-features that are commonly associated with their level of detail in conventional modelling.

Simulation methods such as homogenization can be applied based on the calculated effective properties. However, how to accurately predict the effective properties for MSHLS is a critical issue, since micro cell of this type of structure may dramatically change in a certain local region. RVE (Representative Volume Element) for this type of structure is difficult defined. There is no cell which can strictly satisfy the periodic boundary condition. By directly choosing the smallest unit cell as RVE will cause the significant discrepancy between simulation result and physical testing result.

Another method of mixtures rule is the simplest way. However, the accuracy of this method is limited, since it fails to consider the geometric configuration of the microstructures of the material composition. In most cases, it can only provide a rough estimation of upper or lower bound of effective properties. FEA methods could be used based on decomposition-based models (Section2.1.1). However, the computational cost for a fine-meshed model is too heavy to fulfillment.

For the optimization methods, to efficiently solve the proposed optimization formulation is another critical issue needs to be solved. Comparing to the existing optimization methods for homogenous objects, the proposed optimization formulation contains more design parameters. Moreover, those design parameters may be coupled with each other. For example, during the optimization the stiffness of the heterogeneous structures. The material distribution on macroscale is coupled with the microstructure of materials in different regions. On the one hand, the optimized microstructure can be generated based on predefined material compositions. On the other hand, the effective properties of optimized structures will affect the distribution of material compositions in turn. To solve this complex optimization problem, there are several potential methods which can be adopted in this problem. Given that MSHLS needs to deal with a large number of design parameters, one strategy for MMTO is to introduce multi-resolution method to obtain multiple discretization level. The discretization level for the displacement elements, the density elements, and the design variables can be chosen independently, thus it can realize high-resolution MMTO with relatively low computational cost by employing a coarser discretization for finite elements and a finer discretization for both density elements and design variables. Another strategy is to reduce the number of design parameters by using the heuristic algorithms. Generally, the heuristic algorithm may be applied based on initial analysis of structure to divide the initial design domain on a single design scale into several sub-regions. The material or its microstructure are supposed to be homogeneous in each sub region during the design optimization process. This method can significantly decrease the number of design parameters and computational time for the defined optimization problem. To optimize the structures on different scales sequentially. The method (Section3.2.3) tries to decouple the relation between design parameters on different design scales. The effects of design parameters defined on the lower scale will be pre-estimated during the optimization iteration for the design parameters defined on the upper scale. Additional optimization constraint will be added to limit the discrepancy between the estimated effects and the real effects of the design parameters defined in the lower scales. This strategy has been widely used in the Multidiscinplnary Design Optimization (MDO), and known as IDF (Indiviual Discinplinary Feasible) analysis. It can be used for design and optimization a complex coupled system.

The developed geometric&material composition model for MSHLS is a universal model. It cannot only be used to represent the information of designed and optimized structures for industry and engineering applications. This model can be also used to represent the information of bio-tissues which also has hierarchical complexities such as human's bone. Based on the developed model, researchers from bioengineering and biomedical fields can further investigate the mechanical behaviour of these tissues.

Conclusion and perspectives

The contribution of this review is to provide a comprehensive view of the MSHLS design process, and of the modeling, simulation and optimization strategies currently 2015, "Additive manufacturing of metal cellular structures: design and fabrication," Jom, 67 [START_REF] Bourell | The roadmap for additive manufacturing and its impact[END_REF], pp. 608-615. Permission to reprint from IEEE copyright 2018 [START_REF] Leung | Approximate functionally graded materials for multi-material additive manufacturing[END_REF]. Permission to reprint from ACM copyright 2015 [START_REF] Schumacher | Microstructures to control elasticity in 3D printing[END_REF]. Permission to reprint from Elsevier copyright 2019 [START_REF] Deng | Concurrent topology optimization of multiscale structures with multiple porous materials under random field loading uncertainty[END_REF].

Figure Captions List

Fig. 1

 1 Fig. 1 Heterogeneous lattice structure: a) Heterogeneous lattice with gradient strut dimensions, b) Heterogeneous lattice with variational material compositions.

Fig. 2

 2 Fig. 2 Multiscale structures with hierarchical lattice architectures. Permission to reprint from Springer Nature copyright 2016 [8].

  Many scholars have carried out an in-depth research on the techniques of material composition modeling. It mainly includes decomposition-based models, set-based models, explicit function-based models, control feature-based models, control pointbased models. A brief review of the methods is provided here to clarify the difference and limitations of present works:

Fig. 4

 4 Fig. 4 2D and 3D heterogeneous lattice structures. Permission to reprint from Elsevier copyright 2014 and ASME copyright 2019 [22, 23].

Fig. 5 2 . 1 . 3

 5213 Fig. 5 Material blending with three material source profiles. Permission to reprint from Elsevier copyright 2019 [30].

Fig. 6

 6 Fig. 6 Common material composition variation features in MSHLS. Permission to reprint from Elsevier copyright 2019 [31]. As shown in Fig. 7, distance-based material distribution function, e.g., Piecewise, Polynomial, Power, Gaussian, Cauchy, etc, can realize different material gradients. The material distribution is invariant under some transformations, such as rotation or translation. Due to the effective evaluation of point-to-feature distances, the representation supports efficient query of material composition. Moreover, the user can select the reference features and material distribution function through the predefined function library and can intuitively and flexibly complete the interaction modeling of multi-material objects.

Fig. 7

 7 Fig. 7 Effect of different material functions on material-distribution. Permission to reprint from IEEE copyright 2018 [33].

Fig. 8 1D

 8 Fig. 8 1D heterogeneous features: a) a heterogeneous line, b) a heterogeneous B-spline. Permission to reprint from Elsevier copyright 2005 [37].

Fig. 9

 9 Fig. 9 An example of LODs associated with CAD features. Permission to reprint from John Wiley

Fig. 10

 10 Fig. 10 Lattice structure voxelization, adapted from [74].

Fig. 11 A

 11 Fig. 11 A standard cubic unit element: a) and a cylindrical shape with its cylindrical sector unit element, b) used for stress analysis. Permission to reprint from Elsevier copyright 2013 [75].

Fig. 12 FEA

 12 Fig. 12 FEA-based simulation aided design: a) multi-material fiber network of artificial disc. Permission to reprint from ASME copyright 2019 [98], b) medical tweezers with spatial material distribution. Permission to reprint from ACM copyright 2015 [100].

Fig. 13

 13 Fig. 13 Overview of the digital design and manufacturing framework for soft lattice structures. Permission to reprint from Elsevier copyright 2019 [106].

Fig. 14

 14 Fig. 14 Material microstructures design: a) level-set-based design. Permission to reprint from Elsevier copyright 2016 [108], b) design of architected truss lattice materials with three base materials. Permission to reprint from Elsevier copyright 2020 [109], and c) voxel-wise material distribution. Permission to reprint from ACM copyright 2015 [110].

Fig. 15

 15 Fig. 15 An example of F-P-P-D model. Permission to reprint from Springer Nature copyright 2018

Fig. 16

 16 Fig. 16 Relative density mapping to strut diameter. Permission to reprint from Elsevier copyright 2015 [116].

Fig. 17

 17 Fig. 17 Example of mechanical properties mapping process and its manufacturability verification by 3D printing: a) material properties diagram, b) 3D mapping result, c) 3D printing result. Permission to reprint from Elsevier copyright 2018 [117].

Fig. 18

 18 Fig. 18 Structures used for FEA: a) Solid (SIMP solution), b) Intersected Lattice of D-P, c) Intersected Lattice of BCC, d) Graded Lattice of D-P, e) Scaled Lattice of D-P, and f) Uniform Lattice of D-P [118].

Fig. 19

 19 Fig. 19 Deformation of an object with varying material properties per voxel, and the same object with the material in each voxel replaced with the corresponding pattern. Permission to reprint from ACM copyright 2015 [124].

Fig. 20

 20 Fig. 20 Two-scale hierarchical structure with uniform microstructure. Permission to reprint from Elsevier copyright 2019 [130]. The most commonly applied strategy is designing a uniform material microstructure at the microscopic scale either for a fixed or concurrently changed structure at the macroscopic scale. Obviously, such designs have not yet released the full potentiality of concurrent two-scale designs. A further step has been made by Xia

Fig. 21

 21 Fig. 21 Multi-scale topology optimization framework with microstructure database. Permission to reprint from ACM copyright 2018 [11].

Fig. 22 3

 3 Fig. 22 3-phase inverter result by multiphase SIMP approach. Permission to reprint from ASME copyright 2014 [140].

  where D represents the design domain including all admissible shapes. k Ω denotes the kth material region with positive value of the kth level set function, k Γ is the boundary of the kth material, negative value of the kth level set function signifies the domain not containing the kth material. x represents a point located in D . M is the number of the level set function. An example for the design domain containing three level set functions is illustrated in Fig. 23.

Fig. 23 A

 23 Fig. 23 A schematic diagram of the design domain which includes three level set functions. Permission to reprint from Elsevier copyright 2016 [143].

Fig. 24

 24 Fig. 24 Optimization results of the Michell structure problem with heterogeneous source profiles. Permission to reprint from Elsevier copyright 2019 [30]. The advantages of level set are as followings: level set-based model is one of a few exceptions that could support both aspects of modeling and optimization. The gray area problem is solved by the unique advantages of smooth and clear boundaries of various materials and convenient extraction of topological configurations, and the optimization results can be directly used for However, since the design variables are indirectly linked to optimization problems, it involves the approximation of finite elements that are cut by the horizontal set, thus affecting the optimization accuracy. Level-set-based topology optimization methods suffer from local minima and dependency of the final result on the initial guess. This property is a consequence of the shape optimization character of LSM. Seeding an initial design with a large number of holes or inclusions may lead to numerical issues and does not necessarily mitigate the dependency of the final result on the initial guess.

Fig. 25 3DFig. 26

 2526 Fig. 25 3D crane design subjected to multiple load cases: a) optimized structure for the 3D crane design with three materials, b) printed model using FDM process. Permission to reprint from Springer Nature copyright 2017 [147].

  material information of MSHLS, voxel-based, volumetric mesh-based, truss or beambased and curve feature-based design methods are summarized. Among them, voxel-based modeling, simulation and optimization methods include voxel model (Section2.1.1), physics-based method (Section3.1.1), homogenization method (Section3.1.2), and SIMP method (Section3.3.1); Volumetric mesh-based methods include volume mesh model (Section2.1.1) and FEA-based method (Section3.1.3); Truss or beam-based methods include control feature-based model (Section2.1.3), FEA-based method (Section3.1.3) and DEM (Section3.3.3); Curve feature-based methods include control point-based model (Section2.1.4), IGA-based method (Section3.1.4). Especially, as an implicit function, level set-based method (Section2.1.2)

 6 Acknowledgment

 6 available, which are still lacking. In the end, based on the summaries of different modeling, simulation and optimization methods, it comes out that further researches are needed and several future work perspectives concerning the design of MSHLS have been pointed out:  To develop a multi-scale model which can represent the information of MSHLS. This model should satisfy the following requirements. Firstly, it should be able to store both material and geometry information on multiple design scales. The relations between different design scales are also needed to be clearly defined in this model. Secondly, the developed model should have a compact and efficient data structure, which enable the quick accessing and visualization of data on different design scales. Thirdly, the developed model should be general enough to incorporate different types of MSHLS. Fourthly, this model should also be easily converted manufacturing readable data model, which enable the direct fabrication via selected AM process. To develop a multiscale simulation model of MSHLS. In this model the effective properties of material defined in the bottom scale need to be calculated based on the numerical homogenization methods. Different numerical homogenization methods will be developed or modified to evaluate the different type of microstructures. For example, if the microstructures are lattice structure, numerical homogenization methods can be applied. By carefully selecting the suitable numerical homogenization methods, it can significantly shorten the computational time. Moreover, in order to improve the accuracy of the calculation, the modification of calculated effective properties from numerical homogenization methods also needs to be done. The detailed modification algorithm needs to be further investigated. Based on evaluated effective properties, the numerical simulation model such as FEA model can be applied to evaluate the performance of modeled multiscale MSHLS. To investigate general optimization parameterization formulation. On different design scale, the parameters which can efficiently control the distribution of material compositions and its microstructures should be figured out. For example, the coefficient of exponential function can be used to control the material distribution in design domain, while the implicit level set functions can be used to parameterize the microstructures of materials. The selected design parameters can be directly converted to the generalized units defined on each level. This work is supported by financial support from the Program of the Academic Excellence Foundation of BUAA for PhD Students, China Scholarships Council (No. 201806020095), and National Sciences and Engineering Research Council of Canada (Discovery Grant RGPIN-2018-05971).

Fig. 1

 1 Fig.1Heterogeneous lattice structure: a) Heterogeneous lattice with gradient

Fig. 3

 3 Fig. 3Multi-material-property configuration[START_REF] Liu | Material-Unit Network for Multi-Material-Property and Multiscale Components[END_REF]

Fig. 5

 5 Fig. 5Material blending with three material source profiles. Permission to

Fig. 6

 6 Fig. 6Common material composition variation features in MSHLS. Permission

Fig. 7

 7 Fig. 7Effect of different material functions on material-distribution.

Fig. 8 1D

 8 Fig. 8 1D heterogeneous features: a) a heterogeneous line, b) a heterogeneous

Fig. 9 Fig. 10

 910 Fig. 9An example of LODs associated with CAD features. Permission to reprint

Fig. 11 A

 11 Fig. 11 A standard cubic unit element: a) and a cylindrical shape with its

Fig. 12 FEA

 12 Fig. 12 FEA-based simulation aided design: a) multi-material fiber network of

Fig. 13

 13 Fig. 13 Overview of the digital design and manufacturing framework for soft

Fig. 14

 14 Fig. 14 Material microstructures design: a) level-set-based design. Permission to

Fig. 15

 15 Fig. 15 An example of F-P-P-D model. Permission to reprint from Springer

Fig. 16 Fig. 17

 1617 Fig. 16 Relative density mapping to strut diameter. Permission to reprint from

Fig. 18

 18 Fig. 18 Structures used for FEA: a) Solid (SIMP solution), b) Intersected Lattice of

Fig. 19

 19 Fig. 19 Deformation of an object with varying material properties per voxel, and

Fig. 20

 20 Fig. 20 Two-scale hierarchical structure with uniform microstructure.

Fig. 21

 21 Fig. 21 Multi-scale TO framework with microstructure database. Permission to

Fig. 22 3 Fig. 23 AFig. 24

 32324 Fig. 22 3-phase inverter result by multiphase SIMP approach. Permission to

Fig. 25 3D

 25 Fig. 25 3D crane design subjected to multiple load cases: a) optimized structure

Fig. 26

 26 Fig. 26 Michell cantilever optimal designs. Red (3) indicates material

Fig. 27

 27 Fig.[START_REF] Gupta | Heterogeneous object modeling with material convolution surfaces[END_REF] The design process of multi-scale heterogeneous lattice structure.

  Software developers could get insights for the development

	or adjustment of design tools to represent the information of designed and optimized
	structures for industry and engineering applications. Researchers could get a
	comprehensive view on the topic together with indications for developing more
	integrated and evolutionary design methods.
	After the introduction above, this paper is organized as follows: firstly, material
	composition modeling methods and multi-scale geometric for representation of material

and geometry information are separately discussed in Sec.2. Moreover, optimization methods including multi-scale optimization and multi-material optimization design methods suitable for MSHLS are respectively reviewed in Sec.3. In addition, a design process of MSHLS is discussed in Sec.4. Finally, this paper is wrapped up with a conclusion and some prospects for the future research in Sec.5.

2 Modeling methods of multi-scale heterogeneous lattice structure

	Generally speaking, if a heterogeneous object is composed of 2 n different
	materials and attached with 3 n different physical properties, the modelling space of
	heterogeneous object can be denoted as the tensor product of 1 n dimensional
	geometry space 1 n E , 2 n dimensional material space	M	2 n	and 3 n dimensional
	property space	3 n

Table 1 Comparison of existing models for heterogeneous materials Methods Including Advantages Disadvantages

 1

	Decomposition-	Voxel-based	Voxel model is	Inexact and
	based models	model	suitable for	computationally expensive
		volume mesh-	visualization	Non-trivialness to
		based model	Volume mesh model	manually
			is easy for FEA	manipulate the material
				distribution
	Set-based	R m -set	R m -set model is easy	Incompatible with other
	models		to integrate	modeling formats and
			constructive	cannot store topology
		F-Rep model	operation	information
			F-Rep model
			facilitates topology
			optimization
	Control feature-	Single and	More favorable in	Heavily based on the
	based models	multiple	capturing the	features with predefined
		feature-based	designers' intents	material distributions
		model	Suitable for FGM
			representation
	Control point-	B-spline and	Compact in both	Heavily based on spatial
	based models	radial basis	geometry and	parameterizations
		function	material
			representations

Table 2 Comparison between polygon mesh modeling and voxel-based modeling approaches for MSHLS

 2

	Methods	Advantages	Disadvantages
	Polygon mesh	Computational efficiency	Not suitable for complex geometrical
	modeling		shapes
		LOD concept is applicable
	Voxel-based	Handles high complexity	Not computationally efficient
	modeling
		No need for resolution higher	Introduce distortion when modeling on
		than required for printing	low resolution
		Efficient Boolean operations
		Feature recognition is well
		applicable

Table 3 Comparison of existing simulation methods for MSHLS

 3

	Methods	Advantages	Disadvantages
	Physics-based method Realistic visual effects in real time	Low accuracy
			Able to simulate multiple interspersed
			materials of varying properties
			Computationally efficient
	Homogenization	Computationally efficient	Satisfying local periodicity
	method		Easy to conduct multiscale simulation	hypothesis owing to scale
			through repeated unit cell behavior	effect
	FEA-	3D	Suitable for arbitrary geometry	High computation cost and
	based	element	High accuracy and heterogeneity	time consuming
	method	1D	Suitable for truss-like structure	Unable	to	simulate
		element	Low heterogeneity	multiple	interspersed
			Relatively low computation cost	materials at joints
				Low heterogeneity
	IGA-based method	Suitable for curves frame
			The parameters control points of the
			model can be directly optimized, and
			the boundary of the optimization results

Table 4 Comparison of existing optimization methods for MSHLS

 4

		Methods	Advantages	Disadvantages
		SIMP-based method	The algorithm has good	Complex
			convergence, is suitable	interpolation
			for any design domain, and	expresses and large
			is easy to integrate with	computational cost
			commercial FEA software
	Multi-material	LSM-based method	The boundary of	Final results
	optimization		optimization results is	strongly depend on
	methods		clear, which can be	starting guess
			directly used in
			manufacturing
		DEM-based method	It is suitable for truss-like
			MSHLS, and can deal with

Table Caption List

 Caption Table 1 Comparison of existing models for heterogeneous materials Table 2 Comparison between polygon mesh modeling and voxel-based modeling approaches for MSHLS Table 3 Comparison of existing simulation methods for MSHLS Table 4 Comparison of existing optimization methods for MSHLS