A numerical approach for designing functionally stiff triply-periodic-minimal-surface structures

. Biomedical implants require stiffnesses matching those of the surrounding bone to avoid stress shielding. Triply-periodic-minimal-surface (TPMS) structures have shown promising characteristics in preventing stress shielding; however, they are limited in allowing for anisotropic stiffness properties that are typically inherent in bone. This paper presents an approach to simplifying a TPMS structure so that common geometries can be used to approximate it. The wall thickness is varied until the desired functional and directional stiffnesses are determined through finite element modelling. Validation of the finite element models is provided through the compression testing of laser powder bed fusion (LPBF) produced specimens. The displacement response to compression testing is presented along with the final structure, which closely matches the stiffness of bone. This approach has the potential to increase implant longevity and improve the lives of implant recipients.


Introduction
In the medical field, bone implants are used to great effect when treating patients with bone trauma. Such trauma can vary from minor fractures to complete bone loss due to severe accidents or bone cancers. In the latter case, when a large portion of bone is lost, a bone-like replacement structure is needed to perform the structural function of the lost bone. This can be achieved using titanium three-dimensional (3D) printed bone implants. The use of 3D printing allows for intricate implants to be created on demand according to a patient's specific needs. It also allows for purpose-built interconnecting pores which improve the growth of healthy bone throughout the structure of the implant [1]. This bone ingrowth is termed osseointegration.
One of the factors that must be considered during the design of such implants is the stiffness of the implant. Human bone consists of two main types, namely cortical (dense) and cancellous (spongy). The exact stiffness of each type of bone varies from person to person, but studies have shown that it averages around 18.2(±0.85) GPa and 1.08(±0.86) GPa in the longitudinal direction, respectively, and 11.7(±1.01) GPa and 1.08(±0.86) GPa in the lateral directions, respectively [2]. It is necessary to create functionally and directionally stiff implants that match the stiffness properties of bone to mitigate the stress shielding phenomenon, in which the surrounding bone degrades due to the high stiffness of the implants. By lowering implant stiffness to match that of the surrounding bone, this stress shielding phenomenon can be reduced [3].
Current solutions to lower implant stiffness use lattices consisting of nodes and struts with relatively sharp corners which often result in stress concentrations [4]. These stress concentrations can act as crack initiation sites which lower fatigue performance and in turn lower implant longevity. Triply-periodic-minimal-surface (TPMS) structures have shown promising stress reduction and osseointegration characteristics due to their continually curving structure [5].
The complex mathematical nature of TPMS structures makes it difficult to modify the structure to achieve functional and directional stiffness. The aim of this article is to present a numerical approach for approximating a TPMS structure using simplified geometries and modifying the resulting structure to achieve functional and directional stiffness similar to bone stiffness. This was done by varying the structure's wall thickness using optimised values obtained from finite element analysis (FEA). The resulting model was then validated by manufacturing specimens from Ti-6Al-4V using laser powder bed fusion (LPBF) and conducting compression tests to obtain the force vs displacement response.

TPMS selection
TPMS structures range from extremely complex surfaces with patterns that are difficult to identify to simpler surfaces with more easily defined patterns. The more complex TPMS structures such as the Fischer-Koch S more closely match the intricate structures observed in bone, which are characterised by small undulating voids. These structures, however, are challenging to modify and manufacture and prove difficult in removing un-sintered titanium powder when compared to more basic TPMS structures. Structures with more easily defined patterns are easier to model and manufacture to create functionally stiff implants.
The Schwarz Primitive TPMS is one such example of a more easily defined surface. It has large openings connecting interior voids, and consistent wall thickness. Large openings between voids increases blood flow through the structure which promotes osseointegration while the consistent wall thickness more evenly distributes stresses throughout the structure. It was for these reasons that the Primitive TPMS was chosen for this research.

