A topology optimisation-based design method for 3D Voronoi porous structures and its application for medical pillows

ABSTRACT This study introduces a novel approach to design non-uniform porous structures with gradient density through the integration of the Topology Optimisation (TO) method and the Voronoi porous structure design technique. With the homogenisation method of Voronoi structures, the density data derived from the TO process is converted into seed point distribution for Voronoi diagrams. The porous structure with controlled mechanical properties is constructed based on Voronoi diagrams using the surface mesh superposition method. Compared with uniform Voronoi porous structures, TO Voronoi porous structures exhibit improved strength and stability. The proposed method for generating non-uniform Voronoi structures in this study exhibits notable advantages in terms of simplicity of implementation and robustness. The surface mesh superposition method has advantages in model generation efficiency and accuracy. In addition, the TO Voronoi porous structure design method is applied to design medical pillows, showing significant advantages in shape retention, weight reduction, and personalisation.


Introduction
Additive manufacturing (AM) technology has been rapidly developed over the past decade [1].In virtue of its layer-by-layer manufacturing process, AM can create novel structures with complex geometries [2,3].Structural light-weighting is among the pioneering technologies in modern industry, with the goal of achieving structures with minimal mass and maximum performance [4].Topology optimisation (TO) technique [5] and porous structure are two representative light-weighting strategies.
Porous structures have the potential to replace traditional materials and are increasingly used in engineering applications, such as aerospace engineering [6] and biomedical sciences [7].Specially designed porous structure possesses various advantages, such as a high strength-to-weight ratio [8], negative Poisson's ratio [9] and enhanced energy absorption capability [10].The Voronoi-based porous (VP) structure is a kind of typical random porous structure without obvious repeating units [11].The utilisation of commercial software is an effective approach in the VP structure modelling method.By using the Voronoi tessellation function in commercial software, such as RHINO, processing models with complex shapes, such as bionic bones and shoe soles, is possible to obtain VP structures [12][13][14].In addition, VP structures can also be obtained using the signed distance field (SDF) technique [15].In VP structures, the density distribution is closely related to the selection of seed points in the Voronoi diagram.With this property, the mechanical properties can be directly controlled by changing the seed point distribution [16].
Similar to porous structures, TO techniques present merits in structural lightweight design.Based on the given load cases, constraints, and performance metrics, TO can provide the best material distribution solution in the design space.In the field of TO, numerous TO methods have been developed with varied merits, such as homogenisation based method [17], solid isotropic material with penalisation (SIMP) [18], level-set method [19] and bi-directional evolutionary structural optimisation [20].Recently, Liu et al. [21,22] proposed a subdomain level-set method (VCUT) for the optimal solution of porous structures, each microcell in the subdomain is optimised to improve the structural performance.Based on this method, further exploration is carried out in terms of expanding the degree of freedom of microstructure design [23] and improving the optimisation performance [24].
Various types of porous structures have been successfully applied in TO methods [25,26].In the lattice structure-based TO method, the basic parameters of the lattice cells are used to regulate their relative densities to realise the TO method of the structure [27][28][29].Alternatively, multiple lattice types can be used as the basic lattice to obtain an optimised structure by exploiting the mechanical advantages of different lattice types [30].However, the adjacent cell joints of the lattice can cause abrupt changes in the truss size of the transition part because of density variations that cause local nonsmoothness.The topology based on triply periodic minimal surfaces (TPMS), benefiting from a mathematically controlled topological form and its demonstrated ability to modulate mechanical properties, provides TPMS with a natural advantage in TO applications [31,32].As a type of porous structure, the VP structure demonstrates significant advantages in lightweight design.In comparison to other types of porous structures, the outstanding advantage of the VP structure is reflected in the strong shape-following ability and selfadaptive ability.The advent and development of the restricted Voronoi diagram [33] algorithm along with the centroidal Voronoi tessellation algorithm [34] lead to a uniform and beautiful Voronoi mesh as well as unbroken edges.Voronoi structures exhibit excellent adaptive characteristics, allowing for flexible adjustments of the structure's shape, size, and density based on data distribution.In comparison to traditional porous structures, Voronoi structures are better suited for designing non-uniform structures.Yadav et al. [35] proposed a stress-weighted centroid method to optimise 2.5D and 3D Voronoi porous structures, which eliminate the stress concentration caused by structural inhomogeneity by repositioning materials and pores.Feng et al. [36] proposed a TO algorithm based on differentiable Voronoi representation.The method encodes the discrete Voronoi diagram into a continuous density field in Euclidean space, and the differentiable Voronoi representation can be effectively integrated into the TO pipeline.Do et al. [37] employed a density-mapped processing flow for macroscopic structure optimisation using an optimisation iteration method at the macroscopic level, and 2D Voronoi structure generation and performance analysis at the microscopic level.Similarly, Chen et al. [38] implemented natural frequency optimisation in 2D Voronoi structure optimisation.However, it is worth noting that these studies are limited to the 2D form.Laccone et al. [39] proposed an innovative heterogeneous VP design method that provides excellent mechanical properties in addition to its aesthetic appearance.This heterogeneous VP design method is expected to be a reliable alternative solution to traditional topological connections such as diagonal grids in novel high-rise building designs.
In this study, a porous structure construction method is proposed using 3D Voronoi tessellation and TO theory, resulting in structures with optimised mechanical properties.The topology optimisation Voronoi porous (TO-VP) structure construction method is implemented in four steps, as shown in Figure 1.First, the target structure-optimised density data are obtained using the SIMP-based TO method.The algorithm can handle complex 3D models through secondary development based on the TO function in Abaqus.Second, the optimised density data are transformed into the distribution results of the seed points of the target structure.The mapping relationship between the relative density distribution of the structure and the number of seed points in the Voronoi tessellation is obtained through the modelling experiments of the uniform VP structure.Third, the 3D Voronoi clipping algorithm is used to obtain the Voronoi diagram based on the seed point data.Lastly, the final VP structure is obtained by modelling the basic geometry mesh based on the framework of the Voronoi diagram.The surface mesh superposition method is used to avoid Boolean operations between models, thereby simplifying the modelling difficulty and facilitating the implementation for AM.
The study explores the modelling effect of TO-VP structures under function-controlled seed point distribution and density distribution and demonstrates the processing efficiency and effect of TO-VP structure modelling method in dealing with complex 3D models.The superiority of the TO-VP structure is substantiated in terms of strength and stability through comparative experiments.Additionally, the medical pillow design case demonstrated the potential for personalised customisation with the TO-VP pillow.

