A novel bio-inspired design method for porous structures: Variable-periodic Voronoi tessellation

This paper introduces a novel approach, namely Variable-Periodic Voronoi Tessellation (VPVT), for the bio-inspired design of porous structures. The method utilizes distributed points deﬁned by a variable-periodic function to generate Voronoi tessellation patterns, aligning with a wide diversity of artiﬁcial or natural cellular structures. In this VPVT design method, the truss-based architecture can be fully characterized by design variables, such as frequency factors, thickness factors. This approach enables the optimal design of porous structures for both mechanical performance and functionality. The varied, anisotropic cell shapes and sizes of VPVT porous structures provide signiﬁcantly greater design ﬂexibility compared to typical isotropic porous structures. In addition, the VPVT method not only can design micro-macro multiscale materials, but is also applicable for the design of meso-macro scale truss-based porous structures, such as architecture constructions, biomedical implants, and aircraft frameworks. This work employs a Surrogate-assisted Diﬀerential Evolution (SaDE) method to perform the optimization process. Numerical examples and experiments validate that the proposed design achieves about 51.1% and 47.8% improvement in compliance performance and damage strength, respectively, than existing studies.


Introduction
Porous or lattice structures, also known as cellular structures, are ubiquitous and widely used in nature and many engineering fields, due to their high specific strength and multi-functionality [1][2][3].Some of the examples of porous structures in nature and engineering systems are presented in Fig. 1: (a) A two-phase soil porosity model in nature and its pores and hard particles in phases increase with the soil depth, which is analogous to functionally graded materials with unidirectional variation in pore density; (b) Artificially manufactured porous aluminum foam material, with extremely high specific strength and load-carrying capacity; (c) Human bone trabecular possesses a multi-layer porous structure with small bone branches, resulting in high strength and impact resistance; (d) A gymnasium building constructed by variableperiodic porous structures exhibits superior deformation and buckling resistance; (e) Dragonfly wings naturally exhibit variable-periodic porous structural frameworks, which provide extensive benefits for airflow controlling and stability maintenance.Regarding functionalities at the engineering level, porous structures are known to have many excellent properties such as lightweight, vibration absorption, wear re-sistance, electrical shielding and multiple load-carrying capacity [4][5][6][7].In addition, the material-transferring property of porous structures is also particularly important at the level of biological systems, such as trabecular bone, where pores provide space and passages for vessels or plant stems to transport nutrients [8].To date, a variety of ultralight bio-inspired porous materials and structures have been proposed, designed and fabricated, which have attracted great interests from both academic and industrial communities [9][10][11].To achieve optimal mechanical performance and behaviour for porous structures, researchers have attempted to apply Topology Optimization (TO) methods to design the layouts and configurations of porous structures [12,13].The pioneering work on this topic can be traced back to Bendsøe, who described the design domain in terms of infinitely small square cells with rectangular holes and developed a numerical homogenization method to assign equivalent mechanical properties to these cells [14].This groundbreaking work paved the way for subsequent research using the homogenization method based TO to design porous structures.Up to now, many different TO methods have been developed and successfully applied to design porous structures for various engineering applihttps://doi.org/10.1016/j.matdes.2024  cations in mechanical, aerospace, biomedical and architectural fields, etc. [15][16][17].
The advent of 3D printing technology achieves substantially large design freedom for porous structures, in which the number of layers and the complexity of configurations increase significantly to an unprecedented extent [18].Notably, the discussed 'porous structure' in this work is a general concept that contains an inner complex structure rather than a pure porous structure.Although the traditional monolithic TO method is theoretically capable of designing porous structures with diverse geometric patterns across multiple scales, its practical application is hindered by high computational costs [19].To this end, there have been numerous efforts to develop advanced multiscale optimization methods for the design of 3D-printed porous structures [15,19,20].In their methods, the original problem is carried out with respect to micro and macro scales, respectively, using a two-step optimization procedure.The macro-scale optimization process determines the optimal material distribution in the design domain with a set of discretized macroscopic elements.While, the micro-scale optimization is applied to obtain specific porous or lattice architectures whose effective material property matches with that of macroscopic element.Finally, by tessellating the determined micro-structures into each macro-scale element, an integral design for a multi-layer porous structure with optimal performance can be achieved.The homogenization of the micro-structure plays a key role in the multiscale optimization process, which directly couples the computation between macro-and micro-scales.Up to now, a large number of elemental homogenization and de-homogenization techniques have been developed, which can be based on relative density mechanism [21], bulking modulus [22], effective elastic tensor mechanism [23,24], or strain-driven [25].Besides the above studies derived by an analytical computation model or experiment data, many studies have recently begun to utilize the data-driven or machine learning algorithms [26][27][28] to complete the homogenization processing.With this approach, the implicit, complicated relationship between microarchitecture and macroscopic elemental property can be derived and applied easily in the design.Such methods have presented a natural advantage in the multiscale design of meta-materials or irregular structures [29].However, the homogenization method often suffers from a "size effect" issue in practical applications.Most methods exhibit decreased accuracy as the length difference between two physical scales decreases, as highlighted in existing studies [30][31][32].In the design of porous structures with meso-macro scales, such as bridge supports or gliding wings, the homogenization method may inaccurately evaluate mechanical properties or the elastic tensor, leading to a failure of multiscale optimization.The conceptual framework of macro-, meso-, and micro-scales usually lacks strict numerical boundaries.In the design context of Voronoi-based porous materials, it is necessary to let the number of voids in the unit cell surpass a certain threshold to ensure the design problem at the 'micro-macro' scale.This threshold ensures alignment of the designed model with the Gibson-Ashby model [6,33,34], therefore, it guarantees that the homogenization of each unit cell can achieve a relatively accurate isotropic elasticity property.The recognition of this "size effect" underscores the significance of developing innovative design methods, which are capable of accommodating various scales and providing extensive geometry tailoring capabilities for material topology [35,36].
The early framework of the multiscale TO method is generally in a 'varied-density parameterized unit cell' mode, which adopts predetermined lattice or pore unit cells (e.g.BCC, Hexagon, Kagome types) as the micro-architecture [37].Due to the fixed structure, these lattice units only require offline homogenization prior to the optimization process.The tunable relative density of the lattice unit is used as the only kind of design variable, which generally exhibits an exponential mapping function to the effective modulus of the lattice.This feature can achieve very fine geometric details with a minimum computational cost comparable to the mono-scale topology.However, the unique fixed microstructure also greatly reduces the solution space of the microstructures to one dimension.In addition, many typical square-and X-shaped cells have been proved with a theoretical upper bound on their effective modulus according to the Hashin-Shtrikman bounds [19].Later, several methods with a 'multi-parameterized unit cell' mode have been developed.The parameterized unit cell is defined with multiple parameters e.g.[17,23], and provides a large design freedom for optimization.This kind of unit cell allows, to some extent, an adaptation of the microstructural anisotropy to the local stress direction during the optimization process.Consequently, the optimized design will also be close to the theoretical limit of effective modulus.In the last decade, "concurrent multiscale optimization" becomes much more popular, which can have a hierarchical or concurrent optimization framework at both two length levels [15,28,38].The method performs an additional local optimization to determine the microscopic spatial distribution of solid and void phases.Compared to the above two predefined microstructure modes, the micro-architecture of this mode is much more diverse, which can greatly improve both the mechanical performance and functionality of the optimal structure.In the meanwhile, ensuring good connectivity between different microstructures or unit cells has also become a critical and difficult point in this method.In addition, a large number of numerical homogenization operations are usually required in concurrent multiscale optimization, resulting in a relatively cumbersome computational cost.However, regardless of the fixed or variable microstructure lattice design, all the above optimization methods in context belong to a 'restricted unit cell' mode, where the macroscopic elements always present a fixed size and shape.As shown in Fig. 2, this feature usually leads to a relatively inflexible and artificial pattern in the final structure, whereas the elements in a natural porous structure, i.e., dragonfly wings, often exhibit a much greater diversity in terms of shapes and sizes.Over the past three years, several studies have suggested that the periodic design mode of porous or lattice structures significantly limits the design freedom and the optimum solutions.As a result, more and more researchers resort to the development of new methods with 'unrestricted unit-cell design' [9].
Evolutionary algorithms (EA) have become increasingly popular in the field of structural optimization due to the rapid development of computational power in recent years [39,40].Compared with gradientbased optimization methods, the EA-based methods avoid the use of explicit expressions of sensitivity information and apply crossover and mutation operations to update design variables.Although the convergence speed of EA-based methods is slower than that of gradient-based algorithms, it is suitable for solving optimization problems with strong non-convexity and complexity and usually can avoid trapping in local optima.In this work, a revised surrogate-assisted difference evolution  [38].The design of artificial porous materials appears monotonous due to the influence of the predefined macro-scale element size on the final design pattern.On the other hand, real biological porous tissues exhibit immense diversity, with lattices shaped like rectangles, pentagons, or even filled.method (SaDE) [41,42], as an advanced variant of EA-based algorithms, is applied to perform the optimization design of porous structures.The SaDE algorithm uses a sequential approach with global and local search phases, resulting in a balance between exploitation and exploration.In doing so, the SaDE method simplifies the process of solving objective functions through using surrogate models, and thus significantly reduces the computational cost.
This work proposed a novel variable periodic Voronoi tessellation (VPVT) method to design and optimize bio-inspired porous structures.It will demonstrate that this proposed VPVT design concept can provide a generic means to create bio-inspired porous patterns without any obvious limitations in the above context, while naturally ensuring multiscale tailoring capability and connectivity.The porous pattern mimics irregular cellular structure in nature, i.e. dragonfly or trabecular bone, achieving advanced material property and performance.In the proposed method, mesh interpolation and function-driven Voronoi seeding techniques are applied to create the grains (elements) with varying sizes and shapes for porous structures.The varied patterns of Voronoi lattices not only provide extensive design freedom for structural optimization, but also possess natural connectivity along the entire design domain.The VPVT method transforms the original topology optimization problem for porous structures into a parametric optimization problem, whose objective can be expressed in terms of various structural performances including compliance, vibration, and buckling.Compared to the conventional multiscale TO methods, since only the frequency parameters of mesh grids, without micro-scale topological elements, are needed to generate VPVT structures, the number of design variables is significantly reduced.The proposed method is applicable to both micro-macro or meso-macro multiscale optimization problems, e.g., biomedical implants, bridge brackets and architectural supports.
The paper is organized as follows.In Section 2, the VPVT design concept and modelling methodology is first introduced.Subsequently, the superior tailoring ability and optimization framework to support aimed design are presented.Section 3 presents the numerical example and experiment to validate the advantages and effectiveness of the proposed method.Finally, Sections 4 and 5 summarize the key findings and conclusions drawn from this study.