Structure simplification
Due to the complex mathematical nature of TPMS structures, it was necessary to develop a simplified model using common geometries that could accurately represent the stiffness response of the Primitive structure. Variations of the simplified model could then be used in conjunction with each other to create a structure that would mimic three-dimensional functional grading of the original TPMS structure.
The degree of simplification is inversely proportional to the accuracy of a model. Therefore, it was necessary to find the model with the best trade-off between accuracy and simplicity to effectively represent the Primitive TPMS structure. As a point of departure, Primitive unit cells ranging from 10 % to 50 % density in steps of 10 % were created using 370, 01002 (2022) https://doi.org/10.1051/matecconf/202237001002 MATEC Web of Conferences 2022 RAPDASA-RobMech-PRASA-CoSAAMI Conference the MSLattice [6] software which allows unit cells of TPMS structures to be created and exported as .stl files. These .stl files were used to take measurements of various aspects of the Primitive structure. An example measurement set is given in Table 1 for the 30 % density Primitive TPMS. These measurements were subsequently used in the simplified model design process as is seen in Figure 1. Throughout this text, the word 'density' refers to the relative density rather than the material property. Similarly, this paper will refer to 'combining densities' which is equivalent to varying the wall thickness according to the measurements taken from the Primitive TPMS structure. All dimensions are in millimetres (mm) unless stated otherwise. From the numerous simplified models that were created during initial development, two models were selected. These were selected based on their simplicity and ability to be modified to achieve functional stiffness. These two models have been called the low complexity (LC) and high complexity (HC) models.
The LC model was created by starting with a 50 mm square to act as a framework. Two concentric circles were drawn on each corner with radii corresponding to the 'C' and 'A + C' measurements shown in Table 1. This framework is shown in Figure 1. The quarter circles created by the intersection between the concentric circles and square were extruded in the third dimension and this process was repeated in three planes. The excess material was removed which resulted in the final LC model. The HC model used the same technique, but rather than extruding the shape in the third dimension, the shape was revolved around the centreline of the square. This process was once again repeated in all three planes and extra material was removed. These two models, along with the intermediate design steps, are shown in Figure 2 with the LC model corresponding to Figure 2(a) and the HC model corresponding to Figure 2(b). Each model is oriented such that the y axis represents the longitudinal direction, and the x and z axes represent the lateral directions. This convention is shown in Figure 2(c) and has been carried throughout the paper.

Finite element analysis of TPMS models
The overall research aim was to achieve anisotropic stiffness. Prior to developing the anisotropic model, a single isotropic model needed to be chosen which best represents the TPMS structure. As a point of departure, two TPMS models, as shown in Figure 2, were created with isotropic stiffness. The two simplified models along with the original TPMS structure were evaluated using MSC Apex (Hexagon Manufacturing Intelligence, North Kingstown, Rhode Island, United States) to determine the predicted stiffness values of each. Material properties for LPBF produced grade 23 Ti-6Al-4V were used when analysing the TPMS models, as this is a material commonly used in biomedical implants. These material properties are shown in Table 2.

Model analysis and selection
Finite element meshes were applied to the TPMS models. All models were constrained to mimic a compression testing rig with a force of 10kN applied in the negative y direction to the entire upper surface of each model. The lower surface nodes were constrained in the y direction and were thus allowed to move freely in the x and z directions, which mimics a physical test setup. The deflection values obtained from the FEA were used to calculate the effective stiffness of the structures using Hooke's law. The effective stiffness of each of the models, along with the deviation of model stiffness from the TPMS stiffness, is plotted in Figure 3(a) and (b) respectively. When analysing the model performance, the aim was to identify the model that most consistently represented the real Primitive TPMS. As expected, the deviation from the stiffness of the real TPMS structure is proportional to the simplicity of the model. Neither of the models appeared to match the Primitive structure's stiffness when considering the plot of the effective stiffness vs relative density, hence the need for Figure 3(b) where the two model's effective stiffness values have been divided by the real Primitive TPMS stiffness values to give the deviation of each model's stiffness from that of the real TPMS. In Figure  3(b), a lower value represents a closer correlation to the Primitive TPMS.
As seen in Figure 3(b), both models show a recognisable trend, where for lower densities, the accuracy and consistency of the models is poor, but as density increases, these metrics improve. Between 30 % and 50 % densities, the stiffness relative to TPMS stiffness remains roughly constant for both models. It might seem of concern that the stiffnesses of the models are more than 50 % greater than that of the TPMS structure. This is likely because the models were created based on the average wall thickness which is not an exact representation of the TPMS structure. This, however, was not of concern since the consistency of the model was of more importance than the exact stiffness values. Cortical bone, as indicated previously, has a longitudinal stiffness of 18.2(±0.85) GPa and a lateral stiffness of 11.7(±1.01) GPa [2]. From Figure 3(a), we note that these stiffness values correspond to Primitive TPMS relative density values of 36 % and 29.5 % while cancellous bone stiffness of 1.08(±0.86) GPa [2] corresponds to a relative density of 9 %. Therefore, the model that would be selected needed to accurately represent the Primitive structure for those densities. As seen in Figure 3(b), both the LC and HC models represent the Primitive structure reasonably well for densities corresponding to cortical bone, but the LC model shows a very large rate of deviation for densities near that of cancellous bone. The HC model remains essentially constant for densities corresponding to cortical bone with a much more gradual rate of deviation at lower densities which correspond to cancellous bone. It should also be noted that the HC model is only marginally more complex than the LC model and thus the trade-off between complexity and accuracy was deemed acceptable. In addition, when considering manufacturability, the large overhangs of the LC model could lead to print failure. The HC model avoids large overhangs due to its more continuously curved surfaces. It is for these reasons that the HC model was selected for further development and testing.