Homogenisation of Voronoi structures
A Voronoi diagram is a division of a plane into regions close to a given set of objects.In each Voronoi cell, there is only one unique seed in the set, and any point within the cell has a shorter distance to that site than to other sites in the set.The Voronoi cell corresponding to the seed is defined as follows: where d( • , • ) denotes Euclidean space.Each Voronoi cell V i is the intersection of a set of half-spaces, delimited by the bisectors of the Delaunay edges incident to the site x i .And the clipped Voronoi diagram for the sites x with respect to a connected compact domain V is the intersection of the Voronoi diagram and the domain, denoted as: Each clipped Voronoi cell V i|V is the intersection of the Voronoi cell V i and the domain V.For a more detailed description of the Voronoi diagram clipping method, please refer to the literature [40].
To investigate the mechanical properties of VP structures comprehensively and scientifically, it is necessary to account for the non-uniform features of the structure.Therefore, a study on the mechanical properties of representative volume elements (RVE) at a microscopic scale is conducted based on homogenisation theory.todetermine the macroscopic equivalent properties of the structures.The Gibson-Ashby equation [41] can be used to predict the elastic modulus and mechanical properties of porous materials.For porous structures with low density (0.04 < ρ < 0.5), the following expression is available where r b and r denote the density of solid material and porous material, respectively.E is the relative elastic modulus.The constants C and t can change with the foam group [42].Martinez et al. [43] proved that C = 0.85 and r = 1.95 could be applied for VP structure.Relative density is defined as the ratio of cellular material density to solid material density: where r is the relative density, which represents the average density of the porous solid, V is the actual volume of solid part, and V b is the volume of solid and void part.
In uniform VP (uniform VP) structures, the major factors that affect the relative density of the structure are the basic geometric cell design parameters and the number of Voronoi cells.The radius is the key design parameter for the basic geometry (cylinder and sphere), and it can effectively control the relative density of the structure.The number of Voronoi cells is the same as the number of seed points in the region.The number of Voronoi cells varies with the number of seed points, thus the relative density of the structure is changed.The specified number of random 3D sampling points can be obtained in the cubic space with the Latin hypercube sampling (LHS) method.The cube of length L is used as the sampling space to explore the relative density control method of the VP structure.After obtaining the seed points, the cube is used to perform the clipping operation with the Voronoi diagram and construct the VP structure.
In the modelling of the basic geometry, the manufacturability of the structure will be reduced when the radius is too small, thereby requiring the use of betterperforming equipment to achieve successful manufacturing.If the radius is too large, then the fusion part of adjacent geometric units is too much, and the Voronoi structure characteristics and advantages are lost.The radius of the basic geometric unit is fixed at 0.06L, which makes it easy to achieve manufacturing while reducing the fusion of geometric units.After each acquisition of random sampling points, the volume of the generated porous structure is counted to obtain its relative density.In the design process, the radius of the rods in the cell is initially determined based on the accuracy of the manufacturing equipment and the application scenario.The radius of the rods allows further determination of the cell size.
To reduce the effect of scale effects on the relative density of uniform VP structures, a 3L × 3L × 3L square area is constructed, and the corresponding porous structure is obtained.The relative density versus the number of seed points in a square space of length L is obtained by averaging in a VP structure of length 3L.The experiment has been repeated several times for each group of data, and the mean value is counted.
Under each fixed seed point number, five sets of random seed points are generated to obtain Voronoi porous structures.And the relative density along with the standard deviation of the structure is calculated.As shown in Figure 2(a), the standard deviation values remained below 1.3%, indicating that the statistical analysis of relative density exhibits a small error.As shown in Figure 2(b), the fitted relationship between the relative density of the Voronoi structure and the number of seed points is obtained by the polynomial fitting method.N = 14.9 r + 119.8 r 2 + 63.0 r 3 (0.05 ≤ r ≤ 0.5) (R 2 = 0.99) (5) where r represents the relative density and N denotes the number of seed points.