Methods and materials property for VPVT structure
To enable efficient design and optimization of porous structures, the VPVT method establishes a one-to-one mapping framework from a set of pre-defined parameters (design variables) to the porous geometry.The basic design domain is first divided into several basic unit cells, then the proposed method utilized an interpolation mesh grid method, named "frequency mesh grid", to construct the porous structures of each unit cell.The four corner nodes of each mesh grid will be assigned four pre-defined frequency values, respectively, which then will serve as the basis to generate a set of seeding points in the design domain.Afterwards, these seeding points are used to create a porous structure by introducing Voronoi tessellation at their locations.Furthermore, the number and the locations of seeding points can be controlled by manipulating the frequency values associated with each frequency mesh grid.As a result, the Variable-Periodic Voronoi Tessellation (VPVT) patterns are generated.The frequency values at each mesh grid corner are used as the design variables for the parameterized optimization procedure.

The VPVT design method
The whole design domain is decomposed into M × N elements (named as Region of Interest, ROI) by the "frequency mesh grid", as shown in Fig. 3(a).At each corner of the element in a mesh grid, the corresponding 'frequency factor' is first predefined.As such, the frequency factor  (, ) at any location (, ) in the design domain can be defined and determined through a linear interpolation procedure.

|sin(𝐶
Alternatively, Eq. ( 1) can be expressed in the following two trigonometric functions as where  0 ∈  + is the magnification factor,  (, ) denotes the frequency factor distribution in each "frequency mesh".Herein, Eq. ( 2) can be analytically calculated by simplifying it into a one-variable cubic equation with a linear interpolation frequency mesh grid.For a detailed computation of the above Equations, please refer to the Appendix A. Applying the above procedures, a variable-periodic Voronoi tessellated (VPVT) scaffold structure is generated as shown in Fig. 3(c) based on the Voronoi seeds in Fig. 3(b), however, it is inherently a onedimensional wireframe model.For practical engineering applications, this wireframe model is transformed into a 2D solid model by assigning thickness properties to each line segment (or a 3D model by introducing an additional height scale to the 2D model).A mesh grid named "thickness mesh" is defined to achieve the above process, as shown in Fig. 3(d).Similar to the "frequency mesh" method, the thickness factor at any position within the design domain is determined through an interpolation operation with respect to a "thickness mesh".For each line segment within the design domain, a thickness factor at its midpoint is calculated, as depicted in Fig. 3(e), which, subsequently, will be assigned as the thickness of the whole line segment, as shown in Fig. 3(f).Fig. 3 demonstrates the construction process of a VPVT cellular structure within a single ROI.By appropriately adjusting the design variables of 'frequency factor' and 'thickness factor', this proposed VPVT method is capable of creating diverse porous structures with different geometric characteristics and structural performances over the entire design domain, as illustrated in Fig. 4. Fig. 4(a) shows the design capability of the VPVT method to create a rhombus-based cellular structure, which is analogous to the uni-axial functionally graded material.In  In this cellular structure, the size of the porosity decreases along the radius direction from the origin, and the unit cells in some ROIs exhibit typical honeycomb hexagonal structures.Fig. 4(c) shows a typical square cellular structure, which is obtained by assigning constant frequency factor in the VPVT design method.Fig. 4(d) shows a novel diagonally concentrated porous structure, in which the materials tend to cluster in the diagonal region, resulting in dense and small cavities at the sites.The corresponding ROI image (in the lower right box of Fig. 4 shows that the design domain is enriched by long slender-shape cellular structures, in which the truss structures are mainly aligned in a bi-direction mode, i.e., parallel to either 45 • or -45 • directions.With this cellular structure pattern, the whole structure shown in Fig. 4(d) exhibits a strong orthogonal heterogeneous characteristic.
Fig. 4(e) and (f) demonstrate two VPVT porous patterns generated from completely-random or partially-random frequency factors, respectively.From a bionic perspective, these two randomly generated VPVT porous structures show good bionic and similarity to the lattice patterns in the dragonfly wing.(e) features a dense pentagram cellular structure, while (f) combines a large rectangle with an irregular porous structure.In addition, the pentagram cellular structure in (e) shows complex variation that resembles a stochastic porous structure.Nevertheless, we declare this is a 'pseudo-stochastic' design process, in which although the frequency factors are completely selected using a random process, the generated VPVT porous pattern is unique and reproducible provided that the same frequency factors are given.On the contrary, real stochastic porous structures may exhibit a certain amount of variation in their final designs even when all design parameters (i.e.porous density and regularity) remain the same [6,43,44].The high repeatability of the VPVT method can ensure consistency and interchangeability of part designs in actual manufacturing, and so this "pseudo-stochastic" yet deterministic design method can provide considerable beneficial in in-dustrial applications.After understanding the relationship between the frequency mesh grid (design variables) and the resultant VPVT geometric pattern, a designer is able to intuitively design the desired porous structure pattern.The VPVT porous structure of case (f) is an example of simulating the dragonfly wing structure, as denoted by Fig. 4(h), which is tessellated by large horizontal structures with aligned rectangular cells and small, irregular, curved cells.To obtain such "dragonfly wing"-inspired porous structures, we first place an uni-axially varying variable frequency grid in the design domain, whose value decreases along the -axis, and then select (0, 30) and (60,30) as two center points to further reduce the frequency factor along the radius direction.Finally, considering the stochastic nature of biological structures, a small random increment is added to each existing frequency factor within the region 0 <  < 40 to enhance the diversity of the porous pattern.In this way, the lattice size, shape and alignment of a typical pattern (h) in the dragonfly wing can be mimicked using the VPVT method (f), correspondingly.In addition, the independent thickness mesh grid has the same design flexibility as the frequency mesh grid.Herein, a radiusvarying thickness configuration has been applied for the cases of (a), (b), (c) and (d), while a constant thickness is used for cases (e) and (f).The thickness design feature provided by the VPVT method further enhances the design flexibility of porous structures.
It is worth noting that, a complete VPVT pattern may simultaneously integrate several types of these basic units with different characteristics, such as squares, hexagons and rhomboids, as shown in the entire design domain.Depending on their specific locations and manifested functionalities, our proposed optimization routine in the later section will select appropriate types of cellular structures to form an optimal configuration to meet the structural design objectives.Throughout the optimization process, this deliberate adaptability of the lattice geometry not only facilitates accurate emulation of biological structures, but also provides a From a bionic perspective, the two ROI units of (e) and (f) exhibit significant bionic similarity to the lattice pattern found in the dragonfly wing (g) and (h), respectively.solid foundation for achieving enhanced material properties and superior mechanical functionalities.
Based on the proposed design framework, the topology of a complete VPVT structure can be expressed using a level-set function [45].The coordinates of the  ℎ truss bar in VPVT structure can be numerically solved and marked as a Set {( 1, ,  1, ), ( 2, ,  2, )}.With this definition, the material at location (, ) in a specific design domain D  , is defined as follows where the region Ω defined by (x) ≥ 0 is the final topological geometry of the VPVT structure.x = [, ] is the coordinate vector.The term i corresponds to the truss component ID in the VPVT structure.The default topological shape of the i ℎ truss component consists of the region containing a rectangular body and two circular endpoints [45], which are represented by the geometric functions  , (, , ),  ,1 (, , ) and  ,2 (, , ), respectively.  is the thickness of the rectangular body and the diameter of the two circular regions.( ,0 ,  ,0 ) is the center of rectangular body of i ℎ truss component, while ( ,1 ,  ,1 ) and ( ,2 ,  ,2 ) denote the end points of truss and the center of two circular regions.
In particular, the functions  ,1 (, , ),  ,2 (, , ) and  , (, , ) represent the distances of the location x to those three center points, respectively.Obviously, the location x is inside the rectangular or the circular regions if  ,1 (, , ) or  ,2 (, , ),  , (, , ) is greater than 0. The complete VPVT structure is finally attained by using a Boolean merging operation of the above rectangular, circular regions of all truss components, which is achieved by a Max operation of terms   in Eq. (3).

Mechanical property of VPVT structures
This subsection discusses the mechanical property of VPVT porous structures and their corresponding computational method.In this work, a full-scale finite element method (FEM) is used to accurately calculate the mechanical response of the structure under given boundary conditions.Specifically, the original VPVT geometry model is meshed with typical beam FEM elements instead of solid elements to reduce compu-tational costs.The method used to derive the constitutive relationship for the FEM element is based on the Euler-Bernoulli beam theory [46].Once the constitutive matrix D of the beam element is determined, the global stiffness matrix K is assembled using a standard FEM process as given by Eq. ( 4), with which the displacement field U of the VPVT structure is also computed.
where D  is the constitutive matrix of the beam element e, B  is the strain-displacement matrix based on the linear shape function,   is the volume of the element e. F is the corresponding external loading vector,   is the total element numbers.The above FEM-based computation process was completed and implemented in ABAQUS (Dassault.Co.Ltd, ver.19) with a Python script, which is used to control the parameters and the entire computational workflow.The 'B32' beam element with 0.4 unit length mesh density was selected to construct the FE volume mesh model, while the default width of the beam section was set to 0.1 mm.On the other hand, it is worth for designers to consider the effective modulus of the VPVT ROI, as it clearly characterizes the local mechanical behaviour of each VPVT cellular structure (lattice pattern).The effective modulus of some typical VPVT ROIs are evaluated to further demonstrate varying material properties provide a large design space for the proposed VPVT method.The macro-scale effective elasticity tensor,    , for the cellular structure under periodic boundary conditions [34,47] is computed as follows, where  is the volume of a unit cell (ROI) of a porous structure. 0()

𝑖𝑗
is the prescribed macro-scale strain term with respect to periodic boundary conditions (PBC), meanwhile    is the local varying elasticity tensor corresponding to the FEM beam elements.With Eq. ( 5), the compliance tensor   is determined as the inverse matrix of the elasticity tensor   (   in the matrix format).The effective material modulus along each axial direction of the cellular structure is given by Herein, Eq. ( 5) is utilized to determine the effective modulus for some typical ROIs of VPVT structures, and then the results are compared with that of other two recent Voronoi-based porous structures introduced by Do [33] and Lu [6], as depicted in Fig. 5. Considering that the effective modulus of both Do and Lu's Voronoi-based porous structures are derived from the exponential function for open-cell foams [48], the computational results given in Fig. 5 only concern about structures with a volume fraction in 0.1 to 0.5 where the open-cell foams keep a relative high computational accuracy.
Thanks to the capability of generating diverse patterns for porous structures, the VPVT design approach provides a much larger design space than that given by other existing methods.As illustrated in Fig. 5, the Voronoi structures designed by Do and Lu exhibit deterministic, isotropic material properties within a specific volume fraction.This limitation primarily arises from their stochastic and quasi-homogeneous designs on Voronoi patterns.On the contrary, VPVT structures (ROIs), with different values of frequency factors '  ', exhibit remarkable diversity among the shape and size of their Voronoi cells, resulting in a distinct effective Young's modulus even under the same volume fraction.In various designs of VPVT structures, not only the isotropic properties can be readily achieved but also some of them exhibit distinct anisotropic characteristics.For structures like 'VPVT-slender' and 'VPVT-rand', there is a notable disparity between their effective Young's modulus along two axial directions,   ∕  and   ∕  (  is the Young's modulus of basic material utilized in homogenization).For 'VPVT-rhombus' and 'VPVT-honeycomb', their structural patterns approximately exhibit homogeneous, isotropic properties, thus their effective modulus   ∕  and   ∕  are close to each other.For the standard 'VPVT-square' structure, it keeps a strict orthotropic and so its effective modulus values along and directions are equal with each other.Finally, the grey semi-transparent region highlighted in Fig. 5 corresponds to the total design space of the effective modulus values obtained using the VPVT design approach.This is attained through a Monte-Carlo sampling of the frequency and thickness parameters   and   in feasible space.The expansive grey region implies a broad design space provided by the VPVT method, which substantially enhancing the probability of discovering a porous structure with desired performance.

The tailoring ability of VPVT method
The plate with a hole is a typical mechanics problem, which exists widely in many engineering fields, e.g., civil, aerospace, biomedical, and manufacturing etc. Herein, the VPVT method is applied to a structure configuration considering a square plate domain with an inner hole Ω, as shown in Fig. 6(a).The problem can easily be addressed by introducing an extra additional constraint in the original periodic equations, i.e.,  (, ) = 0 over the hole region (, ) ∈ Ω.In doing so, the Voronoi seeding points in the Ω region will be removed, which results in a very sparse VPVT tessellation with few thin rods in this hole region, as shown in Fig. 6(b).These remaining structures within the region Ω can be further removed by applying  = 0.After artificially defining a circular boundary for the region Ω, an appropriate VPVT porous structure design is formed for this square plate with an inner hole problem.Different from conventional element-based topology optimization methods that apply an arrangement of small square elements to fit the circular boundary, an analytical and accurate circular boundary Ω can be defined directly in the point-based VPVT method.Note, although the ROIs near the cutout edges in the VPVT method cannot maintain a complete rectangular shape, it does not affect the Voronoi tessellation, as well as the modelling and optimization process.Therefore, when applying the VPVT method in the design of structures with cutouts, it only requires to introducing additional constraints over the void regions, while the tailoring process for the cutout domain and boundary is analogous to the process of defining an outer boundary.
Fig. 6(d) demonstrated a case of applying the VPVT method to design porous structures within an irregular domain.Herein, an airfoil profile is selected as the outer boundary of the VPVT porous structure.There are two methods to address this problem.The first one is to transfer the Voronoi seeding points generated in the original rectangular domain (e) into the irregular domain (d) using an isoparametric mapping method, and then generate the VPVT porous structure in the aerofoil domain (d).However, this approach may suffer from the issues of large void regions in the final design.Alternatively, in the second method, the VPVT porous structure is generated in the original rectangular domain (e), and then transfer each rod in the design (e) into the aerofoil domain (d) using the conformal mapping method.In doing so, a new Voronoi tessellation well filled the aerofoil domain is generated.The geometric patterns of the three ROIs in close-ups are jointly influenced by both the original structure and the mapping function.Due to the stretching process, these three ROI patterns, which are originally uniform in the rectangular domain (e), exhibit strong heterogeneity in the new irregular domain (d).The detailed computational process is presented in Appendix B.
VPVT also has the tailoring ability to construct metastructures or metamaterials with specific mechanical properties.For example, the VPVT method allows us to create a porous structure or specific ROIs that exhibit negative Poisson's Ratio (NPR) properties.According to the features of re-entrant honeycomb structures [49][50][51], a unit cell ROI with NPR can be generated by introducing additional constraints to frequency  (, ) and thickness (, ) functions in the VPVT method, as shown in Fig. 7(a) ∼ (c) step-by-step.Note, ,  are the normalized coordinates of ,  with respect to the ROI unit.An initial VPVT pattern with a uniform lattice in the ROI is firstly defined, as shown in Fig. 7(a).After applying an additional constrain of the frequency function, i.e.,  (, ) = 0 over a defined region (outside a circle with-out four corners of the ROI), arrow-shaped internal folding components for the re-entrant auxetic structure are emerged, as demonstrated in Fig. 7(b).Subsequently, by performing a thickness tailoring (, ) = 0, as depicted in Fig. 7(c), redundant truss members are removed and eventually a desired unit cell ROI with NPR property is obtained.Note, although the tailoring steps (a) ∼ (c) are separately demonstrated, an NPR-VPVT structure can be directly and rapidly achieved using integrated frequency and thickness constraints in practical applications.
Moreover, the VPVT constructed NPR structure can be modified and optimized, straightforwardly, by applying different design parameters in constraints.The complete VPVT structure can comprise various ROI configurations with different NPRs (or other functionality) while still maintaining a strong connectivity.As illustrated in Fig. 7(d), the Poisson's ratio   in this axuetic structure gradually decreases from −3.01 to −1.35.This example suggests the strong geometric tailoring capability of VPVT method at the ROI level.As such, through the proposed VPVT method, we can systematically design specialized ROIs in various regions to accommodate their unique functional requirements, ultimately assembling them to create a hybrid structure.This design approach is analogous to the concept of digital modular material assembly, which strongly complies with the concept of contemporary 'advanced manufacturing' [52].Fig. 7(d) demonstrates an example.In the human spine, trabecular vertebrae typically require high stiffness to withstand axial compression, while the intervertebral disc (IVD) imparts flexibility to enable spinal bending.However, under axial compression, a lumbar disc herniation (LDH) may occur, compressing the adjacent spinal cord and nerves, eventually leading to severe back pain or dis- ability.The use of an auxetic material is expected to mitigate this effect due to its negative Poisson's ratio [53].In this study, VPVT method is employed to incorporate a functionally graded material (FGM) into the vertebral portion, while integrating auxetic material in Fig. 7(d) for the IVD section.The 'shrink' vertebral is a graded porosity structure for providing sufficient compressive load-carrying capacity, while the IVD region is composed of a re-entrant auxetics structure for preventing LDH.The related FEM simulation result can be referred to Appendix B. Challenges such as the 'curved outline of vertebral bone' and the 'connectivity of adjacent modules' both can be readily addressed by using the VPVT tailoring process.Furthermore, since both types of structures (FGM and auxetic material) are modelled within the VPVT design framework, it is possible to perform the optimization design for different modular parts of the entire implant, simultaneously.As such, the coupling effects between different modules can be automatically tackled during the optimization process.This case study highlights the significant potential of this proposed VPVT method to perform integrated structure design for achieving multiple functionalities.

Inner structure reinforcement
To further demonstrate the benefits of applying variable periodic cells in the design of lattice structures, the VPVT structure is compared with two state-of-the-art porous structures studied in 'restricted cell' mode [38,54] in later Section.3.2.Since both their lattice units contain the inner-layer structure, an inner X-shape architecture is also added in the pure VPVT structure for consistency, as shown in Fig. 8.This will further enhance the lattice performance in terms of stiffness and strength-to-weight ratio.To avoid extremely high computational costs in the later optimization process, herein the angles of X-shaped bars are set to be the same for all Voronoi cellular structures, and therefore, only two design variables are needed for X-shaped bars.To ensure structural stability, the centre points of the X-shaped bars can be positioned either at the Voronoi geometric centre (seeding points) or the gravity centre.However, since this work adopts a bounded Voronoi configuration for the outer-layer framework construction, the centre points of some original Voronoi cells may locate at or even outside of the boundaries of the design domain.Under such a circumstance, only a partial of the internal X-shaped bars is presented in the cellular structure.To avoid this issue, the intersection centre of X-shaped bars is placed at the gravity centre rather than the geometric centre of each cellular structure.Furthermore, for computational analysis and comparison among 'Pure VPVT', 'VPVT with single-orient inner-truss' and 'VPVT with multipleorient inner-truss', please refer to Appendix D.

Optimization process
Although gradient-based methods have high efficiency in solving optimization problems with smooth, differentiable objective and constraint functions, they often suffer from the convergence issues of local optima, particularly when dealing with highly non-convex problems.The SaDE method, on the other hand, provides a global search ability without any requirements for explicit computation of gradient information.In addition, to some extent, the SaDE method can avoid being trapped in local optima.By applying the VPVT method, the original complicated topological design problem for porous structures is transferred into a parametric optimization process, in which three types of design variables, e.g., "frequency mesh grid", "thickness mesh grid" and "angles of X-shaped bars" are considered.These two or three kinds of design parameters will simultaneously influence the VPVT topology in the design domain, which leads to an implicit, non-smooth, non-convex optimization problem.Consequently, in this work, the SaDE method instead of the conventional gradient-based methods is applied to solve the optimization problems for the design of VPVT porous structures.An advanced SaDE algorithm, namely Surrogate-Assisted Classification-Collaboration Differential Evolution (SaCCDE) method [42], is applied due to its low computational costs and good convergence rates.The details of the SaCCDE optimization process are presented in Appendix C.
As the strength-to-volume ratio is one of the most interesting aspects of porous material design, herein the work focus on minimizing compliance performance of the designed structure.The typical optimization problem can be described as follows: find  , , where ( , , ) is the objective function for the compliance of a designed VPVT structure with the parameters of  , , . ( , , ) is the volume fraction of the VPVT structure and   is its upper bound.  and  D are the design domain and corresponding total volume of the VPVT structure.F, U, K denote the external load vector, displacement vector, and structural stiffness matrices of VPVT structure.

𝜙(𝑥)
is the topological design function, as given in Eq. (3).Herein,  = [ 1 ,  2 , ⋯ ,   , ⋯ ,   ] and  = [ 1 ,  2 , ⋯ ,   , ⋯ ,   ] are the frequency factor set and thickness factor set, respectively. = [ 1 ,  2 ] corresponds the two angles of inner X support.t, t, f , f are the threshold of above design variable set.In this work, a 'single-orientation' mode is applied, with which the orientation angles of all the inner support rods in different Voronoi cells share the same values with two design variables,  1 ,  2 .This will significantly reduce the dimension of the optimization problem, and relatively migrate the 'local minimum optimum' issue that often encounters in other angle optimization studies.VPVT cellu-lar structures with diverse-angled inner supports are demonstrated in Appendix D.

MBB numerical example
In this section, a numerical example and the corresponding experiment were carried out to demonstrate the advantages of the VPVT method.A Messerschmitt-Bölkow-Blohm (MBB) beam with a concentrated force loading at its middle point is studied, as shown in Fig. 9. Since this is a symmetry problem, only a half domain of the MBB beam is defined, i.e., 100 mm × 60 mm domain divided into 10 × 6 ROI regions.The total number of design variables of this problem is 156, including 11 × 7 frequency mesh grid factors  , 11 × 7 thickness mesh grid factors , and two angles for the supporting X-shaped bars .Regarding the material properties,  = 1000 MPa,  = 0.3,  0 = 0.05.The volume limitation is set as   = 0.4, while design variable thresholds as t = 0.2, t = 0.4, f = 0.5, f = 4.5.It is worth highlighting that the vol- ume fraction  ( , , ) is a highly non-linear function with respect to all the design variables  , , .Under the above settings, the MBB example leads to much more complex VPVT porosity patterns, which give rise to significant challenges for the global optimization to achieve convergence with respect to the given volume constraint,   .In view of this, after the SaCCDE optimization procedure achieves a relatively stable state approaching but not yet arrived at the final convergence, the obtained design values of  ,  will be fixed and thickness variables  are solely optimized in the later process.After a few more iterations, the final optimization solution with full usage of the volume fraction limit is achieved.As a result, the material volume threshold defined by   can be fully utilized without the need for numerous additional populations and generations applied in the global optimization procedure.This two-stage optimization procedure is marked as dashed line in Fig. 9 convergence history.A population size of 60 per generation is used in the SaCCDE algorithm.
The convergence curve for this optimization design process is illustrated in Fig. 9.During the optimization process, the compliance of the MBB beam showed a remarkable reduction from an initial value of 2.07 to 0.76.Notably, the optimal lattice structure for the MBB example exhibits significant variability and diversity.The pore size ranges from 1.03 to 8.21, and the lattice types in the optimal VPVT structure include square, rectangle, and hexagon configurations, etc.
Furthermore, in this MBB example, a clear material concentration trend was observed, wherein several dense pore cells of diminutive size align along the diagonal axis of the design domain, as depicted in Fig. 10.It appears that the materials are strategically distributed to form a bridging path from the loading point to the fixed vertex at the bottom right corner.This configuration can effectively resist the compressive loading and deformation of the whole structure.The pore shape in this VPVT design also exhibits noticeable changes along the diagonal path.Pore cells near the loading and fixed boundary mainly appear in pentagon or hexagon shapes, as shown by the dashed circle in Fig. 10(a).In contrast, a cluster of lattices in the middle part of the diagonal path is similar in shape to a rectangle, as shown by the solid circle in Fig. 10(a).From an engineering perspective, the material within this dashed circle exhibits relatively diverse principal stress and principal strain directions.As a result, the lattice cells in these areas are predominantly in the form of pentagonal or hexagonal shapes, providing enhanced load-carrying capacity from various directions.In the central section of the diagonal path, where the local deformation involves primarily unidirectional compressive strain along the diagonal direction.As a result, the lattice shapes in this region consist mainly of rectangular units arranged along the diagonal direction, which increases the structural stiffness.Fig. 10(b) shows the von Mises stress distribution of the optimal VPVT porous structure.The primary regions of high stress are highlighted in green, indicating a clear alignment of the high-stress trace along the diagonal path.In line with this observation, the two inner truss rods are optimally angled at −46.5 • and −8.9 • , respectively.
The average of the two angles closely approximates the diagonal direction of the MBB beam, as depicted by red arrows in Fig. 10(a), which provides significant resistance to the high stress and deformation along those directions.

Comparison with other optimization works
To demonstrate the advantages of VPVT, we compare its design performance with the design results obtained from two other optimization methods.One method utilizes Gao's open-access codes [38], while the other one employs Hoang's codes [54].These two methods employ the multiscale and full-scale techniques, respectively, to design 2D porous structures with a 'restricted unit cell' mode.Gao's method applies a simultaneous optimization procedure to design porous structures at two distinct scales.The macro-scale optimization determines the desired material distribution within the design domain, while micro-scale optimization finds the corresponding micro-architectures that can match the desired elastic properties of macroscopic elements.In his work, the actual volume fraction   =   ⋅   , where the macro-volume fraction   and micro-volume fraction   are the parameters manually predefined before the optimization process.In contrast, in Hoang's program, instead of applying the homogenization process, an adaptive geometric component was used to project onto the background element to obtain a continuous density field, which is then used to calculate mechanical response and sensitivity information of the designed model [54].
To ensure an objective comparison, the same 2D deformation model as Gao and Hoang is applied to the optimal VPVT beam structures.Some adjustments had also been made to certain parameters in their methods.Firstly, the 'TRI' element with a 0.25 mm mesh size instead of the 'B32' element is applied to recreate the FEM model for VPVT structures.This allows us to retain the same geometry while omitting the width property.The detailed FEM modelling process for the VPVT structures is presented in Appendix E. Secondly, the design configurations in Gao's and Hoang's codes have been altered to ensure that their lattice sizes match with the variable cavity sizes (lattice sizes) of the VPVT structures, which range from approximately 1 mm to 8 mm.Additionally, Gao's code, which has a multiscale framework, requires the recalculation of volume constraint parameters, including micro-volume fraction and macro-volume fraction.The VPVT structure has an average lattice porosity of approximately 0.5 ∼ 0.6.Correspondingly, we

Table 1
The comparison of the optimal solutions between the VPVT method, Gao et al.'s method [38] and Hoang's method [54].* The 'geometry' column demonstrates the macro-scale optimal structure on the left side with the corresponding magnified lattice unit on the right side.
set the micro-scale lattice porosity in Gao's codes as 0.5 or 0.6, resulting in corresponding micro-volume fractions   = (1 − 0.5) = 0.5 or (1 −0.6) = 0.4, respectively.As the total volume fraction for the MBB example is fixed at  = 0.4, the macro-volume fraction in Gao's code can be calculated as   =  ∕  = 0.8 or 1.For Hoang's case, the optimized structure was also set as an 'X-shape support' type equal to the VPVT method.The detailed values of these new parameters are listed in Table 1.
Overall, the porous structures obtained by the VPVT method show a better compliance performance than those fixed periodic porous structures obtained by Gao and Hoang's methods, as presented in Table 1.The compliance value of MBB.VPVT is about 7.0% -30.4% lower than most other cases, such as MBB.Gao.1,2,34,5 and MBB.Hoang.2.Considering practical application such as implant design, a frequency limit f = 0.5 and density filter is set on the frequency grid of the MBB.VPVT case in section 3.1 to avoid local blank area, but some compliance performance is sacrificed.Nevertheless, the MBB.VPVT still achieve a similar compliance as MBB.Gao.6 and MBB.Hoang.1 which allow for a large local blank area.A better compliance performance can be attained by further decrease f and refine the meshgrid in VPVT method.The difference in performance between VPVT and other methods mainly results from 'porous architecture' and 'cellular (lattice) size'.The porous architecture of different method as shown in the 'Geometry' column of Table 1.Thanks to the irregular and diverse porous architecture at the lattice level, the VPVT method allows the materials to have more flexibility in aligning with the local principal stress direction, consequently enhancing its load-carrying capacity.In terms of cellular size, the VPVT method benefits from a variable Voronoi cell size that is optimized with respect to the set-up, whereas the cellular size for the two methods mentioned is constant and determined by their macro element size, which is manually selected and non-optimized.As demonstrated in Table 1, although the lattice size of the 'restricted unit cell' method is altered to a series of different values, the fixed-periodic structure consistently exhibits lower performance than the variable periodic structure.
The VPVT method requires fewer design variables than the two 'restricted unit cell' methods.The VPVT method is capable of generating a fine and highly diverse porous microstructure by adjusting the design parameters,  , ,  within each macro-scale element (ROI region).Conversely, the 'restricted unit cell' TO methods typically require not only macro-scale but also numerous extra micro-scale topological elements as design variables to determine the micro-architectures.From an optimization perspective, the limited number of design variables used in the VPVT method relatively enhances the convergence efficiency of the optimization process and reduces computational costs.

Experimental validation: 3D printing fabrication and three-point bending test
To validate the effectiveness and accuracy of the optimal VPVT porous structures, 3-point bending experimental testing was conducted on the MBB beam example in this subsection.It is important to note that four optimal VPVT designs obtained in Section 3.1, denoted by points (a) to (d) in Fig. 9, were modified due to limitations of the available 3D printing and 3-point bending test equipment.Specifically, their design domain and truss thickness were reduced to half of their original values, and their height was increased to 2.5 mm.Additionally, the case of 'Gao.MBB.3' in Table 1 was also printed and tested for the purpose of comparison.Finally, five specimens were printed using a stereolithography apparatus (SLA) on an AUTOCERA-R 3D printer under the same conditions, as depicted in Fig. 11(a).The preparation of the model file, including build supports and component slicing, is completed in the 3DXpert software (OQTON Lo.ctd).The material of specimens is PLA resin with an elasticity modulus of 669 MPa and a Poisson's ratio of 0.3.The machining processing parameters utilized are the same as the default configuration given by the printer manual, including layer thickness of 100 μm laser power of 500 mW, scan speed of 10 m/s, and high fill density of 100% mm.The specimens are aligned horizontally on the building plate with an additional 0.5 mm high rectangular support placed underneath.The support partition is later removed by sanding with 40-and 100-grit sandpaper.
Three-point bending tests were performed on Shimadzu AGS-X universal testing machine.Firstly, a pre-loading force 0.1 N is applied at the midpoints of MBB beams, and then each SLA specimen is loaded at a rate of 2 mm/min until failure.Two clamp points are respectively loaded at a 1∕6 distance of the total length from the side edge of specimens.In the meanwhile, a 'B32' FEM simulation with the same boundary conditions, as presented in Appendix E, had been completed for the purpose of comparison.The corresponding linear stiffness is calculated as the external F divided by the displacement U of the loading pin.
The experimental load-displacement curves for the four VPVT specimens, along with the Gao.MBB.3 specimen, have been obtained and presented in Fig. 11(b).A steeper curve during the linear elastic stage indicates a higher stiffness for each structure.The experimental results further confirm the numerical findings presented in Section 3.3, i.e., the stiffness of the optimal design case (d) exceeds that of the other four cases.It can be observed that the stiffness curves obtained from FEM simulation and experiment show relatively high agreement in linear period, as denoted by the star and green lines in Fig. 11(b).The minimal discrepancy between the simulation and the experimental testing implies a high computational accuracy of using 'B32' beam elements instead of a 3D volume model for the FEM simulations presented in Section.3.
Furthermore, the experimental data also demonstrate the superior mechanical performance of the VPVT-designed structures in the nonlinear stage until failure.It is noteworthy that all four variable-periodic porous structures (the four optimal VPVT cases) exhibit higher ultimate damage strength than the fixed-periodic structure (Gao.MBB.3), with case (d) demonstrating the highest strength.The final damage threshold of VPVT(d) is about 47.8% higher than Case Gao.According to the 'redundant load paths' theory [9] and engineering intuition, when a local area of the VPVT-designed material is damaged, the stress within the broken region is redistributed through the nearby intertwined truss network, supported by adjacent 'backup' truss members.This stress-transformation mechanism is advantageous in maintaining the maximum stress at a consistent level before and after initial damage, thus preventing cascading failure or complete collapse of the structure.The load-displacement curves of cases (a) and (d) support this structural design concept to some extent, as the curves display a noticeable inflection point during the middle stages of damage.This phenomenon suggests a temporarily increased local stiffness at the crack tip, likely resulting from the stress transformation mechanism discussed above.This property positions the VPVT method as a promising candidate for the development of advanced fracture-toughness materials and structures.

Discussions
VPVT provides a fast and efficient modelling approach to construct irregular, nature-inspired porous structures through a function-driven mechanism.By varying the 'frequency factor' design variables at four corners of each ROI unit, this novel method is capable of producing a wide variety of lattice structures, covering the most common structures found in nature.While many bio-inspired design methods primarily rely on neural network training or imaging processing techniques [28,55], the proposed VPVT method offers a comparable design space using only a few number of design variables.The method is also highly efficient, as the seeding positions for VPVT are determined from 'analytical solutions of mathematical equations' rather than any 'data-training' process.
From a design perspective, the frequency design variables have a clear physical meaning and exhibit a strong correspondence to the resultant geometric shapes.This enables one to rapidly design desired structures through intuitive understanding.According to our study in this work, a few simple and intuitive cases can be directly generated: (1) a square or rectangular lattice generally needs equal frequency factor pairs on corners; (2) a standard honeycomb requires certain initial phase difference in Eq. ( 1) along  or  direction in ROI.This can be approximated by utilizing a radius-varied frequency field.With the above understanding, the dragonfly wing pattern in Fig. 4 can be directly derived.On the other hand, from the optimization perspective, VPVTbased design framework is also capable of achieving appropriate lattice structures with desired material properties and structural performance by optimizing the design variables.The SaCCDE optimization will automatically search for the appropriate frequency variables to match ideal porous shape and size to improve the performance.
The function-driven seeding mechanism presented in Appendix A can be flexibly adjusted to accommodate practical design requirements.For example, a radial basis function (RBF) interpolation that is often used in meshless methods can be employed to replace the 4-node interpolation, allowing for the design of more complex, curved boundaries within the design domain.Alternatively, additional constraint functions or appropriate phase shift terms can be incorporated into the trigonometric functions of Eq. ( 2) to make specific seeding arrangement for desired functionality.For instance, Section 2.3 presents an example of designing auxetic honeycomb structures by introducing extra tailoring constraints in Eq. ( 2).
Mimicking nature's features to design engineering materials and structures with unique and unprecedented properties has attracted great interest.As the biological structure usually processes certain natural randomness, many state-of-the-art bionic works [4,56,57] applied a stochastic generation procedure to design porous structures with vivid biological patterns.Although the stochastic property can bring a great increase in the geometric diversity of porous structures, it can also lead to problems such as "random defects" and "low reproducibility", which are detrimental in actual industrial manufacturing.In contrast, the VPVT method proposed in this work offers a novel perspective to efficiently parameterize irregular porous cellular structures, thus enabling a deterministic yet pseudo-stochastic approach to bio-inspired structure design.This method can, to a certain extent, avoid the above-mentioned problems of stochastic-based methods.
The proposed VPVT method exhibits robust structural connectivity without additional constraints.Through the application of VPVT, the porosity of the structure is translated into a 'G0' continuous frequency function  (, ) across the design domain (, ).Consequently, adjacent Regions of Interest (ROIs) share identical frequency values  along their interface boundaries, facilitating a seamless transition of Voronoi shapes between ROIs.In contrast, many conventional design methods require a manually seeding operation on the boundaries of ROIs to ensure geometric connectivity between adjacent microstructures [9,30,38,[58][59][60].Furthermore, by further increasing the continuity of the function  (, ) among ROIs, the smoothness of porosity variation and geometric structure in the final model can be automatically improved at the architectural level.
Many of previous works applied a 'SIMP-based' two-level method to design porous [4,6,33] or multi-scale structures [15,34,38], as discussed in Section 2.2 or Section 1.Initially, they match the elastic properties of ROI units with the relative density , which will be optimized with gradient-based methods.Subsequently, the optimized  is de-homogenized to the actual cellular geometry based on design vari-ables, such as the number of Voronoi seeds  and the truss thickness  [4,6,33,34].This optimized  serves as the basis for designing the cell geometry.To ensure a relatively accurate matching relationship between the ROI elasticity tensor and de-homogenized architecture, a very fine design scheme for the cellular structures must be applied in the ROI region [6,34].Although their optimized result may contain different cell shapes such as rectangle and hexagon, these shapes are not directly involved in the optimization process.For example, different shaped cells with the same  of ROIs may be considered as the same elastic property [4,33].In contrast, the VPVT method proposes a novel full-scale approach to straightforwardly construct and optimize the desired cellular architecture.Unlike previous studies that used a single design variable  for each ROI, the VPVT method uses four frequency values  at ROI corners and their varying combinations is capable of creating a wide diversity of lattice structures.This feature also enables a full controllability on the ROI geometry.As  was utilized as basic design variables directly participating in the optimization process, the proposed method can take into account the shapes of lattice cells during the optimization process to achieve desired performance, as demonstrated in Section 3.1.In addition, this design framework also allows us to model porous structures with fewer Voronoi seeds, resulting in greater adaptability to meso-macro scale problems.For example, the design of truss-based buildings, scaffolding where the amount of voids is relatively limited.

Conclusions
This paper proposes an innovative variable periodic Voronoi tessellation (VPVT) method for designing bio-inspired porous structures.The resulting designs exhibit a remarkable resemblance to natural materials and structures, reminiscent of structures like dragonfly wings.This efficiency meets the material requirements for a wide range of mechanical and multifunctional demands, such as porous bone implants and architectural supports.
With this VPVT method, porous structures are designed through a function-driven mechanism, wherein the porosity and lattice pattern are generated from periodic functions and the associated frequency and thickness parameters.In this design manner, a VPVT porous structure may contain a wide variety of lattice types, such as, honeycomb, rhombus, square, and slender, all of which are able to maintain varying lattice size and high-quality connectivity.Furthermore, the functiondriven mechanism provided by the VPVT method promises a strong material tailoring capability, making it suitable for irregular boundaries and modular designs.This forms the foundation for an integrated and comprehensive approach to apply in practical engineering designs, e.g., spinal implant, impact protectors and bridge structures.The variable-periodic porous design enabled by the VPVT method substantially extends the design freedom compared to fixed periodic methods, ultimately resulting in structures with superior performance.
The future work will mainly focus on the following three aspects: (1) The VPVPT method will be further extended for efficient multiscale optimization at the micro-macro scale by incorporating homogenization methods.(2) The VPVT method will be extended to design three-dimensional problems by introducing a three-dimensional seeding mechanism.(3) The VPVT method can be further extended by adjusting the periodic equations or introducing specific additional constraint equations.This modification enables the rapid design and optimization of porous structures with tailored morphological characteristics, precisely meeting the functional requirements of porous structures, such as heat propagation and conductivity.These future works will enable this proposed VPVT design method to further enhance the performance of porous structures and attract wide applications.
(A.1) and (A.3) as To further simplify the equations system in Eq. (A. ) According to Cardano formula [61], the solutions can be divided into 4 cases and expressed as follows , where If   = 0,   = 0,   = 0, it has   = 0.For attaining the complete solution set, the MATLAB program will iterate through all possible pairs of  1 and  2 , which are then substituted into Eq.(A.6) to obtain the potential solutions for  and the corresponding values of  within a given ROI.Subsequently, ,  can be determined by substituting the solutions of  and  into Eq.(A.5).The final solutions for each pair of (, ) in a Cartesian coordinate will be used as the seeds for Voronoi tessellation.Notably, the previous calculation only remains valid under the assumption ' 1 ≠ 0'.In cases where ' 1 = 0,  2 ≠ 0', the corresponding calculation process to compute the final solutions is analog to the one for ' 1 ≠ 0'.When ' 1 = 0,  2 = 0', the solutions are (, ) = (0, 0) or  (, ) = 0. Note, in this work, only non-negative values for  (, ) are considered, due to the physical significance of frequency factors  .Further, we artificially specify that the points (, ) satisfying  (, ) = 0 can never be used as the seeds of Voronoi cells.The 'pore density' is an important parameter for the design of porous structures [6].For a typical VPVT structure, the number of voids in a single ROI equals to the number of Voronoi seeds and thus depends on the corresponding frequency function  .Regarding a location (, ) in the ROI region, the interval of the pore density  in a ROI can be estimated based on Eq. ( 2), as follows, where the definition of above parameters is equal to Eq. ( 2).The, the seed number N can be approximately estimated using the following integral function as,  >> 1. Regarding the examples presented in Section 3, with the threshold of  ∈ [0.5, 4.5] and   = 0.05, the estimated number of pores  will be in the range of [1,31].Notably, this number is significantly smaller than the number of pores in other porous structure [33,57,62], which suggests the design cases presented in this paper are all of the meso-macro scales.

Appendix B. Computation of VPVT tailoring operation
By introducing additional constraints in the original periodic seeding equations, the VPVT framework can be extended to design structures with complex boundaries and geometric features, at both macroscale or micro-scale (a single ROI) levels.This advantageous feature of the VPVT method allows engineers to efficiently perform preliminary bionic designs based on their intuition and experience.Herein, a detailed information is presented to demonstrate the extended tailoring capability provided by the proposed VPVT method.Voronoi seeds in the Ω region are removed after applying the zero frequency constraint over this domain, i.e.,  (, ) = 0 in Ω, while the seeds outside Ω are unaffected.Traditional topology optimization methods often rely on complex topological elements to determine the presence or absence of materials.Therefore, the tailoring normally requires a properly refined mesh to fit exactly within the desired boundaries.On the contrary, the VPVT method surpasses this limitation by tailoring at the point level instead of the element level.For the cases when the size of a ROI substantially larger than the Ω region with inner curved boundaries, the efficient tailoring process of the VPVT design framework will not be affected by adjusting the  ,  functions.This adaptability makes the VPVT method particularly advantageous for industrial designs with curved boundaries or profiles, such as bone implants or advanced vehicle battery electrodes.
For the design example with an irregular boundary as shown in the main manuscript Fig. 6(d

𝑃 𝑖
′ point is placed in a proper location with  1 continuity.To transfer the endpoint coordinates of every rod from the square domain to the airfoil domain, a shape mapping function utilising a high-order element interpolation technique [63] as expressed in Eq. (B.1) is applied: Herein, ,  are the normalized parametric coordinates of the original Cartesian coordinates (, ) within the entire design domain (Note, ,  are not normalized with respect to a single ROI).  ,   ,   ,   are the lower and upper bounds of the entire design domain along  and  directions, respectively.
In above 'Negative Possion's Ratio' example, a key step in the 'frequency tailoring' process is to ensure the existence of the Voronoi seeds at four corners of the ROI, while simultaneously eliminating other Voronoi seeds located near the four ROI boundaries.This results in a truss structure with a 'folded arrow' pattern, which can later be cut into a re-entrant honeycomb structure by appropriate 'thickness tailoring'.In the above process, the length of the 'folded arrow' truss plays a crucial role in determining the Poisson's ratio of the final structure where where  is the interval distance of neighbouring Voronoi seeds in the parametric coordinate space, as shown in where In practical design processes, the negative Poisson's ratio of the structure under small strains can be expressed as a function of .A quick and straightforward approximate method for the calculation of Poisson's ratio , without solving complex equations, is provided.Given that the overall strain is very small, axial deformations of the rods within the system can be neglected.Referencing the work by Gibson [64], the overall deformations and strains of the structure are expressed as follows: In this analysis, rather than employing the effective Young's modulus and the mechanical equilibrium equations, we directly resort to the deflection formula for a cantilever beam, to estimate the relationship between   0 ,  1 ,  2 .Considering of truss members including  1 ,  0 with the 'P1' point clamped at the base under a specific moment, we deduced the deflection ratio:   0 ∕  1 ≈  0 ∕ 1 ⇒   0 ≈  0 ( 1 ) 2 (B.7)Similarly, we assume  2 ∕ 1 =  2 ∕ 1 .Substituting above equations into (B.6), the Poisson's ratio is as: The discrepancies in the prediction of Poisson's ratio between the above analytical method and the results from FEM simulation (with   = 0.01 of the whole ROI) are presented in Table B.1.The results suggest good accuracy of the analytical method when the values of  are small, while the accuracy is decreased for large .This reduction is primarily attributed to the amplified bending effects of an elongating  0 as  increases, which render our equation (B.8) less accurate.Nevertheless, the examples we have presented highlight the potential of the VPVT tailoring technique in designing auxetic honeycombs.By selecting appropriate values for the feature parameter, 'k', it becomes feasible to achieve a desired negative Poisson's ratio.In the future, our emphasis will be on developing a more diverse range of auxetic structures and improving the computational method of their mechanical properties.
The proposed VPVT method is able to construct a more comprehensive, hybrid structure in a single space by operating ROIs as packaged modular units.The coupling effect and connectivity requirement of different ROI modules can be considered and easily fulfilled by this VPVT method.A FEM test of the designed implant is proposed in Fig. B.2(f), it is obvious that the hybrid structure exhibits a concave curved profile at the location of the intervertebral disc.This profile effectively alleviates lumbar disc herniation in patients.Although, in this example the aforementioned procedure remains primarily at the design stage, the  ,  mesh factor can still theoretically function as optimization variables for subsequent SaCCDE optimisation calculations.

Appendix C. SaCCDE optimization
In this work, a global optimization method, named Surrogate-Assisted Classification-Collaboration Differential Evolution (SaCCDE), is applied to solve the VPVT design problems.Herein, the detailed workflow of this SaCCDE method is presented in the following Algorithm 1.In the SaCCDE method, the population is divided into two subpopulations, i.e., good and poor populations, based on the feasibility rules.Subsequently, a classification-collaboration mutation operation is applied to generate promising mutant solutions from both good and poor populations.Then, the surrogate model is used to identify the most promising offspring solutions, which can speed up the convergence of the optimization process.To address the reduced population diversity caused by excessive greedy information from classified solutions, an adaptive global search algorithm is applied.This algorithm dynamically adjusts the mutation operation based on iterative information, enhancing its ability to thoroughly explore the entire search space.As such, a balance between local and global searches is achieved by applying this SaCCDE algorithm.
The parameter 'population size per generation' in the SaCCDE method significantly influences the effectiveness of the optimization process.Herein, an initial test is conducted to assess the impact of population size on the optimization results, as presented in Table C.1.The final optimal compliance is compared for configurations with different population sizes per generation, while maintaining an approximate total population of 1600.It was observed that 60 populations per generation can be a balanced choice considering effectiveness and computational cost.

Appendix D. Computational analysis of various types of VPVT structures
In this study, we employed a Variable-Periodic Voronoi Tessellation (VPVT) structure with uniform truss orientation across all Voronoi cells.However, it is noteworthy that our framework allows for the generation of alternative designs, such as a pure VPVT structure without inner truss or VPVT with varying truss orientations among different cells.In terms of optimization, lattice units (ROI) with higher anisotropic properties are expected to yield superior solutions, particularly concerning compliance objectives and single loading boundaries [19].
Nevertheless, the application of 'varying truss orientation' lattice will significantly lead to a significant increase in the number of design variables in SaCCDE, potentially causing non-convergence issues in the optimization process.A comparison of the three VPVT mesh types is presented in Fig. D.1.Due to constraints on the number of design variables in SaCCDE (less than 180), we reduced the design domain to a 4x4 mesh grid.For attaining more distinct and comparative optimized solutions, we apply a 'cantilever' boundary conditions rather than 'MBB' boundary here (for small square design domain under MBB boundary condition, material for three different configurations will all concentrate in diagonal region and induce high similarity).The 'multi-

Fig. 1 .
Fig. 1.Porous materials and structures in nature and civil constructions: (a) Gradation model of soil porous material, exhibiting a gradual variation in soil particle density and porosity from the top to the ground; (b) Artificial aluminum foam, a typical functionally graded material with a high strength-to-density ratio.(c) Variable-periodic porous structure of bone trabeculae in the skeletal system, providing high strength and impact resistance.(d) Variable-periodic porous structures construction for the well-known Bird's Nest national stadium in Beijing.(e) Variable-periodic porous structural frameworks in dragonfly wings.

Fig. 2 .
Fig. 2. Comparison of artificial (a) and natural (b) porous materials.(a) was obtained by rerunning an advanced multiscale optimization code, as referenced in[38].The design of artificial porous materials appears monotonous due to the influence of the predefined macro-scale element size on the final design pattern.On the other hand, real biological porous tissues exhibit immense diversity, with lattices shaped like rectangles, pentagons, or even filled.

Fig. 3 .
Fig. 3.The process of designing a VPVT cellular structure.(a)∼(c) demonstrate the seeding and Voronoi forming process in a single VPVT ROI unit (indicated by the black dashed square) refer to design variables 'frequency meshgrid': (a) the basic design domain and corresponding frequency interpolation mesh grid; (b) the solutions of the seeding equations in a single ROI; (c) a Voronoi scaffold (porous structure) generated from seeding points in (b).(d)∼(f) illustrate the thickness assignment of VPVT structure refer to 'thickness meshgrid':(d) the thickness mesh grid; (e) thickness interpolation calculation at truss middles points; (f) the thickness assignment of the wireframe model.

Fig. 4 (
Fig.4(b), a radially varying functionally graded material is designed using the VPVT method.In this cellular structure, the size of the porosity decreases along the radius direction from the origin, and the unit cells in some ROIs exhibit typical honeycomb hexagonal structures.Fig.4(c) shows a typical square cellular structure, which is obtained by assigning constant frequency factor in the VPVT design method.Fig.4(d)shows a novel diagonally concentrated porous structure, in which the materials tend to cluster in the diagonal region, resulting in dense and small cavities at the sites.The corresponding ROI image (in the lower right box of Fig.4shows that the design domain is enriched by long slender-shape cellular structures, in which the truss structures are mainly aligned in a bi-direction mode, i.e., parallel to either 45 • or -45 • directions.With this cellular structure pattern, the whole structure shown in Fig.4(d) exhibits a strong orthogonal heterogeneous characteristic.Fig.4(e) and (f) demonstrate two VPVT porous patterns generated from completely-random or partially-random frequency factors, respectively.From a bionic perspective, these two randomly generated VPVT porous structures show good bionic and similarity to the lattice patterns in the dragonfly wing.(e) features a dense pentagram cellular structure, while (f) combines a large rectangle with an irregular porous structure.In addition, the pentagram cellular structure in (e) shows complex variation that resembles a stochastic porous structure.Nevertheless, we declare this is a 'pseudo-stochastic' design process, in which although the frequency factors are completely selected using a random process, the generated VPVT porous pattern is unique and reproducible provided that the same frequency factors are given.On the contrary, real stochastic porous structures may exhibit a certain amount of variation in their final designs even when all design parameters (i.e.porous density and regularity) remain the same[6,43,44].The high repeatability of the VPVT method can ensure consistency and interchangeability of part designs in actual manufacturing, and so this "pseudo-stochastic" yet deterministic design method can provide considerable beneficial in in-

Fig. 4 .
Fig. 4. Various VPVT structures under different frequency and thickness grid set-ups.For clarity, the line width of the images is appropriately scaled.(a) a unidirectional functionally graded structure; (b) a radius-directional functionally graded structure; (c) a square cell structure; (d) a diagonal anisotropic structure; (e) and (f) are porous structures generated by random mesh grid configurations.ROI units can be classified as: (a) rhombus, (b) honeycomb, (c) square, and (d) slender, or as irregular ones (e) and (f).From a bionic perspective, the two ROI units of (e) and (f) exhibit significant bionic similarity to the lattice pattern found in the dragonfly wing (g) and (h), respectively.

Fig. 5 .
Fig. 5.The effective modulus of typical VPVT structures is compared to the other two Voronoi-based isotropic porous structures [6,33].Do's and Lu's designs exhibit isotropic properties, so their effective modulus remains constant along different directions.Meanwhile, the VPVT structure may possess either isotropic or anisotropic properties.The data line for each type of VPVT structures is attained by gradually increasing the thickness values (, ) while fixing its frequency values  (, ).The total effective modulus range for the VPVT structure is determined using a Monte Carlo sampling process of  (, ) and (, ) in design space.This is shown as the semi-transparent gray area in the image.

Fig. 6 .
Fig. 6.The tailoring process of the VPVT method by adjusting  (, ) and (, ) functions for the designs that involve domains with either internal or external irregular boundaries.(a)∼(c) show the tailoring process of the VPVT method in steps to the design of a square domain incorporating an inner circle boundary Ω.(d) presents the transforming of the original VPVT cellular structure into an irregular aerofoil by quadratic eight-node isoparametric mapping.

Fig. 7 .
Fig. 7.The 'negative Poisson's ratio' tailoring process and modular design for ROI units.(a)∼(c) shows the tailoring process by the VPVT method in steps to build an auxetic honeycomb structure; (d) Revised structure with distinct negative Poisson's ratios; (e) A design of human spinal implant with distinct VPVT modular structure, e.g., FGM and NPR modules.The close-up on right side presents the details of the connection between two distinct ROI modules.The human spine demonstration in (e) is re-graphed from [53].

Fig. 8 .
Fig. 8.The process of adding inner X-shape supporting in pure VPVT porous structure.(a) the gravity center of each lattice cell; (b) inner X-shaped bars reinforcement at the gravity center of each lattice cell; (c) the final design of a modified VPVT porous structure in a full domain.

Fig. 9 .
Fig. 9.The optimization convergence curve of the VPVT design process for the semi-symmetry MBB beam under 4320 population (4320/60 =72 iteration in SaCCDE).Four marked points in the convergence curve are selected and the corresponding FEM displacement results are illustrated in (a) ∼ (d), respectively.(e) and (f) are the optimized frequency variable field and thickness variable field of case (d).

Fig. 10 .
Fig. 10.The truss alignment and stress field of the VPVT optimized result for the MBB beam problem.(a) the geometry model; The direction of the trunk in the TO solution and the inner truss angles  are denoted in (a).(b) illustrates the von Mises stress distribution in the optimal design of VPVT porous structure.The primary direction of the high-stress region aligns with the diagonal of the design domain.

Fig. 11 .
Fig. 11.Five 3D-Printed experimental specimens: first four are the VPVT designs for the MBB beam example, as illustrated by cases a∼d in Fig. 9, and the last one is of the periodic lattice structure design based on Gao.MBB.3 in Table 1.The upper row of image (a) shows the overview of specimens before experimental testing, while the images in the bottom-hand show the failed specimens after the three-point bending test; (b) Loading-Displacement curves from the three-point bending tests; (c) is the corresponding linear FEM simulation model and result for MBB.VPVT.d.
Fig. B.1(a) and (b) illustrate the distributions of frequency factor  (, ) and thickness factor (, ), respectively, which are applied in the VPVT method for the design of a porous structure with a central hole (Fig. 6(a)∼(c) in manuscript).As demonstrated in Fig. B.1(a), the ), the corresponding relationship between the original and post-mapping domains is illustrated in Fig. B.1(c).Since three points at least are needed to define an edge, eight vertices named   and   ′ are used to define the original and post-mapping domains, respectively.A smooth airfoil boundary can be generated once each

Fig. B. 1 .
Fig. B.1.The tailoring process on  (, ) and (, ) for the design domains with internal or external irregular boundaries.(a) and (b) show the frequency function  (, ) and the thickness function (, ) for the tailoring of a rectangular plate with an inner hole example, respectively.(c) shows the relationship between the original frequency function and the post-mapped one for the design with an airfoil outer boundary.The corresponding control nodes along the boundary are marked and denoted as [ 1 ,  2 ,  3 .. 8 ]. design.By adjusting the parameters during the 'tailoring process' over different ROI units, it is possible to create 'folded arrow' trusses with varying lengths, resulting in auxetic structures with different Poisson's ratios, as shown in the Fig. 7(d).Herein, further details are presented of this modelling process.Fig. B.2(a) suggests the basic square cell for the tailoring operation.In particular, the initial frequency function  (, ) is increased to 6, different from the previous value of 3.This adjustment results in a sufficient density of Voronoi seeds in the ROI domain, thus providing a highly variable tailoring flexibility for the ROI unit.Fig. B.2(b) highlights the difference in Voronoi seed distribution and truss member placement in three types of ROI units after the frequency  (, ) and thickness (, ) tailoring.The blue dots are the remaining Voronoi seeds after 'frequency tailoring', while the red bar symbolises the corresponding Voronoi pattern generated by these seeds.The black, bold bar represents the final re-entrant honeycomb material after 'thickness tailoring'.For the purpose of clarity, we introduced a feature index  with values 2, 3 and 4 to denote three cases.Thus, the location of the important geometric element in the ROI can then be parameterised in terms of , as shown in Fig. B.2(c).In view of this, the fitted frequency function  (, ) can be expressed as follows: Fig. B.2(c).The definition of parametric thickness function (, ) is as follows: { (, ) = 0.2, in (, ) ∈ Ω ′ ; (, ) = 0, the other region.

B. 5 )Fig. B. 2 .
Fig. B.2.The detailed information for constructing auxetic honeycombs with varying negative Poisson ratios.(a) the basic ROI unit built by the VPVT method for tailoring; (b) Three kinds of ROI units with distinct negative Poisson's ratios, where the index 'k' is used to differentiate and identify each case; Red and black truss members compose the revised Voronoi pattern after  (, ) tailoring, while the black truss members denote the final auxetic structure for different ROI units after the thickness tailoring (, ) based on revised Voronoi pattern; (c) The parameterized geometry information for the structure for defining tailored function  (, ), (, ); (d) The demonstration of the name, location of each beam in Equation (B.5); (e) The deflection angles of the beam components in the three types of ROI units under a 0.01 tensile strain in the -direction are presented.The deflection of middle bridge beam in '=4' case is quite large due to its big spanning; (f) shows the FEM testing for the hybrid structure for human spinal implants.The concave curved profile suggests that the auxetic component can effectively alleviate lumbar disc herniation (LDH).

Table B . 1
The error of negative Poisson' ratio between above estimated computation and FEM simulation result.

Table C . 1
The effectiveness and time costing for different 'population size' configurations in SaCCDE.