Achieving anisotropic stiffness
Combining multiple densities allowed for anisotropy to be achieved. When considering cortical bone, the stiffness along the x and z directions is the same while the stiffness along the y direction is higher [2]. Since the x and z stiffnesses are kept the same, they are both represented by x from now on. By analysing the stress pattern in the structure, it was found that if a load was applied in the y direction, most of the stress was experienced by the x neck regions. It was suggested that to achieve higher stiffness in the y direction, the neck area between the openings in the x direction needed to be maximized and therefore the x density would need to be higher than the y density.
When multiple densities are combined, the effective stiffness values are altered depending on the densities chosen. Therefore, to determine the density combinations required to achieve the desired stiffnesses, a set of computer-aided design (CAD) models had to be made and analysed. This process is discussed in the following section.

Determining density combination for desired stiffnesses
To determine the exact density combination that would result in the desired cortical stiffness values, models with x density ranging from 32.5 % to 42.5 % in steps of 2.5 % and y density ranging from 17.5 % to 32.5 % in steps of 2.5 % were made. This resulted in a dataset consisting of 35 models each with its own unique x and y stiffness. These datasets are shown in Figure 4 with interpolation between datapoints. For each combination of x and y density values, the structure exhibits a specific x and y stiffness which is represented by the two angled surfaces. Intersecting each surface with a horizontal plane representing the stiffness of cortical bone and projecting these intersections to the horizontal plane resulted in two datasets of density values that achieve the desired stiffness values in either the x or y directions. This is shown in Figure 5.   By plotting the line of best fit for both datasets, it was found that the lines intersect at the x and y densities of 42.15 % and 20.25 % respectively. This represents the combination of density values that satisfy both the x and y stiffness values that were desired for the cortical bone model. Of the models which were used in creating Figure 4, the nearest model is the 42.5 % / 20 % combination which exhibits an x stiffness of 12.3 GPa and y stiffness of 18.5 GPa. This model was therefore used to represent cortical bone. Since cancellous bone is isotropic, the original stiffness data from Figure 3(a) could be used to determine the appropriate density value. The nearest model to the desired 1.08 GPa stiffness is the 10 % density model with a stiffness of 1.38 GPa. This model was selected for the cancellous structure. Figure 6 shows the final cortical and cancellous models along with an example of the printed structure which is described further below.

Model validation
Validation of the finite element model was necessary to ensure the model was created and analysed correctly, and accurately reflected the response of the physical structure. LPBF produced grade 23 Ti-6Al-4V specimens were printed using an EOS M290 (Electro Optical Systems, Krailling, Germany) and compression tested to determine the stiffness of the structure. All specimens were created using 5 mm unit cells and consisted of 448 unit cells making the structure 202040 mm. Specimens were printed with the 2020 mm area forming the base, and thus, the models were rotated such that some of the specimens were printed with the y directions corresponding to the build direction while others were built with the x direction corresponding to the build direction. This ensured that all material properties were measured along the build direction. An example of such a titanium printed structure is shown in Figure 6(c). A total of six specimens were built. After printing, the specimens were annealed at 920 ˚C for two hours in an argon atmosphere followed by a furnace cool. This heat treatment is similar to that recommended by the manufacturer of the M290 for this specific material [9]. The upper and lower specimen surfaces were lightly machined using a lathe to ensure parallel, flat surfaces for compression testing.
Computed tomography (CT) scans were performed to confirm dimensional accuracy. The thin-walled nature of the structure was expected to result in slight print inaccuracies. Specimen dimensions were largely found to be accurate to within 0.05 mm of the computer model. However, deviations of up to 0.2 mm were observed in a few locations which is significant considering that the wall thickness of the specimens ranged from 0.279 to 1.19 mm. Stereoscopic imaging of the specimen surface showed partially sintered globules of Ti-6Al-4V and sharp build faults which may have acted as stress risers during testing.
Testing was performed using an MTS Criterion 44.304 (MTS Systems Corporation, Eden Prairie, Minnesota, United States) fitted with a 30 kN MTS 661.20 force transducer. This testing setup resulted in a force reading error of less than 0.16 % of rated capacity which corresponds to roughly 48 N. A strain rate of 0.15 mm/min was used for testing. This did not follow an ASTM or ISO standard since such a test standard does not exist for this type of geometry. All specimens were loaded and unloaded three times within the elastic region to prove repeatability of the force vs displacement results.
Stereo digital image correlation (DIC) was used for deflection measurement to determine strain. Test specimens were painted with a base layer of white spray paint with the speckle pattern being applied using black spray paint. Two 5 MP LaVision M-lite cameras (LaVision Incorporated, Göttingen, Germany) with 75 mm Double Gauss lenses were used. Loads ranging between 1 and 5 kN were required to sufficiently deform the specimens for DIC measurement. Using the equipment mentioned above, a field of view of 34.9829.18 mm was observed. A noise floor of 0.17 µm was measured. This was done prior to specimen testing to ensure that the testing deflections were sufficiently large as to provide useable results. Test data was analysed using the LaVision Strainmaster 1.3.0 software. The test images were analysed using a subset size of 59 pixels and step size of 19 pixels to limit computational effort while ensuring that any slight changes in speckle pattern from specimen to specimen would not greatly affect the resulting accuracy. These optimum parameters were found during the noise floor tests.