Structural density optimisation method
TO is a computational design method that automatically generates a structural configuration with maximum performance under design specifications.In this section, the TO technique is used to obtain the density distribution results of the structure.Subsequently, density distribution is used to guide the distribution of seed points in the VP structure.By optimising the density distribution in the design domain, an optimal material layout that satisfies a given performance index can be obtained.The typical design problem is as follows: where U and F are the global displacement vector and loading vector, respectively.K is the global stiffness tensor.K e is the element stiffness tensor after density interpolation.r e is the element density, and V e is the material volume of a solid element.V target is the upper bound of the total material volume.The 3D TO method is improved by MATLAB on the basis of the research of Liu and Tovar [44].After the optimisation analysis, the relative density results are transformed into a distribution of seed points using a sampling algorithm to obtain TO-VP structure.However, the TO procedure method based on MATLAB language has difficulties in dealing with complex models.The commercial finite element software, such as ABAQUS, OptiStruct, and COMSOL, have powerful pre-processing functions that can effectively and stably perform analysis and optimisation tasks on complex 3D models.
In this article, to achieve more flexible model optimisation settings and enhance the optimisation analysis capabilities, the topology optimisation module in the Abaqus software is employed for conducting topology optimisation calculations.The seed point generation function is integrated into the optimisation process using a secondary development program written in Python.

Tetrahedron-based random sampling method
In pre-processing with commercial finite element software, the 3D model needs to be divided into a body mesh, and then the target is analysed and solved.Compared with other types of meshing methods, tetrahedral meshing has better adaptability and robustness in dealing with complex geometries and can provide more accurate and reliable numerical simulation results (Figure 3(a)).To enable the random sampling algorithm to be applied to the weighted random point generation of 3D models, a random point sampling method based on tetrahedral cells is obtained with a tetrahedral mesh as the target.
During the sampling process, the random points generated are required to ensure that they lie within the tetrahedron.As shown in Figure 3(b), in a tetrahedron, an internal point satisfies the following equation where p i is the vertex of the current tetrahedron with the 3D coordinate information (x, y, z). a i is the weight of the vertex p i , which is used to obtain random points.
In the process of generating random points inside the tetrahedron, data are obtained using LHS, and then the random data are processed by normalisation method.The process presented is as follows: where a ′ i indicates the original random data obtained using the LHS, and factor a denotes the normalised scaling factor.Using the random point generation algorithm based on a tetrahedral element, one sampling can  obtain a random 3D point located within the current tetrahedron.
In the TO process based on the secondary development of commercial software, the model is composed of tetrahedral elements.After completing the optimisation analysis, any tetrahedral element has a definite volume and relative density.
As shown in Figure 4(a), the topologically optimised relative density data of the tetrahedral element T e in the target model are r e , and its volume is V e .Based on Equation ( 5), the number of sampling points of the standard volume at the current density is obtained as N r e .As shown in Figure 4(b), the size of the RVE is set to L, and the standard volume is V basic = L 3 .The parameter L is determined by the manufacturing equipment and application requirements.The values of the sampling points are converted according to the ratio of the current tetrahedral volume V e to the V basic of the RVE.
N e denotes the number of seed points that should be distributed in the current tetrahedron.
After obtaining the number of seed points N e of the current tetrahedron T e , the random point generation method is used in the tetrahedron to obtain the specified number of seed points.
N e consists of an integer part and a decimal part.N int e denotes the integer part of N e and N double e denotes the decimal part of N e .As shown in Figure 4(c), in tetrahedron T e , the random point generation algorithm of the tetrahedron is used to obtain 3D random points, the number of which is N int e , as the seed points of tetrahedron T e .For the decimal part N double e , the LHS method is used to generate random data samp double e from 0 to 1.When 0 ≤ samp double e ≤ N double e , a random sampling point is obtained within the tetrahedral element T e using the tetrahedral random point generation algorithm and added to the current sequence of tetrahedral random seed points.
For all tetrahedral elements in the model, the corresponding 3D seed points are obtained using their relative density and volume values, respectively.Within the overall model, the distribution of seed points corresponds to the topologically optimised relative density data.Then, based on the seed points and the original input model, by implementing the 3D Voronoi clipping algorithm reference [40], the clipped Voronoi diagram is obtained, and the TO-VP structure is obtained using the superposition method of the basic geometric mesh.

VP structure design
In the study of Voronoi structures, utilising the signed distance field (SDF) method to obtain the Voronoi structure is a commonly employed approach [30,43].The method has promising robustness, but it requires high computational cost and suffers from model distortion.
To avoid these problems, VP structures are constructed using the surface meshes of basic geometric cells in this paper.Boolean operations, which consume a large amount of computational resources, are not performed during the merging of basic geometric units.The overall structure is obtained by simply summing the surface mesh data of the basic cells.By using this modelling method, the modelling accuracy in porous structures can be guaranteed while reducing the computing resource requirements significantly.
In the truss structure, its basic geometric units include a sphere and cylinder.The surface triangulation mesh results of the cylinder can be obtained by flattening the surface of the cylinder and by using the planar mesh triangulation algorithm.Four typical sphere surface meshing methods can successfully obtain the surface triangular mesh of a sphere [45].As shown in Figure 5(a), the surface triangular meshes of a sphere and a cylinder are shown.The capsule-shaped structure is obtained by combining the surface triangular mesh of a sphere and a cylinder.The structure, although stylistically a uniform model, is still divided into three parts in the surface mesh of the structure, that is, two complete triangular meshes for the spheres and a triangular mesh for the columns.Figure 5(b) illustrates the implementation principle of the surface mesh superposition method.After slicing the model, intersections are computed between the cutting plane and the relevant surface triangles, with the use of triangle topology to establish a closed polyline.Subsequently, the interior and exterior of the polyline region are determined based on the triangle normal vectors, with the interior region being rendered in white.Once all closed polylines have been drawn in the specified region, the correct B/W image processing information is obtained, guiding the photopolymerization equipment in the manufacturing process.In the Voronoi diagram shown in Figure 5(c), each line segment corresponds to a capsule-shaped structure.In Figure 5(d), the results of modelling the porous structure of a Voronoi cell are shown.The cell consists of capsule-shaped units connected in order according to the Voronoi cell topology.The B/W image in Figure 5(d) is a light-cured sliced image of the cell at a specific layer, and the slicing process can be performed in the software ChituBox V1.9.2.
The VP structure can be obtained by constructing all Voronoi units using capsule-shaped structures in the Voronoi diagram.The structure can be manufactured into shape by a light-curing manufacturing process after slicing.When the structure needs to be simulated and analysed, the SDF is constructed using the current structural data, and the unified output model is used to obtain the merged mesh structure required for simulation and analysis.

Additive manufacturing
After completing the model design, the Voronoi model is exported to the Standard Tessellation Language (STL) format and imported into a 3D printer for manufacturing.The manufacturing of the Voronoi structure is performed on a commercial Liquid crystal display (LCD) 3D printer (Photon Mono X, Anycubic, Shenzhen, China).Tough Resin (3D Printing UV Tough Resin, Anycubic, Shenzhen, China) is used to manufacture the structure.In the material parameters of the tough resin, the Young's modulus is 33 MPa, and the Poisson's ratio is 0.31.The density of material is 1.1 g/cm 3 .In detail, the sliced layer height is set to 50 mm, and the exposure time is set to 2 s per layer.

Structure generation effect
In the two classical geometric models of VP structure construction, 100 random seed points are specified.
The VP structure is shown in Figure 6(a), the VP structure accurately reflects the shape of the original structure, the internal mesh transition is smooth and the mesh quality is high.In Figure 6(b), the number of seed points of the VP structure increases with the horizontal coordinate.In this VP structure, the cell distribution is sparse on the left side and dense on the right side.From the left to the right, the number of Voronoi cells shows a uniform gradient variation.In addition, the effect of the density variation on the modelling of the VP structure is shown in Figure 6(c), the input model is a rectangular body with a size of 10L × 10L × 5L.The density control function is r = −0.25 sin p 10L (x + y) − p 2 + 0.25.Figure 6(c) shows that the sparse region of the structure agrees relatively well with the density control equation.
In the 3D model generation process of the VP structure, the combined structure is obtained by using the surface mesh superposition method based on the basic geometry to avoid complex Boolean operations, which can effectively save computational resources.To explore the advantages offered by the geometric the surface mesh superposition approach qualitatively, it is compared with the voxel-based SDF modelling approach.In the following algorithm tests, the comparison is made in terms of generation efficiency and model accuracy.
A specific number of seed points is generated within a square area to obtain a clipped Voronoi diagram, and the VP structure is obtained using a surface mesh superposition method.In the SDF method, a uniform SDF is constructed using data from the basic geometry, and the VP structure is obtained by extracting the iso-surface.As shown in Figure 7(a), the surface mesh superposition method exhibits almost linear time growth as the number of Voronoi seed points increases.For the same number of Voronoi units, the surface mesh superposition approach generates the model faster than the SDF-based modelling approach, thereby reducing the computational time by over 100 times.
Similarly, in the square model, the VP structures with 500 seed points are obtained using the surface mesh superposition method, and the VP structure is obtained   The structure obtained by the surface mesh superposition method and the structure based on the modelling method with SDF (voxel resolution at 0.5) are sliced separately to obtain a single-layer slice map. Figure 7(c) shows the images of the co-located slices of both structures.In the sliced images based on SDF, more fusion is observed in adjacent print regions, which loses model details and leads to distorted structural printing.In the VP structures obtained by the surface mesh superposition method, the slicing patterns are more accurate and can be correctly manufactured by the AM equipment.
In the 3D design task domain, the merging of various components is a common process.The merging of porous and solid structures, especially the interface transition between solid and porous structures, requires intensive attention.The surface mesh superposition method proposed in this paper avoids Boolean operations and saves time and computational resources.By utilising mesh-based slicing techniques, slicing and path planning can be efficiently executed to successfully generate manufacturing files for the merged structure.
The effectiveness of the modelling and slicing methods is demonstrated by interfacial transition experiments with the cylinder model.As shown in Figure 8(a), the cylinder is divided into two distinct regions, with the outer region filled with a Voronoi structure and the inner region with a solid structure.The results of merging solid and Voronoi structures using the surface mesh superposition method are shown in Figure 8(b).At the interface combination of the solid structure and the porous structure, the cylinder rods of the Voronoi structure fit closely to the surface of the solid structure.In Figure 8(c), the sliced image clearly demonstrates the effect of merging the structures at the transition interface.Half of the rods contacting the inner structure are embedded in the solid structure and half are located on the outer side of the solid structure.Finally, model samples are successfully manufactured using the sliced manufacturing files of the structure (Figure 8(d)).