Results and discussion
After testing, the force and displacement data were exported and analysed in MATLAB 2020a (MathWorks Incorporated, Natick, Massachusetts, United States). Figure 7 shows example results for the cortical-based structure with strain measured in the y direction. The measured stiffnesses from the experimental tests have been compared to the calculated stiffnesses and bone stiffnesses in Figure 8.   As shown in Figure 8, the results for both cortical x stiffness and cancellous stiffness closely match the expected stiffness. When studying the data from the cortical y stiffness tests, a large stiffness range is noted.
A few important observations were made during the specimen preparation, testing, and analysis phases. Firstly, the higher stiffness of the cancellous specimen relative to the FEA model stiffness can be attributed to the manufacturing process. Microscope observations showed that the cancellous model had a wall thickness roughly 25 % higher than designed. This is due to the inability of the printing process to achieve the thin wall thickness [9]. Nonetheless, this stiffness value remains within the upper range limit of 1.94 GPa specified in literature for cancellous bone [2] and is therefore deemed to be a successful representation. Specimen ends were uneven due to the cutting process used to remove specimens from the LPBF build plate. An attempt was made to correct this by lightly machining each specimen to flatten the end surfaces. A marked improvement in surface flatness and parallelism was observed; however, it was not possible to ensure perfect alignment of the top and bottom surfaces since the specimen was remounted between machining the upper and lower surfaces. The authors recommend that future researchers ensure both upper and lower surfaces can be machined without the need to remount the specimen. During testing of the cortical specimens with loads applied in the y direction, it was observed that there was relative motion between the upper compression platen and portions of the top surface of the specimens. This occurred in the tests corresponding to the lower stiffness results which suggests that the uneven specimen surface played a significant role in causing the low stiffness measurements. Had these results been removed from the pool of data, the measured cortical y stiffness would have ranged from 16.3 to 17.9 GPa, which is much closer to the desired stiffness. The results presented in Figure 8 clearly indicate that the designed structure was able to achieve anisotropic stiffness similar to that of bone.
Although the structural design process presented in this paper was applied specifically to TPMS structures for use in Ti-6Al-4V biomedical implants, it is not limited to this specific application and may be used for simplification of other structures or for use with different materials and products. One such application may be the design of impact absorbing zones in products such as plastic bicycle helmets. The continuously curving shape and lack of stress rising features would allow for maximum plastic deformation and impact absorption without total structural failure. This would allow impact absorption without the helmet breaking into multiple pieces. This is simply a single example of the wide scope for which this presented work is applicable.

Conclusion
This paper outlined a methodology to simplify the Schwarz Primitive TPMS structure for use in biomedical implants using various geometries and showed that the simplified structure was able to achieve the desired anisotropic stiffness of bone through the combination of multiple densities.
This research contributes to the field of biomedical engineering by increasing the selection size of implant designs to improve medical care. Further research is suggested for the purpose of better understanding the potential applications. Refinement and automation of the TPMS-based structure design would play an integral role for use in custom implants. For the sake of implant longevity, a study into the transition regions between various structure densities such as the interface between cortical and cancellous structures should be conducted. Furthermore, characterisation of the fatigue performance of the structure as well as the bending and shear behaviour of the structure should be done.

Acknowledgements
This work is based on the research supported by the Collaborative Program for Additive Manufacturing (CPAM) and National Research Foundation of South Africa (Grant Number 138432).