Topologically optimised porous structure and mechanical properties
In topological optimisation research, numerous topological optimisation theories are implemented with MATLAB tools because of its convenient and powerful computational capabilities [44,46].In this section, topological optimisation functions with optimised density distribution functions are implemented based on MATLAB.The 3D Voronoi structure generation method is utilised to model and analyse classical TO cases and study their mechanical performance.
The processing of the TO-VP structure design method is demonstrated using a regular model.The TO analysis is performed on a rectangular structure with the dimensions of 10 mm × 10 mm × 5 mm.The boundary conditions of the structure, as well as the force conditions, are shown in Figure 9(a).The four corner points at the bottom are fixed, and a vertical downward displacement is applied at the centre bottom surface of the structure.As shown in Figure 9(b), the target structure is divided into a 20 × 20 × 10 hexahedral mesh.An optimisation algorithm is written in MATLAB based on Equation ( 6) to iterate the structure for TO.
Figure 9(b) displays the TO outcomes, in which darker gray shades indicate higher relative density close to 0.5.The TO outcome presents elevated relative density in the four and the central area.During meshing, due to the square shape of all grid cells, 3D random points within the cube are achieved using LHS sampling method.The point distribution outcomes are attained through the combined implementation of topologically optimised density data of cells and LHS sampling.As illustrated in Figure 9(c), in the created Voronoi structures, denser Voronoi cells are observed in the bottom corners and top central regions, thereby exhibiting higher relative densities in the corresponding areas and lower relative densities in other regions.The density distribution features of this Voronoi structure agree quite well with the topology of the optimisation results.
To study the mechanical properties of the structure, the MBB beam is selected as the initial structure for structural optimisation design.The dimensions of the MBB beam are 120 mm × 20 mm × 20 mm.For the boundary conditions, the bottom angle of the beam is fixed, and a vertical downward displacement is applied in the middle of the structure.Based on the boundary conditions of the MBB beam, the density data of the structure are obtained using Equation ( 6).Both structures have the same total relative density of 0.29, and the standard rectangular unit defined is 10 mm × 10 mm × 10 mm with a radius of 0.6 mm for the cylinder.
Based on the input model data and density data of the MBB beam, the uniform VP structure and the TO-VP structure are obtained using the Voronoi clipping algorithm and Equation ( 5), respectively.The uniform VP structure and the TO-VP structure are fabricated separately by LCD printer and tested in a universal testing machine for compression.The modelling results and profiles of the uniform VP structure and the TO-VP structure are presented in Figure 10(a,b), respectively.The uniform VP structure exhibits more uniform cell sizes throughout the structure, whereas the TO-VP structure has a non-uniform cell distribution with significantly sparse cells at the upper side edges.In the cross-sectional diagrams of the two structures, this phenomenon becomes even more pronounced.
To assess the mechanical performance of the two structures, a three-point bending comparative test was conducted at a compression rate of 2 mm/min.The experimental data for both structures are presented in Figure 10(d), where repeated tests are performed for the same set of parameters to ensure consistency.The parameter curves for the same set of parameters exhibit excellent agreement.For TO-VP structure, compared to uniform VP structure, the data curve shows a steeper slope in the compression displacement range of 0-12 mm, indicating that TO-VP structure possesses higher structural strength.The peak force for TO-VP structure is 131 N, while for uniform VP structure, it is 81 N.Under the same relative density conditions, TO-VP structure's strength has increased by 62%.Additionally, TO-VP structure experiences its first failure at a compression of 12 mm, whereas uniform VP only undergoes its first failure at 20 mm compression, demonstrating that uniform VP structure has better deformation capability.
Subsequently, the strength properties of the TO-VP structure are analysed in comparison with other classical structures.In the field of lattice structures with repeating cells, the simple cubic structure and the face centred cubic (FCC) structure are chosen for comparison.The simple cubic structure exhibits excellent strength in the axial direction but relatively weak strength in the transverse direction.On the contrary, the FCC structure exhibits balanced strength in all directions.At the lattice design stage, the radius of the cylinders is set to 0.6 mm, resulting in a relative density of 0.29.In addition, using the topology optimisation function of Abaqus, a classically topologically optimised structure with a relative density of 0.29 is obtained.After the structure is manufactured, the three-point bending test is performed using a universal testing machine.
In the compression test, as illustrated in Figure 11(a), it is evident that the TO-VP structure exhibits strength characteristics that are not inferior to those of the simple cubic structure during the compression process.Furthermore, the strength of the TO-VP structure is significantly higher than that of the FCC porous structure.Additionally, a comparison with traditional topologically optimised structures is conducted, as depicted in Figure 11(b).In the three-point bending test, both structures demonstrated similar strength performance.Notably, in the initial stages of compression, the topologically optimised truss structure outperformed the traditional classical topologically optimised structure in terms of strength.
During TO, the forces applied in simulation tests and the specified boundary conditions have a direct impact on the optimisation results of the structure.However, when the structure is applied as a functional device, it will face different force situations.Therefore, it must possess sufficient mechanical strength to ensure that the structure does not fail or deform significantly.The stability of the structure depends on its strength and stiffness, even when the direction of the force changes.To address the structural stability problem, a design scheme based on Literature [32] is utilised.The structure has dimensions of 16 mm × 40 mm × 80 mm, a basic cylinder radius of 0.6 mm, and an overall relative density of 0.29.The boundary condition is set with a vertical downward concentration force applied in the middle of the top surface of the structure, and the bottom ends are fixed.After the TO operation, the obtained seed points are used to generate a TO-VP structure.Furthermore, to facilitate comparison, a uniform VP structure is also generated, with the same relative density and basic cell radius.
The modelling results of both structures are presented in Figure 12(a).After the fabrication of the structures using LCD method, they are tested in a universal testing machine.A fixed load of 10 N is applied, and the downward pressure displacement of the testing machine is recorded.As shown in Figure 12(b), the test sample is fixed with tiltable flat-nose pliers, and the mechanical properties of the structure are obtained at different inclination angles by adjusting the angle of the flat-nose pliers.Tilt angles are tested from 0°to 50°, spaced 10°apart.
The test results are presented in Figure 12(c).Under the same stress, the deformation of the TO-VP structure is only 1 mm, whereas the deformation of the uniform VP structure is twice as much.Moreover, the deformation of the uniform VP structure increased at a faster rate as the tilt angle increased, with a deformation of 6.18 mm at 50°.On the other hand, the deformation of the TO-VP structure increased at a slower rate.Throughout the entire tilt angle test, the compression displacement of the TO-VP structure is significantly lower than that of the uniform structure, which demonstrates the overall structural stability of the TO-VP structure in contrast to the uniform VP structure.Therefore, the TO-VP structure exhibits superior structural stability.

Voronoi porous structure for non-regular models
For the test example, multiple models with complex profiles, including models named fertility, bunny, kitten, dolphin, and squirtle, are selected for TO operations, as shown in Figure 13(a).The density data are transformed into the distribution data of seed points using an automated program developed in Python.Subsequently, the VP structure is obtained using a modelling-based method for such structures.During the TO operation, the bottom of the model is fixed, and the top and the area to be supported apply a vertical, concentrated downward force.The applied concentrated forces are not described in detail because of space constraints.The VP structures are fabricated using a LCD printer.
Figure 13(b) presents the TO-VP structures that correspond to the input models in Figure 13(a), and Table 1 provides statistics on the processing.The TO-VP structure exhibits clear sparsity of Voronoi cells based on the lower side in fertility model.In all models, the number of cells in the TO-VP structure varies in space in a gradient manner.When the density of the TO-VP structure is higher, more Voronoi cells exist in the region, thereby causing the shape to be closer to the original model.Conversely, when the density of the Voronoi structure in this region is low, some features of the original model will be lost in the corresponding region of the TO-VP structure.The clipped Voronoi diagram method also ensures that each Voronoi cell in the Voronoi porous structure is complete and closed, preserving the overall model's integrity.Benefitting from the interconnectivity of the truss in the TO-VP structure, the TO-VP structure can be printed without support if the original model does not require it, such as in the case of bunny, kitten, dolphin, and squirtle.This feature significantly reduces the amount of material used and sample handling.Moreover, uniform VP structures are obtained based on these models using random point generation.These VP structures are fabricated using a LCD printer, as shown in Figure 13(c).The fabricated samples accurately reflected the design details and features of the design model, and the internal structure can be clearly visible with excellent printing results.The fabrication tests confirm that the TO-VP structures possess good manufacturability characteristics.
In the TO method based on Voronoi porous structures, Equation ( 5) is used to convert the density data into seed point distribution data.However, this equation is derived from the statistical data obtained from uniform VP structures.In TO-VP structures, the internal density is heterogeneous.Thus, the accuracy of Equation ( 5) must be verified.To this end, the uniform VP and TO-VP structure design methods are used to construct structures using the five irregular models discussed in this section.The theoretical densities obtained from both design methods are compared with the actual densities to evaluate the accuracy of Equation ( 5).The theoretical relative density is given by the volume constraint in the TO task and has a value of 0.25.The actual relative density is obtained by counting the volume of the Voronoi porous structure.
The relative densities of the uniform VP structures are illustrated in Figure 14(a).For these uniform VP structures, the deviations between the actual relative densities and the theoretical densities are within 0.01, whereas the errors of the model fertility are relatively large.This finding confirms that the uniform VP design method can obtain the VP structure close to the theoretical relative density, as per Equation (5).However, as depicted in Figure 14(b), the error between the theoretical and actual relative densities increases when modelling the TO-VP structure.Apart from the TO-VP structure of the fertility model, the relative density errors of the other structures are within 0.01.The fertility model has larger error between the two design methods because of its complex shape and has the largest surface area-to-volume ratio among all models.To explore the relationship between the surface area-to-volume ratio of the structure and the density error, the ratio is calculated, and the results are shown in Figure 14(c).The density error of the structure has a strong correlation with the surface area-to-volume ratio of the model.A smaller surface area-to-volume ratio of the input model leads to a smaller relative density error, whereas a larger ratio results in a larger error.Moreover, the more details in the model, the higher the relative density error of the structure.

Application in medical pillow
Neck pain is a common problem that affects millions of people worldwide.Medical pillows are widely used to alleviate neck pain and provide support during sleep.However, traditional medical pillows often fail to provide adequate support and may exacerbate the problem.Porous structure medical pillows have significant benefits over traditional pillows because they can provide customised support and pressure distribution [47].The medical pillows based on VP structures have several advantages.The mechanical properties of the VP structure are controlled by regulating the density of the structure.Furthermore, the VP structure pillow has a smooth surface and fits better to the body due to the rational design of the pillow shape.Customised VP structure can be easily obtained by manufacturing the structure with elastic resin, producing a medical pillow that caters to individual sleeping habits.

Additive manufacturing based on flexible materials
In the field of commonly used pillows, a desirable characteristic is a combination of good elasticity and durability to ensure adequate support and comfort for the body.Among photosensitive resins, flexible resins have proven successful in various applications due to their excellent elasticity.In this study, a flexible resin is used as the material to manufacture Voronoi porous pillow.The incorporation of flexible resin aims to endow the Voronoi porous pillow with exceptional elasticity and long-lasting performance.
The evaluation of pillow comfort requires a clear definition of soft and hard areas.Customised Voronoi models of varying densities are designed and   Based on this information, the pressure applied to the pillow can be estimated to be 10,000 Pa.
In the experiment, compressive testing is conducted on the specimens to obtain definitions for the structural soft and hard regions.As depicted in Figure 15(a), five specimens are manufactured using flexible materials, each having distinct densities.During the compressive testing, a maximum pressure of 25 N is applied, considering the cross-sectional area of the specimens.Each specimen undergoes five compression cycles and their  respective compression displacements are recorded for analysis.
As shown in Figure 15(b), during cyclic compression of a single sample, the sample exhibits a good shape recovery ability, and the compression displacement results of the sample are more consistent and concentrated.When comparing the compression displacements of multiple samples, there is a significant shift in the compression displacements as the relative density of the samples changes.As shown in Figure 15(c), as the relative density increases, the compression displacement becomes smaller, i.e. the stiffness of the structure becomes higher.Thus, as an empirical result, a deformation of 1.5 mm was used as a criterion for the regional soft and hard regions.Then, the structure with deformation less than 1.5 mm belongs to the hard region, and the region greater than 1.5 mm is the soft region.

Design of Voronoi pillow
The individual's pillow firmness needs may vary throughout the life cycle of a medical pillow.Using the non-uniform characteristics of the TO-VP structure, multiple functional zones with different degrees of softness are designed into the pillow to accommodate various usage scenarios.As depicted in Figure 16(a), the pillow's usage area is divided into relatively hard and soft areas.The envelope size of the target model is 200 mm × 120 mm × 40 mm.
In the design of the TO-VP structure pillow, the analysis is performed using the topology optimisation plug-in provided in Abaqus.Considering the practical usage scenario, the structural loading conditions are shown in Figure 16(b), where the high pressure region is set to 15,000 Pa and the low pressure region is set to 5000 Pa.In the parameter settings, the optimised volume target is 0.25 and the density range is restricted to be between 0.1 and 0.4.And the cylinder radius of the obtained structure is set to 1 mm.Additionally, a uniform VP structure pillow is obtained with the same relative density and cylinder radius.The design of the Voronoi pillow was realised using a flexible photosensitive resin (ELASTO 1000, Polymer, Suzhou, China) as the manufacturing material.Figure 16(c,d) shows the design and manufacturing results of uniform VP pillow and TO-VP pillow, respectively.Uniform VP pillow exhibits a uniform density distribution with edges that effectively restore the shape of the original structure, indicating a good shape retention advantage.On the other hand, TO-VP pillows reflect optimised non-uniform results.Regions subjected to higher pressures exhibit higher relative densities, while regions subjected to lower pressures have lower relative densities.In addition, the distribution of Voronoi cells is particularly sparse in regions without applied pressure.

Experiment of pillows
In the assessment of structural softness, as depicted in Figure 17(a), compression testing experiments are conducted using a three-point bending fixture.The applied force is estimated to be 5 N based on the fixture parameters to simulate the pressure conditions experienced in the contact area during actual use.The length of the pillow in the design is set to 200 mm.Tests are conducted at 20 mm intervals with a fixed external compression force of 5 N.
In Section 4.1, a threshold of 1.5 mm compression is established to distinguish between the soft and hard regions.As depicted in Figure 17(b), for uniform TO structure pillow, the compression values at different positions fluctuate around 1.5 mm, indicating a moderately firm property.In contrast, for TO-VP structure pillow, there is significant variation in compression displacement among different regions.The compression displacement is lower at the sides, indicating a harder nature, while in the central portion, the displacement is higher, demonstrating a softer characteristic.By comparing the curves of different structures, it becomes evident that TO-VP structure pillow allows for a wider range of softness and firmness variations.

Conclusion
This paper proposes a design method for 3D Voronoi porous structures with density gradients using SIMPbased TO.On a micro level, the relationship between the mechanical properties of the structure and the seed point distribution of the Voronoi-based structure is obtained by using the homogenisation method, hence benefiting from the randomness of the Voronoi structure.On a macroscopic level, the optimal density distribution of the structure is obtained through TO.Subsequently, the density data are converted into Voronoi diagram site data based on the homogenisation method, and the clipped Voronoi diagram is obtained.The VP structure is assembled by superimposing the meshes of the basic geometry to obtain a model that can be manufactured by LCD printer.Through experiments, the advantages of the TO-VP structure are demonstrated from the perspective of modelling, mechanical performance, and manufacturing performance.The primary conclusions are summarised as follows: (1) The method based on basic geometry mesh superposition has high model generation efficiency and excellent model accuracy performance.The VP structure design method has a good lightweight effect and can achieve low relative density model generation.(2) The homogenisation method provides a strategy to control structural properties using Voronoi sites.The clipping operation process of the structure affects the relative density of the structure, and the surface area-to-volume ratio of the model has an obvious positive correlation with the relative density error.(3) The TO-VP structure design method of the irregular model is realised by utilising the secondary development function of commercial software.The robustness of the method is proven by multiple 3D model tests.(4) The TO-VP structure has advantages in terms of strength and stability compared with the uniform VP structure, as demonstrated by the comparison test between the TO-VP structure and the uniform VP structure.(5) The design method shows advantages in the conformability and personalised customisation of medical pillows.The TO-VP structure achieves a wider range of stiffness variation in the pillow's utilisation area, enabling its adaptability to a broader range of application scenarios.

Figure 1 .
Figure 1.Work flow of TO-VP structure design method.

Figure 2 .
Figure 2. Relative density of the structure with varied seed points: (a) statistical analysis of mean relative density and standard deviation, (b) fitted curve.

Figure 3 .
Figure 3. Tetrahedron-based random point generation method: (a) a tetrahedral mesh of a torus, (b) random points within a tetrahedron.

Figure 4 .
Figure 4. Seed point generation process based on tetrahedral mesh: (a) generate tetrahedral mesh based on the target model, (b)obtain the number of sampling points for the current tetrahedron using density mapping method, and (c) obtain the specified number of sampling points in the tetrahedron using the tetrahedral sampling point algorithm and obtain the overall sampling points.

Figure 5 .
Figure 5. VP structure constructing method: (a) surface mesh of the sphere and the cylinder form the pill-shaped basic element, (b) slicing principle, (c) clipped 3D Voronoi diagram, and (d) unit of VP structure and its slice image.

Figure 6 .
Figure 6.VP structure modelling and seed point control effects: (a) uniform VP structures for classic geometries (sphere and torus), (b) VP structure with seed points controlled, (c) VP structure with density controlled.
using different voxel resolutions (number of grids per unit length) to establish the SDF.The true relative density of the target model is 0.296.As shown in Figure7(b), a significant difference with the true relative density of the model observed at lower voxel resolutions.As the resolution increases, the relative density of the structure converges to the relative density of the true model.However, an excessively high resolution leads to more computation time and larger memory usage.

Figure 7 .
Figure 7. Surface mesh superposition method compared with the SDF method: (a) structure generation efficiency, (b) effect of resolution on the structure volume in the SDF method, and (c) comparison of the slice images.

Figure 8 .
Figure 8. Interface between solid and porous structures: (a) the cylinder component is divided into inner and outer sections, with the outer part utilising Voronoi filling, (b) combined solid structure and porous structure, (c) layered slicing images.(d) actual manufactured sample.

Figure 9 .
Figure 9. Voronoi sample structure based on 3D TO: (a) model and boundary conditions, (b) TO results, (c) 3D TO-VP structure.

Figure 10 .
Figure 10.Comparison of mechanical properties of 3D TO-VP structure and 3D uniform VP structure: (a) 3D uniform VP structure and its profile, (b) 3D TO-VP structure and its profile (c) three-point bending test of samples, and (d) force and displacement curves of 3D TO-VP structure and 3D uniform VP structure.

Figure 11 .
Figure 11.Comparison of mechanical properties of 3D TO-VP structure and other classic structures: (a) three-point bending test of TO-VP structure and classic lattice structures, (b) three-point bending test of TO-VP structure and classic TO structure.

Figure 12 .
Figure 12.Stability testing of structures: (a) modelling results of uniform VP structure and TO-VP; (b) the sample is subjected to a compression test in a tillable fixture, (c) the compression displacement of the structure during the angle change.

Figure 13 .
Figure 13.TO-VP structure modelling and fabrication of 3D models: (a) original models, (b) TO-VP structures, and (c) fabricated samples of TO-VP and uniform VP models.
manufactured using flexible materials.Fixed force compression tests are performed to establish criteria for defining the structural properties of the soft and hard areas.The dimensions of the Voronoi blocks are set to 50 mm × 50 mm × 50 mm.Different density schemes are implemented according to Equation (5) to customise the Voronoi block model accordingly.During the generation of the Voronoi structure, the radius of the inner base cylinder is set to 1 mm.Upon investigation, it has been determined that the weight of the human head is approximately 5 kg, and the area it exerts pressure on the pillow is approximately 0.005 m 2 .

Figure 14 .
Figure 14.Density error of porous structures based on Voronoi: (a) comparison of the actual density of the uniform VP structure with the theoretical density, (b) comparison of the actual density of the TO-VP structure with the theoretical density, and (c) the effect of the surface area-to-volume ratio of structure on the density error.

Figure 15 .
Figure 15.Compressive testing of specimens: (a) specimens printed by flexible resin, (b) cyclic compression experiment, and (c) relationship between relative density and compression displacement.

Figure 16 .
Figure 16.Design of porous medical pillows: (a) Voronoi application of medical pillows, (b) force conditions of medical pillows in the optimisation analysis, (c) pillow design result based on uniform VP structure, (d) pillow design result based on TO-VP structure.

Figure 17 .
Figure 17.Softness compression testing of different regions in the pillow: (a) actual compression testing process, (b) compression displacement in different regions of the pillow.

Table 1 .
TO-VP structure generation for models with non-regular shapes.