# A Modeling Method of Graded Porous Scaffold Based on Triply Periodic Minimal Surfaces

Authors: Zhuangya Zhang; Hui Zhang; Yaosong Li; Mingde Duan; Shikun Qin
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In order to solve the design and regulation problems of graded porous scaffolds, the triply periodic minimal surface (TPMS) that based on implicit function was taken as the pore unit for constructing the porous scaffolds. And the TPMS structure was controlled by selecting different types of TPMS and defining the distance function. The transition of different structures of TPMS was performed by the weight fusion method. And the dual contouring (DC) algorithm and octree algorithm were used to divide the solid model of the defect bone into hexahedral elements to determine the target assembly spatial. Through the coordinate transformation of the isoparametric element method, the pore units were mapped into the hexahedral elements to construct porous structure. Finally, a Boolean operation was used to obtain a bone scaffold model with porous structure. The examples showed that the gradient model of various irregular pores could be constructed, and the pore size, porosity, and specific surface area of porous scaffolds could be controlled by adjusting the deformation, distribution, volume of the hexahedral elements, and the structure of the TPMS pore units, which provides a feasible method for the design and regulation of heterogeneous porous scaffolds.


---

## Body

## 1. Introduction

Bone tissue engineering (BTE) is a biological substitute that aims to restore, maintain, or improve tissue functions, which has a good clinical application prospect [1, 2]. Bionic bone scaffold has played an important role in bone tissue engineering, which functions as the architecture for bone growth. The design of bone scaffolds requires that the macrostructure should be consistent with the surface contour of patients’ autologous bone, and the microstructure requires that the scaffolds have micro-porous structures that provide living space and channels for the adhesion, reproduction, nutrition acquisition, and metabolism of seed cells.Bone scaffold can be manufactured using one of many available conventional techniques, such as particulate leaching, phase method, electrostatic spinning method, and gas method. However, these techniques have limitations with incapable of precisely controlling pore size and shape, porosity, and their spatial distributions within the scaffold [3, 4].To overcome the drawbacks, additive manufacturing (AM) techniques were used to prepare the scaffold properties [5–7]. The scaffold was built by the AM techniques with the three-dimensional digital model. However, due to the irregularity of porosity, pore size, and pore distribution of the micro-porous structure of natural bone tissue, it is difficult to effectively describe them by mathematical methods. At present, the digital model design of porous structure was mainly based on computer-aided design (CAD) method, which used the simple sphere, cylinder, prism, and polyhedron to design porous structure [8–13]. However, due to the limitation of the pore shapes, the designed porous structure lacks diversity, especially when designing gradient porous structures, the boundary fusion problem frequently occurred at the contact surface of the adjacent unit, and it is difficult to construct porous structures with gradient variations in pore shape, porosity, pore size, and spatial distributions.Butterfly wings, weevils, crustaceans, etc., in nature have the skeletal structures similar to triply periodic minimal surface (TPMS) [14]. TPMS has advantages including various geometric shapes and can construct parametric mathematical models, which cause researcher’s attention. Rajagopalan and Robb [15] proposed a tissue scaffold design method firstly, constructed a porous scaffold model based on P surface, and prepared them by additive manufacturing technology. Abou-Ali et al. [16] adopted D surface, I-WP surface, G surface, and CY surface as the basic pore unit and prepared nylon thermoplastic scaffolds by selective laser sintering technique, and the mechanical properties of the porous scaffolds with different types of TPMS were experimentally studied. Montazerian et al. [17]and Castro et al. [18] constructed the porous scaffold model based on the P surface, I-WP surface, and G surface, the permeability of the porous scaffold was analyzed via computational fluid dynamics simulation, and the experimental results showed that the permeability of porous scaffold was closely related to the types of TPMS.The natural bone is heterogeneous, which needs to meet the specific requirements of biomechanics and biocompatibility at a specific position. However, a homogeneous porous structure is difficult to become an ideal scaffold due to its single performance. Zhang et al. [19] found that the gradient porous scaffolds were superior to the homogeneous porous scaffolds in comprehensive performance. Han et al. [20] reported that the gradient porous scaffold could optimize the growth space of bone tissue through personalized design because of its varied pore size and porosity. Therefore, the research on modeling methods of heterogeneous porous structure with functional gradient has always been the focus of this field. In recent years, some scholars have also explored the modeling method of gradient porous scaffolds based on TPMS [21, 22]. Liu et al. [23] constructed a gradient porous scaffold model based on D surface and G surface and fabricated by selective laser melting (SLM), and the results showed that the TPMS method was an effective way to achieve gradients in multiple patterns which were comparable to natural tissue with respect to both continuous topology and interconnectivity. Gao and Xiang [24] proposed a design method for gradient porous scaffolds based on G surface. The mechanical properties of the gradient porous scaffold were analyzed, and the results showed that the mechanical properties of the gradient porous scaffolds were better than the homogeneous porous scaffolds under the same average porosity. Although some achievements have been made, there are still many problems in the construction of gradient porous scaffolds, such as the arbitrary arrangement of different pore structures and the precise control of porosity, pore diameter, and specific surface area. The research of gradient structure design method based on TPMS is still in its initial stage. In this paper, TPMS was used as the basic pore unit to construct the gradient porous scaffolds by controlling the pore unit structure and its spatial morphology. The effects of pore unit structure and its assembly spatial configuration on porosity, pore size, specific surface area, and their spatial distributions were studied, which provides a feasible method for the modeling of gradient porous scaffolds.


## 2. Algorithm Flow

The modeling method and general route of the porous scaffold are shown in Figure1.(1)
CT/MRI was used to measure the defect bone, and the solid model of the defect bone was constructed.(2)
TPMS expressed by implicit function was used as the basic pore unit for constructing porous structure. And the pore unit structure was controlled by selecting different types of TPMS and defining the distance function. The transition of different structures of TPMS pore units was performed by the weight fusion method.(3)
The DC algorithm and octree algorithm were used to divide the solid model of the defect bone into hexahedral elements to determine the target assembly spatial. And the hexahedral elements could be regulated by encryption template to meet the design requirements.(4)
Through the coordinate transformation of the shape function, the pore units were mapped into the space units to construct porous structure.(5)
A Boolean operation was used to obtain bone scaffold model with porous structure.Figure 1
Algorithm flow.


## 3. Pore Unit Based on TPMS

### 3.1. TPMS Unit

TPMS is a minimal surface with periodicity in three independent directions in three-dimensional space (the average curvature of any point on the surface is zero), which is usually expressed by implicit surface [25].(1)ϕr=∑k=1KAkcos2πhk⋅rλk+pk,where Ak is the amplitude factor, λk is the periodic wavelength, r is the position vector of Euclidean space, hk is the k-th raster vector in the reciprocal space, and Pk is the phase deviation. The typical TPMS [26]is shown in Table 1.where X=2πx,Y=2πy,Z=2πz,x∈R,y∈R,z∈R. (1) is an implicit expression of TPMS. And marching cubes (MC) method [27] was adopted to construct TPMS, as shown in Figure 2.Table 1
Typical TPMS.
TPMSImplicit surface expressionP surfaceφpr=cosX+cosY+cosZ=0D surfaceφDr=cosXcosYcosZ−sinXsinYsinZ=0G surfaceφGr=sinXcosZ+sinZcosX+sinYcosZ=0I-WP surfaceφI−WPr=2cosXcosY+cosYcosZ+cosZcosX−cos2X+cos2Y+cos2Z=0FRD surfaceφrF−RD=4cosXcosYcosZ−cos2Xcos2Y+cos2Xcos2Z+cos2Ycos2Z=0L surfaceφrL=0.5sin2XcosYsinZ+sin2YcosZsinX+sin2ZcosYsinY+0.5cos2Xcos2Y+cos2Ycos2Z+cos2Zcos2X+0.15=0K surfaceφrK=0.3cosX+cosY+cosZ+0.3cosX+cosY+cosX+cosY+cosX+cosY−0.4cos2X+cos2Y+cos2Z=0G′ surfaceφrG=sin2XcosYsinZ+sin2YcosZsinX+sin2ZcosYsinY+0.32=0P2-DG surfaceφrG=sin2XcosYsinZ+sin2YcosZsinX+sin2ZcosYsinY+0.1cos2Xcos2Y+cos2Ycos2Z+cos2Zcos2X−0.2cos2X+cos2Y+cos2Z=0Figure 2
TPMS. (a)P surface. (b) D surface. (c) G surface. (d) I-WP surface. (e) FRD surface. (f) L surface. (g) K surface. (h) G′ surface. (i) P2-DG surface. Different types of TPMS constructed by marching cubes method.
(a)(b)(c)(d)(e)(f)(g)(h)(i)For the same type of TPMS, the TPMS structure could also be changed by controlling the distancek from the point of φTPMSx,y,z,k to the φTPMSx,y,z [28] as shown in Figure 3.(2)ϕTPMSx,y,z,k=ϕTPMSx,y,z+k=0.Figure 3
FRD unit under differentk values. (a) k = −0.9. (b) k = 0. (c) k = +0.9. The TPMS structure could be changed by controlling the k from the point of φTPMSx,y,z,k to the minimum surface defined by φTPMSx,y,z.
(a)(b)(c)


### 3.2. Fusion Transition of Different Architecture TPMS Units

The transition of different structures of TPMS (including the same type of TPMS unit under differentk values) was performed by the weight fusion method [29].(3)ϕTPMSx,y,z=∑j=1nαjx,y,zϕjx,y,z,where n is the type of TPMS, φjx,y,z is the type of TPMS unit in region j, and αj (x, y, z) is the weight coefficient for region j.(4)αjx,y,z=∑i=1mαijx,y,z=1,x,y,z∈regionj,0x,y,z∉regionj,where αij (x, y, z) is the weight coefficient of the i-th point in region j, i = 1,2, …, m, which could be constructed by the exponential kernel function.(5)αijx,y,z=wij∑i=1mwij,wij=exp−x,y,z−xij,yij,zij2σ2,where (xij, yij, zij) is the coordinate of the i-th point in region j and σ is the width parameter, which controls the local scope of the exponential kernel function.TakeP surface, FRD surface, K surface, and G surface in Table 1 as an example for multi-region fusion transition. Figure 4 shows the fusion transition results.Figure 4
Fusion structure of different types of TPMS. (a) Axonometric drawing. (b) Top view. TakeP surface, FRD surface, K surface, and G surface in Table 1 as an example for multi-region fusion transition.
(a)(b)Figure5 shows the fusion transition results of the same type of TPMS (P surface) under different k.Figure 5
Porous fusion structure of same type of TPMS. (a) Axonometric drawing. (b) Top view. The fusion transition results of the same type of TPMS unit (P unit) under different k.
(a)(b)


### 3.3. Construction of TPMS Point Unit

In this paper, the algebraic method was used to construct TPMS point unit by intersecting the TPMS with the plane perpendicular to the coordinate axis.(6)ϕTPMSx,y,z,k=ϕTPMSx,y,z+k=0,x+=π,x−=−π,y+=π,y−=−π,z=π,z=−π.The FRD,P, and G point units constructed are shown in Figure 6.Figure 6
TPMS point unit. (a) FRD unit,k = 0. (b) P unit, k = 0. (c)G unit, k = 0. Using algebraic method to construct TPMS unit by intersecting the TPMS unit with the plane perpendicular to the coordinate axis.
(a)(b)(c)


## 3.1. TPMS Unit

TPMS is a minimal surface with periodicity in three independent directions in three-dimensional space (the average curvature of any point on the surface is zero), which is usually expressed by implicit surface [25].(1)ϕr=∑k=1KAkcos2πhk⋅rλk+pk,where Ak is the amplitude factor, λk is the periodic wavelength, r is the position vector of Euclidean space, hk is the k-th raster vector in the reciprocal space, and Pk is the phase deviation. The typical TPMS [26]is shown in Table 1.where X=2πx,Y=2πy,Z=2πz,x∈R,y∈R,z∈R. (1) is an implicit expression of TPMS. And marching cubes (MC) method [27] was adopted to construct TPMS, as shown in Figure 2.Table 1
Typical TPMS.
TPMSImplicit surface expressionP surfaceφpr=cosX+cosY+cosZ=0D surfaceφDr=cosXcosYcosZ−sinXsinYsinZ=0G surfaceφGr=sinXcosZ+sinZcosX+sinYcosZ=0I-WP surfaceφI−WPr=2cosXcosY+cosYcosZ+cosZcosX−cos2X+cos2Y+cos2Z=0FRD surfaceφrF−RD=4cosXcosYcosZ−cos2Xcos2Y+cos2Xcos2Z+cos2Ycos2Z=0L surfaceφrL=0.5sin2XcosYsinZ+sin2YcosZsinX+sin2ZcosYsinY+0.5cos2Xcos2Y+cos2Ycos2Z+cos2Zcos2X+0.15=0K surfaceφrK=0.3cosX+cosY+cosZ+0.3cosX+cosY+cosX+cosY+cosX+cosY−0.4cos2X+cos2Y+cos2Z=0G′ surfaceφrG=sin2XcosYsinZ+sin2YcosZsinX+sin2ZcosYsinY+0.32=0P2-DG surfaceφrG=sin2XcosYsinZ+sin2YcosZsinX+sin2ZcosYsinY+0.1cos2Xcos2Y+cos2Ycos2Z+cos2Zcos2X−0.2cos2X+cos2Y+cos2Z=0Figure 2
TPMS. (a)P surface. (b) D surface. (c) G surface. (d) I-WP surface. (e) FRD surface. (f) L surface. (g) K surface. (h) G′ surface. (i) P2-DG surface. Different types of TPMS constructed by marching cubes method.
(a)(b)(c)(d)(e)(f)(g)(h)(i)For the same type of TPMS, the TPMS structure could also be changed by controlling the distancek from the point of φTPMSx,y,z,k to the φTPMSx,y,z [28] as shown in Figure 3.(2)ϕTPMSx,y,z,k=ϕTPMSx,y,z+k=0.Figure 3
FRD unit under differentk values. (a) k = −0.9. (b) k = 0. (c) k = +0.9. The TPMS structure could be changed by controlling the k from the point of φTPMSx,y,z,k to the minimum surface defined by φTPMSx,y,z.
(a)(b)(c)


## 3.2. Fusion Transition of Different Architecture TPMS Units

The transition of different structures of TPMS (including the same type of TPMS unit under differentk values) was performed by the weight fusion method [29].(3)ϕTPMSx,y,z=∑j=1nαjx,y,zϕjx,y,z,where n is the type of TPMS, φjx,y,z is the type of TPMS unit in region j, and αj (x, y, z) is the weight coefficient for region j.(4)αjx,y,z=∑i=1mαijx,y,z=1,x,y,z∈regionj,0x,y,z∉regionj,where αij (x, y, z) is the weight coefficient of the i-th point in region j, i = 1,2, …, m, which could be constructed by the exponential kernel function.(5)αijx,y,z=wij∑i=1mwij,wij=exp−x,y,z−xij,yij,zij2σ2,where (xij, yij, zij) is the coordinate of the i-th point in region j and σ is the width parameter, which controls the local scope of the exponential kernel function.TakeP surface, FRD surface, K surface, and G surface in Table 1 as an example for multi-region fusion transition. Figure 4 shows the fusion transition results.Figure 4
Fusion structure of different types of TPMS. (a) Axonometric drawing. (b) Top view. TakeP surface, FRD surface, K surface, and G surface in Table 1 as an example for multi-region fusion transition.
(a)(b)Figure5 shows the fusion transition results of the same type of TPMS (P surface) under different k.Figure 5
Porous fusion structure of same type of TPMS. (a) Axonometric drawing. (b) Top view. The fusion transition results of the same type of TPMS unit (P unit) under different k.
(a)(b)


## 3.3. Construction of TPMS Point Unit

In this paper, the algebraic method was used to construct TPMS point unit by intersecting the TPMS with the plane perpendicular to the coordinate axis.(6)ϕTPMSx,y,z,k=ϕTPMSx,y,z+k=0,x+=π,x−=−π,y+=π,y−=−π,z=π,z=−π.The FRD,P, and G point units constructed are shown in Figure 6.Figure 6
TPMS point unit. (a) FRD unit,k = 0. (b) P unit, k = 0. (c)G unit, k = 0. Using algebraic method to construct TPMS unit by intersecting the TPMS unit with the plane perpendicular to the coordinate axis.
(a)(b)(c)


## 4. Model Space Subdivision

The purpose of model space subdivision was to provide assembly spatial for the designed TPMS pore unit. In this paper, the method proposed by Zhang et al. [30] was adopted to conduct hexahedral mesh division for the solid model of the defect bone and the irregular hexahedral elements were obtained. The algorithm was shown as follows.Step 1. Octree algorithm was used to establish the topological relationship of the internal mesh.Step 2. The boundary of the model was extracted using dual contouring (DC) algorithm to construct the boundary mesh.Step 3. Relaxation technology was used to improve mesh quality.Taking the triangular prism as an example, the mesh subdivision results are shown in Figure7.Figure 7
Hexahedral mesh division for the triangular prism.In order to control the size and distribution of hexahedral elements, the encryption template proposed by Zhang [31] was used to encrypt hexahedral mesh. The triangular prism mesh encryption results are shown in Figure 8.Figure 8
Encryption of triangular prism mesh. (a) Full encryption. (b) Surface encryption. (c) Line encryption. (d) Point encryption. Using the encryption template proposed by Zhang and Bajaj [31] to encrypt hexahedral mesh.
(a)(b)(c)(d)


## 5. Pore Unit Mapping and Internal Porous Structure Evaluation

### 5.1. Pore Unit Mapping

In order to construct irregular pore unit, the constructed TPMS pore unit was mapped into the hexahedral elements based on the isoparametric element method [32].In this paper, eight-node hexahedral element was used. The mapping relationship of the isoparametric element method was(7)xyz=∑Nixiyizi,where xi, yi, and zi are the coordinate points of hexahedral element nodes in the Cartesian coordinates. x, y, and z are Cartesian coordinates of points inside hexahedral elements. Ni is node-shaped function, which is an interpolation function expressed by natural coordinates. It could be obtained by Lagrange interpolation method:(8)Ni=181+ξξi1+ηηi1+ζζi,i=1,2,…,8,ξ, η, and ζ are the internal points of the entity in the natural coordinate system, −1 < ξ < 1, −1 < η < 1, and −1 < ζ < 1, respectively. ξi, ηi, and ζi are the coordinates of the i-th node in the natural coordinate system. Taking P and FRD pore unit as example, the irregular pore unit obtained after mapping is shown in Figure 9.Figure 9
Mapping relationship. (a)P pore unit, k = 0. (b) FRD pore unit, k = 0.
(a)(b)Namely, the porous bionic bone scaffold model could be constructed based on the Boolean operation.(9)ϕ3D_scaffold=ϕ1∪ϕ2∪…∪ϕi…∪ϕn,where φi is the pore unit of i-th and φ3D_scaffold is the porous bone scaffold with macro- and microstructure.


### 5.2. Evaluation of Internal Porous Structure

The internal porous structure of the bone scaffold plays the important role in the osteogenesis stage of seed cell. Therefore, after constructing the tissue engineering bone scaffold model, it is necessary to evaluate its porosity, pore size, and specific surface area.


#### 5.2.1. Porosity

Porosity is an important parameter to evaluate the pore structure of bone scaffolds.(10)η=1−VPVHEX,where η is the porosity, VHEX is the volume of the hexahedral element, and VP is the volume of the pore unit. In this paper, the pore unit constructed by TPMS was a surface model that is composed of triangular facets, so the center point of the hexahedral element and triangular facet could form a tetrahedron, and the volume of the tetrahedron Vtetrahedronk could be calculated. The volume of the pore unit VP is obtained by accumulating all tetrahedral volumes.(11)VP=∑k=1nVtetrahedronk,Vtetrahedronk=161111∑i=18Nkaixi∑i=18Nkbixi∑i=18Nkcixi18∑i=18xi∑i=18Nkaiyi∑i=18Nkbiyi∑i=18Nkciyi18∑i=18yi∑i=18Nkaizi∑i=18Nkbizi∑i=18Nkcizi18∑i=18zi,where (xi, yi, zi) are the coordinates of the hexahedron nodes of the i-th hexahedral element, and Nkai , Nkbi, and Nkci are the corresponding shape functions of the three vertices of the k-th triangular facet, respectively, which can be obtained by (10). Vtetrahedronk is the volume of k-th tetrahedron.For an arbitrarily irregular hexahedron element, its volume could be calculated by dividing it into five tetrahedrons.(12)VHEX=∑i=aeVi,Va=161111xi5xi1xi6xi8yi5yi1yi6yi8zi5zi1zi6zi8;Vb=161111xi2xi1xi3xi6yi2yi1yi3yi6zi2zi1zi3zi6;Vc=161111xi4xi1xi8xi3yi4yi1yi8yi3zi4zi1zi8zi3Vd=161111xi7xi3xi8xi6yi7yi3yi8yi6zi7zi3zi8zi6;Ve=161111xi1xi3xi8xi6yi1yi3yi8yi6zi1zi3zi8zi6,where xi1,yi1,zi1,xi2,yi2,zi2,xi3,yi3,zi3,xi3,yi3,zi3,xi4,yi4,zi4,xi5,yi5,zi5,xi6,yi6,zi6,xi7,yi7,zi7,xi8,yi8,zi8 are the coordinates of the eight nodes of the i-th hexahedral element. Vi is the volume of tetrahedron.


#### 5.2.2. Pore Size

The pore size is usually defined as the maximum inner sphere diameter of the pore. Therefore, for the pore size calculation, the voxel method [33] was adopted to divide the pore unit equidistant, as shown in Figure 10.Figure 10
Schematic diagram of equidistant plane intersects with the porous unit. (a) Equidistant plane intersects with theP unit, k = 0. (b) Top view. For the pore size calculation, the voxel method [30] was adopted to divide the pore unit equidistant.
(a)(b)The maximum inscribed circle diameterd of the pore unit was approximately solved by the Euclidean distance [34].(13)d=2maxminrxvi,yvi,zvi,xvj,yvj,zvj,rxvi,yvi,zvi,xvj,yvj,zvj,xvi−xvj2+yvi−yvj2+zvi−zvj2,where r is the distance from the center of the i-th voxel marked as “1” to the center of the nearest voxel marked as “0.” xvi,yvi,zvi are the center coordinates of the i-th voxel marked as “1,” and xvi,yvi,zvi are the center coordinates of the j-th voxel marked as “0.”


#### 5.2.3. Specific Surface Area

The specific surface area is a description of the pore space, the adsorption capacity would be stronger, and the space for cell adhesion, migration, and proliferation would be wider with the specific surface area larger.(14)δ=∑i=1nSiVP=∑i=1nSiVHEX⋅η,where Si is the area of the triangle in surface model of the TPMS unit, which could be calculated by Helen’s formula.


## 5.1. Pore Unit Mapping

In order to construct irregular pore unit, the constructed TPMS pore unit was mapped into the hexahedral elements based on the isoparametric element method [32].In this paper, eight-node hexahedral element was used. The mapping relationship of the isoparametric element method was(7)xyz=∑Nixiyizi,where xi, yi, and zi are the coordinate points of hexahedral element nodes in the Cartesian coordinates. x, y, and z are Cartesian coordinates of points inside hexahedral elements. Ni is node-shaped function, which is an interpolation function expressed by natural coordinates. It could be obtained by Lagrange interpolation method:(8)Ni=181+ξξi1+ηηi1+ζζi,i=1,2,…,8,ξ, η, and ζ are the internal points of the entity in the natural coordinate system, −1 < ξ < 1, −1 < η < 1, and −1 < ζ < 1, respectively. ξi, ηi, and ζi are the coordinates of the i-th node in the natural coordinate system. Taking P and FRD pore unit as example, the irregular pore unit obtained after mapping is shown in Figure 9.Figure 9
Mapping relationship. (a)P pore unit, k = 0. (b) FRD pore unit, k = 0.
(a)(b)Namely, the porous bionic bone scaffold model could be constructed based on the Boolean operation.(9)ϕ3D_scaffold=ϕ1∪ϕ2∪…∪ϕi…∪ϕn,where φi is the pore unit of i-th and φ3D_scaffold is the porous bone scaffold with macro- and microstructure.


## 5.2. Evaluation of Internal Porous Structure

The internal porous structure of the bone scaffold plays the important role in the osteogenesis stage of seed cell. Therefore, after constructing the tissue engineering bone scaffold model, it is necessary to evaluate its porosity, pore size, and specific surface area.


### 5.2.1. Porosity

Porosity is an important parameter to evaluate the pore structure of bone scaffolds.(10)η=1−VPVHEX,where η is the porosity, VHEX is the volume of the hexahedral element, and VP is the volume of the pore unit. In this paper, the pore unit constructed by TPMS was a surface model that is composed of triangular facets, so the center point of the hexahedral element and triangular facet could form a tetrahedron, and the volume of the tetrahedron Vtetrahedronk could be calculated. The volume of the pore unit VP is obtained by accumulating all tetrahedral volumes.(11)VP=∑k=1nVtetrahedronk,Vtetrahedronk=161111∑i=18Nkaixi∑i=18Nkbixi∑i=18Nkcixi18∑i=18xi∑i=18Nkaiyi∑i=18Nkbiyi∑i=18Nkciyi18∑i=18yi∑i=18Nkaizi∑i=18Nkbizi∑i=18Nkcizi18∑i=18zi,where (xi, yi, zi) are the coordinates of the hexahedron nodes of the i-th hexahedral element, and Nkai , Nkbi, and Nkci are the corresponding shape functions of the three vertices of the k-th triangular facet, respectively, which can be obtained by (10). Vtetrahedronk is the volume of k-th tetrahedron.For an arbitrarily irregular hexahedron element, its volume could be calculated by dividing it into five tetrahedrons.(12)VHEX=∑i=aeVi,Va=161111xi5xi1xi6xi8yi5yi1yi6yi8zi5zi1zi6zi8;Vb=161111xi2xi1xi3xi6yi2yi1yi3yi6zi2zi1zi3zi6;Vc=161111xi4xi1xi8xi3yi4yi1yi8yi3zi4zi1zi8zi3Vd=161111xi7xi3xi8xi6yi7yi3yi8yi6zi7zi3zi8zi6;Ve=161111xi1xi3xi8xi6yi1yi3yi8yi6zi1zi3zi8zi6,where xi1,yi1,zi1,xi2,yi2,zi2,xi3,yi3,zi3,xi3,yi3,zi3,xi4,yi4,zi4,xi5,yi5,zi5,xi6,yi6,zi6,xi7,yi7,zi7,xi8,yi8,zi8 are the coordinates of the eight nodes of the i-th hexahedral element. Vi is the volume of tetrahedron.


### 5.2.2. Pore Size

The pore size is usually defined as the maximum inner sphere diameter of the pore. Therefore, for the pore size calculation, the voxel method [33] was adopted to divide the pore unit equidistant, as shown in Figure 10.Figure 10
Schematic diagram of equidistant plane intersects with the porous unit. (a) Equidistant plane intersects with theP unit, k = 0. (b) Top view. For the pore size calculation, the voxel method [30] was adopted to divide the pore unit equidistant.
(a)(b)The maximum inscribed circle diameterd of the pore unit was approximately solved by the Euclidean distance [34].(13)d=2maxminrxvi,yvi,zvi,xvj,yvj,zvj,rxvi,yvi,zvi,xvj,yvj,zvj,xvi−xvj2+yvi−yvj2+zvi−zvj2,where r is the distance from the center of the i-th voxel marked as “1” to the center of the nearest voxel marked as “0.” xvi,yvi,zvi are the center coordinates of the i-th voxel marked as “1,” and xvi,yvi,zvi are the center coordinates of the j-th voxel marked as “0.”


### 5.2.3. Specific Surface Area

The specific surface area is a description of the pore space, the adsorption capacity would be stronger, and the space for cell adhesion, migration, and proliferation would be wider with the specific surface area larger.(14)δ=∑i=1nSiVP=∑i=1nSiVHEX⋅η,where Si is the area of the triangle in surface model of the TPMS unit, which could be calculated by Helen’s formula.


## 5.2.1. Porosity

Porosity is an important parameter to evaluate the pore structure of bone scaffolds.(10)η=1−VPVHEX,where η is the porosity, VHEX is the volume of the hexahedral element, and VP is the volume of the pore unit. In this paper, the pore unit constructed by TPMS was a surface model that is composed of triangular facets, so the center point of the hexahedral element and triangular facet could form a tetrahedron, and the volume of the tetrahedron Vtetrahedronk could be calculated. The volume of the pore unit VP is obtained by accumulating all tetrahedral volumes.(11)VP=∑k=1nVtetrahedronk,Vtetrahedronk=161111∑i=18Nkaixi∑i=18Nkbixi∑i=18Nkcixi18∑i=18xi∑i=18Nkaiyi∑i=18Nkbiyi∑i=18Nkciyi18∑i=18yi∑i=18Nkaizi∑i=18Nkbizi∑i=18Nkcizi18∑i=18zi,where (xi, yi, zi) are the coordinates of the hexahedron nodes of the i-th hexahedral element, and Nkai , Nkbi, and Nkci are the corresponding shape functions of the three vertices of the k-th triangular facet, respectively, which can be obtained by (10). Vtetrahedronk is the volume of k-th tetrahedron.For an arbitrarily irregular hexahedron element, its volume could be calculated by dividing it into five tetrahedrons.(12)VHEX=∑i=aeVi,Va=161111xi5xi1xi6xi8yi5yi1yi6yi8zi5zi1zi6zi8;Vb=161111xi2xi1xi3xi6yi2yi1yi3yi6zi2zi1zi3zi6;Vc=161111xi4xi1xi8xi3yi4yi1yi8yi3zi4zi1zi8zi3Vd=161111xi7xi3xi8xi6yi7yi3yi8yi6zi7zi3zi8zi6;Ve=161111xi1xi3xi8xi6yi1yi3yi8yi6zi1zi3zi8zi6,where xi1,yi1,zi1,xi2,yi2,zi2,xi3,yi3,zi3,xi3,yi3,zi3,xi4,yi4,zi4,xi5,yi5,zi5,xi6,yi6,zi6,xi7,yi7,zi7,xi8,yi8,zi8 are the coordinates of the eight nodes of the i-th hexahedral element. Vi is the volume of tetrahedron.


## 5.2.2. Pore Size

The pore size is usually defined as the maximum inner sphere diameter of the pore. Therefore, for the pore size calculation, the voxel method [33] was adopted to divide the pore unit equidistant, as shown in Figure 10.Figure 10
Schematic diagram of equidistant plane intersects with the porous unit. (a) Equidistant plane intersects with theP unit, k = 0. (b) Top view. For the pore size calculation, the voxel method [30] was adopted to divide the pore unit equidistant.
(a)(b)The maximum inscribed circle diameterd of the pore unit was approximately solved by the Euclidean distance [34].(13)d=2maxminrxvi,yvi,zvi,xvj,yvj,zvj,rxvi,yvi,zvi,xvj,yvj,zvj,xvi−xvj2+yvi−yvj2+zvi−zvj2,where r is the distance from the center of the i-th voxel marked as “1” to the center of the nearest voxel marked as “0.” xvi,yvi,zvi are the center coordinates of the i-th voxel marked as “1,” and xvi,yvi,zvi are the center coordinates of the j-th voxel marked as “0.”


## 5.2.3. Specific Surface Area

The specific surface area is a description of the pore space, the adsorption capacity would be stronger, and the space for cell adhesion, migration, and proliferation would be wider with the specific surface area larger.(14)δ=∑i=1nSiVP=∑i=1nSiVHEX⋅η,where Si is the area of the triangle in surface model of the TPMS unit, which could be calculated by Helen’s formula.


## 6. Example Verification

In order to illustrate the effectiveness of the modeling method, a bone scaffold modeling program was developed in Visual Studio 2008 development environment, combined with OpenGL technology. And the cranial prosthesis was used as an example for analysis and discussion, and the modeling process is shown in Figure11.Figure 11
Modeling process of porous bone scaffold. The process includes extracting the surface model, reconstructing the solid model, subdividing the reconstructed model into hexahedrons, obtaining irregular pore unit, and constructing the porous scaffold model with the Boolean operation.Firstly, the surface model of bone defect was extracted from CT or MRT DICOM data by Mimics software. Secondly, the solid model of defect bone was reconstructed by CATIA software. Subsequently, the reconstructed model was subdivided into hexahedrons, and 84 hexahedral elements were generated in this case. Then, the TPMS unit was mapped into hexahedral elements and irregular pore unit could be obtained. Finally, the porous bionic bone scaffold model could be constructed with the Boolean operation.


### 6.1. Design of Single Pore Structure

In this paper, the P and FRD surfaces in Table1 were taken as an example for modeling of porous scaffold. The distance function k is selected as +0.9, +0.7, +0.5, +0.3, +0.1, 0, −0.1, −0.3, −0.5, −0.7, and −0.9, respectively. For the representation of the spatial configuration of the hexahedral elements, the deformation rate and volume were selected to describe it, where the deformation rate was calculated by(15)Q=lminlmax,where lmin is the shortest edge of the hexahedral element and lmax is the longest edge of the hexahedral element.The distribution of the volume and deformation rate of hexahedral elements is shown in Figure12.Figure 12
The distribution of the volume and deformation rate of hexahedral elements. (a) The distribution of the volume. (b) The distribution of the deformation rate.
(a)(b)This paper took theP and FRD unit under k = +0.9, k = 0, and k = −0.9 as an example to introduce, and the constructed porous scaffold model is shown in Figure 13.Figure 13
Design example one: (a) porous scaffold model withP unit, k = +0.9. (b) Porous scaffold model with P unit, k = 0. (c) Porous scaffold model with P unit, k = −0.9. (d) Porous scaffold model with FRD unit, k = +0.9. (e) Porous scaffold model with FRD unit, k = 0. (f) Porous scaffold model with FRD unit, k = −0.9.
(a)(b)(c)(d)(e)(f)


#### 6.1.1. Porosity Distribution

The porosity distribution of P and FRD unit of porous scaffold underk = +0.9, 0, −0.9 is shown in Figure 14.Figure 14
Porosity distribution. It could be seen from Figure14 that the porosity of each pore unit of the porous scaffolds constructed by P and FRD unit showed different distributions. However, for the same type of TPMS, the porosity of each pore unit was a certain value under the same k values.It can be seen from Figure14 that the porosity of each pore unit of the porous scaffolds constructed by P and FRD unit showed different distributions. However, for the same type of TPMS, the porosity of each pore unit was a certain value under the same k values. Under k = −0.9, 0, and +0.9, the porosity of P unit was 79.44%, 56.46%, and 30.95%, and the FRD unit was 33.87%, 55.68%, and 74.84%. Above examples show that the porosity is related to the TPMS pore unit structure (TPMS type and k value), while the spatial configuration of hexahedral elements has no effect on the porosity.To further illustrate the influence ofk value on porosity, the porosity distribution of porous scaffolds under k values of +0.9, +0.7, +0.5, +0.3, +0.1, 0, −0.1, −0.3, −0.5, −0.7, and −0.9 is shown in Figure 15.Figure 15
Porosity distribution under differentk values. For the porous scaffolds constructed by P unit, the porosity decreases with the increase of k value, while for the FRD unit, the porosity increases with the increase of k value.It could be seen from Figure15 that, for the porous scaffolds constructed by P unit, the porosity decreases with the increase of k value, while for the FRD unit, the porosity increases with the increase of k value, and the change of porosity tends to be linear.The fitted linear equation is as follows:(16)ηP=−0.27462∗k+0.55794,ηFRD=0.22618∗k+0.55171,where ηP is the porosity of the pore unit constructed by P unit and ηFRD is the porosity of the pore unit constructed by FRD unit.


#### 6.1.2. Pore Size

The pore size underk = +0.9, 0, and −0.9 is shown in Figure 16.Figure 16
Pore size distribution. (a) Pore size distribution withP unit. (b) Pore size distribution with FRD unit. The pore size of porous scaffold constructed by P and FRD unit showed different distribution rules with different TPMS types and k value.
(a)(b)It could be seen from Figure16 that the pore size of porous scaffold constructed by P and FRD unit showed different distribution rules with a different TPMS type and k value. However, compared with the pore size distribution of the same type of pore unit under different k values, it could be seen that although the pore size distribution varies with the k value, the change trend of the pore size distribution was similar, because the final shape of the pore unit was also controlled by the volume, deformation rate, and spatial distribution of the mapping hexahedral elements.


#### 6.1.3. Specific Surface Area

Figure17 shows the specific surface area distribution of P and FRD unit that under k = +0.9, k = 0, and k = −0.9.Figure 17
Specific surface area distribution. (a) Specific surface area distribution with P unit. (b) Specific surface area distribution with FRD unit. Similar to the porosity distribution, the specific surface area of each pore unit also presents different distribution trends with the different pore unit structures (TPMS type andk values).
(a)(b)It could be seen from Figure17 that, similar to the porosity distribution, the specific surface area of each pore unit also presents different distribution trends with the different pore unit structures (TPMS type and k values). In order to explain the change rule of specific surface area, several pore units were selected for the explanation, as shown in Table 2.Table 2
Structural parameter of pore unit.
The serial number of pore unitVolume of hexahedral (μm3)Deformation rateTPMSkPorosity (%)Surface area (μm2)Specific surface area (μm−1)396.663100.62P−0.979.443.01175.688E − 4816.658100.483.19176.032E − 4315.801100.592.77576.021E − 4191.035110.594.01774.885E − 4535.889100.642.79375.969E − 4645.883100.652.79475.978E − 4396.663100.62FRD+0.974.847.24473.205E − 3816.658100.487.67473.411E − 3315.801100.596.66473.383E − 3191.035110.599.68672.760E − 3535.889100.646.70473.352E − 3645.883100.656.70773.370E − 3It could be seen from Table2 that at 39 and 81 the pore unit structure and the volume of hexahedral elements were similar (6.663E10 μm3 and 6.658E10 μm3), but due to the deformation rate of hexahedral elements (0.62 and 0.48), the specific surface areas of pore unit were different (P unit was 5.688E − 4 μm−1 and 6.032E − 4 μm−1; FRD unit was 3.205E − 3 μm−1 and 3.411E − 3 μm−1), and the specific surface area of the pore unit increases with the decrease of the deformation rate of hexahedral elements, while at 31 and 19 hexahedral elements, the volume deformation rate of hexahedral elements was the same (0.59), but due to the volume of hexahedral elements (5.801E10 μm3 and 1.035E11 μm3), the specific surface area of pore unit was also different (P unit was 6.021E − 4 μm−1 and 4.885E − 4 μm−1; FRD unit was 3.383E − 3 μm−1 and 2.760E − 3 μm−1), and the specific surface area of the pore unit decreases with the increase of the volume of hexahedral elements. However, at 53 and 64 hexahedral elements, due to the similarity of the hexahedral mesh volumes and deformation rates, the specific surface areas of the pore unit were closer (P unit was 5.969E − 4 μm−1 and 5.978E − 4 μm−1; FRD unit was 3.352E − 3 μm−1 and 3.370E − 3 μm−1).In conclusion, the porous scaffold constructed by TPMS unit could not only construct the internal pores with various shapes, but also adjust pore size and shape, porosity, specific surface area, and their spatial distributions through TPMS type,k value, and spatial configuration of hexahedral elements.


### 6.2. Gradient Pore Structure Design Based on Hexahedral Mesh Encryption

In this paper, the gradient structure was designed by changing the spatial configuration of the hexahedral elements. The initial hexahedral mesh was encrypted using fully encrypted template. The hexahedral mesh after encryption is shown in Figure18.Figure 18
Schematic diagram of hexahedral mesh. There are 648 hexahedral elements in the encrypted region, 312 hexahedral elements in the transition region, and 36 hexahedral elements in the sparse region.There are 648 hexahedral elements in the encrypted region, 312 hexahedral elements in the transition region, and 36 hexahedral elements in the sparse region. The volume distribution is shown in Figure19(a), and the deformation rate distribution is shown in Figure 19(b).Figure 19
Distribution of the volume and deformation rate of encrypted hexahedral elements. (a) The distribution of the volume. (b) The distribution of the deformation rate.
(a)(b)The design results are shown in Figure20 with the P unit as the basic pore unit for constructing micro-pores, and k value is 0.Figure 20
Design example two. (a) Gradient porous scaffold withP unit, k = 0. (b), (c) local amplification. The pores of the scaffold are gradient distribution, and the pores in different regions are excessively smooth.It could be seen from Figure20 that the pores of the scaffold are gradient distribution, and the pores in different regions are excessively smooth.


#### 6.2.1. Porosity Distribution

The porosity distribution is shown in Figure21.Figure 21
Porosity distribution. The porosity of the hexahedral elements in transition region, encryption region, and sparse region is all 56.46%.It could be seen from Figure21 that the porosity of the hexahedral elements in transition region, encryption region, and sparse region is all 56.46%. This was because the porosity was determined by the structure of the pore unit, and changing the spatial configuration of the hexahedral elements had no effect on the porosity.


#### 6.2.2. Pore Size and Distribution

The pore size and distribution are shown in Figure22.Figure 22
Pore size distribution. The minimum pore size was reduced by nearly half compared with that before hexahedral mesh encryption, and the distribution of pore size also changed.It could be seen from Figure22 that the minimum pore size was reduced by nearly half compared with that before hexahedral mesh encryption, and the distribution of pore size also changed. This was because the final shape of the pore unit was controlled not only by TPMS type and k value but also by the spatial configuration of the hexahedral elements. The volume, deformation rate, and spatial distribution of hexahedral elements were changed due to the hexahedral mesh encryption, thus indirectly changing the pore size and distribution of the porous scaffold.From Figure23, it could be seen that the specific surface area distribution is more concentrated in the encrypted region and sparse region, while in the transition region the specific surface area distribution is more scattered.Figure 23
The distribution of specific surface area. (a) Specific surface area in encrypted region. (b) Specific surface area in transition region. (c) Specific surface area in sparse region. The specific surface area distribution more concentrated in the encrypted region and sparse region, while in the transition region the specific surface area distribution more scattered.
(a)(b)(c)This is due to the similar spatial configuration of the hexahedral elements in encrypted and sparse region, where the mean volume of the hexahedral elements in encrypted and sparse region was 3.764E9 μm3 and 3.290E10 μm3, and the standard deviation 9.386E8 μm3 and 5.896E9 μm3. The average deformation rate was 0.776 and 0.576, and the standard deviation was 0.067 and 0.063. However, in transition region, the mean volume of the hexahedral elements was 4.828E9 μm3, the standard deviation was 2.604E9 μm3, the mean deformation rate was 0.389, and the standard deviation was 0.131. The spatial configuration of the hexahedral elements in transition region was more variable compared with the encryption and the sparse region, so the specific surface area distribution of the pore unit in this region was more dispersed, mainly concentrated in 1.97E − 3 μm−1 ∼ 7.12E − 3 μm−1, and in encrypted and sparse region was mainly concentrated in 2.5E − 3 μm−1 ∼  3.78E − 3 μm−1 and 9.8E − 4 μm−1 ∼ 1.34E − 3 μm−1.In conclusion, the porous scaffold constructed by TPMS unit, without changing the TPMS unit type andk value, realized the gradient pore structure design only by mesh encryption.


### 6.3. Gradient Pore Structure Design Based on Different k Values

In this section, gradient pore structure design was realized only by changing the distance functionk value. The P unit was still used as the basic pore unit for constructing microscopic pores. k1 of P unit was set as +0.9 in region A, k2 of P unit was set as 0 in region B, and region C was the transition region. The design results are shown in Figure 24.Figure 24
Design example three. (a) Schematic diagram of pore unit distribution. (b) Gradient porous scaffold model. The pore structure of designed scaffold was smooth, and the pore structure present was in a gradient distribution.
(a)(b)It could be seen from Figure24 that the pore structure of designed scaffold was smooth, and the pore structure present was in a gradient distribution. The distribution of the porosity, pore size, and specific surface area are shown in Figure 25.Figure 25
The distribution of porosity, pore size, and specific surface area. (a) The distribution of porosity. (b) The distribution of pore size. (c) The distribution of specific surface area.
(a)(b)(c)It could be seen from the porosity distribution in Figure25(a) that the porosity distribution varies in different regions because the porosity was controlled by the pore unit structure. For the gradient pore structure scaffold constructed by the same type of TPMS, the pore structure in local areas could be controlled by k value. In this case, the porosity in region A and region B was 30.95% and 56.46%, and the transition region C was 43.90%.For the pore size and distribution, it could be seen from Figure25(b) that the pore size and pore size distribution were different compared to the porous scaffold constructed from a single k, due to the different pore structures within each region.Similarly, the specific surface area distribution is shown in Figure25(c), the specific surface area showed different distributions in different regions, and due to the similar spatial configuration of the hexahedral elements, the specific surface area distribution was concentrated in each region. In region A, the specific surface area of each pore unit concentrated 1.38E − 3 μm−1 ∼ 1.82E − 3 μm−1, in region B the specific surface area of each pore unit concentrated 9.8E − 4 μm−1 ∼ 1.34E − 3 μm−1, and in the transition region C the specific surface area of each pore unit concentrated 1.25E − 3 μm−1 ∼ 1.65E − 3 μm−1.Therefore, for the same type of TPMS unit, when constructing gradient pore structure, the porosity, pore size, pore distribution, and specific surface area of bone scaffold could be effectively controlled by adjusting thek value in the local areas.


### 6.4. Gradient Pore Structure Design Based on Different TPMS Types

In this case, the gradient pore structure design was realized by TPMS type without changing the spatial shape of the hexahedron mesh element andk value. FRD and P unit were used as basic pore unit. FRD unit was filled in region A, and k1 is 0. P unit was filled in region B, and the k2 is 0. The design result is shown in Figure 26.Figure 26
Design example four. (a) Gradient porous scaffold model. (b), (c) Local amplification. The pore structure of designed scaffold was smooth, and the pore structure present was in a gradient distribution.It could be seen from Figure26 that the pore structure of designed scaffold was smooth, and the pore structure present was in a gradient distribution. The distribution of the porosity, pore size, and specific surface area is shown in Figure 27.Figure 27
The distribution of porosity, pore size, and specific surface area. (a) The distribution of porosity. (b) The distribution of pore size. (c) The distribution of specific surface area.
(a)(b)(c)The distribution of porosity is shown in Figure27(a). Due to the different types of pore unit in region A and B, the pore structure and porosity were also different. The porosity in region A and B was 55.68% and 56.46%, respectively, and the porosity in transition region C was 67.62% which was higher than the other two regions.As for the pore size and distribution, it could be seen from Figure27(b) that the pore size and pore size distribution were different compared to the porous scaffold constructed from a single type of TPMS, due to the different pore structures within each region.Similarly, the specific surface area distribution is shown in Figure27(c), the specific surface area showed different distributions in different regions, and due to the similar spatial configuration of the hexahedral elements, the specific surface area distribution was concentrated in each region. In region A, the specific surface area of each pore unit concentrated on the range of 1.65E − 3 μm−1 ∼ 2.16E − 3 μm−1, in region B the specific surface area of each pore unit concentrated in the range of 9.8E − 4 μm−1 ∼ 1.34E − 3 μm−1, and in the transition region C the specific surface area of each pore unit concentrated in the range of 1.38E − 3 μm−1 ∼  1.87E − 3 μm−1.It could be seen from the above experiments that the transition of different types of TPMS unit could not only enrich the geometry of pore unit in the scaffold, but also effectively control the porosity, pore size, and distribution of bone scaffolds.


### 6.5. SLM Preparation of Porous Scaffolds

The Ti6Al4V titanium alloy powder produced by EOS company was used as raw material, the porous support designed in example 4 was prepared by dimetal-100 h laser selective melting molding equipment, and the preparation process was carried out under the protection of argon gas. The parameters of the preparation were as follows: scanning speed was 1.2 m/s, the laser power was 280 W, the layer thickness was 35μm, the prepared porous gradient scaffold sample is shown in Figure 28(a), and Figures 28(b)–28(e) show the corresponding SEM photo.Figure 28
Preparation of porous bionic bone scaffold. (a) Porous bionic bone scaffold of example 4. (b), (c), (d), (e) SEM images of example 4. The molding quality of the parts was well, the transition between different types of TPMS unit was relatively smooth, and the prepared porous scaffold was consistent with the design model.It could be seen from Figure28 that the molding quality of the parts was well, the transition between different types of TPMS unit was relatively smooth, and the prepared porous scaffold was consistent with the design model, which verifies the feasibility of the method proposed in this paper.


## 6.1. Design of Single Pore Structure

In this paper, the P and FRD surfaces in Table1 were taken as an example for modeling of porous scaffold. The distance function k is selected as +0.9, +0.7, +0.5, +0.3, +0.1, 0, −0.1, −0.3, −0.5, −0.7, and −0.9, respectively. For the representation of the spatial configuration of the hexahedral elements, the deformation rate and volume were selected to describe it, where the deformation rate was calculated by(15)Q=lminlmax,where lmin is the shortest edge of the hexahedral element and lmax is the longest edge of the hexahedral element.The distribution of the volume and deformation rate of hexahedral elements is shown in Figure12.Figure 12
The distribution of the volume and deformation rate of hexahedral elements. (a) The distribution of the volume. (b) The distribution of the deformation rate.
(a)(b)This paper took theP and FRD unit under k = +0.9, k = 0, and k = −0.9 as an example to introduce, and the constructed porous scaffold model is shown in Figure 13.Figure 13
Design example one: (a) porous scaffold model withP unit, k = +0.9. (b) Porous scaffold model with P unit, k = 0. (c) Porous scaffold model with P unit, k = −0.9. (d) Porous scaffold model with FRD unit, k = +0.9. (e) Porous scaffold model with FRD unit, k = 0. (f) Porous scaffold model with FRD unit, k = −0.9.
(a)(b)(c)(d)(e)(f)


### 6.1.1. Porosity Distribution

The porosity distribution of P and FRD unit of porous scaffold underk = +0.9, 0, −0.9 is shown in Figure 14.Figure 14
Porosity distribution. It could be seen from Figure14 that the porosity of each pore unit of the porous scaffolds constructed by P and FRD unit showed different distributions. However, for the same type of TPMS, the porosity of each pore unit was a certain value under the same k values.It can be seen from Figure14 that the porosity of each pore unit of the porous scaffolds constructed by P and FRD unit showed different distributions. However, for the same type of TPMS, the porosity of each pore unit was a certain value under the same k values. Under k = −0.9, 0, and +0.9, the porosity of P unit was 79.44%, 56.46%, and 30.95%, and the FRD unit was 33.87%, 55.68%, and 74.84%. Above examples show that the porosity is related to the TPMS pore unit structure (TPMS type and k value), while the spatial configuration of hexahedral elements has no effect on the porosity.To further illustrate the influence ofk value on porosity, the porosity distribution of porous scaffolds under k values of +0.9, +0.7, +0.5, +0.3, +0.1, 0, −0.1, −0.3, −0.5, −0.7, and −0.9 is shown in Figure 15.Figure 15
Porosity distribution under differentk values. For the porous scaffolds constructed by P unit, the porosity decreases with the increase of k value, while for the FRD unit, the porosity increases with the increase of k value.It could be seen from Figure15 that, for the porous scaffolds constructed by P unit, the porosity decreases with the increase of k value, while for the FRD unit, the porosity increases with the increase of k value, and the change of porosity tends to be linear.The fitted linear equation is as follows:(16)ηP=−0.27462∗k+0.55794,ηFRD=0.22618∗k+0.55171,where ηP is the porosity of the pore unit constructed by P unit and ηFRD is the porosity of the pore unit constructed by FRD unit.


### 6.1.2. Pore Size

The pore size underk = +0.9, 0, and −0.9 is shown in Figure 16.Figure 16
Pore size distribution. (a) Pore size distribution withP unit. (b) Pore size distribution with FRD unit. The pore size of porous scaffold constructed by P and FRD unit showed different distribution rules with different TPMS types and k value.
(a)(b)It could be seen from Figure16 that the pore size of porous scaffold constructed by P and FRD unit showed different distribution rules with a different TPMS type and k value. However, compared with the pore size distribution of the same type of pore unit under different k values, it could be seen that although the pore size distribution varies with the k value, the change trend of the pore size distribution was similar, because the final shape of the pore unit was also controlled by the volume, deformation rate, and spatial distribution of the mapping hexahedral elements.


### 6.1.3. Specific Surface Area

Figure17 shows the specific surface area distribution of P and FRD unit that under k = +0.9, k = 0, and k = −0.9.Figure 17
Specific surface area distribution. (a) Specific surface area distribution with P unit. (b) Specific surface area distribution with FRD unit. Similar to the porosity distribution, the specific surface area of each pore unit also presents different distribution trends with the different pore unit structures (TPMS type andk values).
(a)(b)It could be seen from Figure17 that, similar to the porosity distribution, the specific surface area of each pore unit also presents different distribution trends with the different pore unit structures (TPMS type and k values). In order to explain the change rule of specific surface area, several pore units were selected for the explanation, as shown in Table 2.Table 2
Structural parameter of pore unit.
The serial number of pore unitVolume of hexahedral (μm3)Deformation rateTPMSkPorosity (%)Surface area (μm2)Specific surface area (μm−1)396.663100.62P−0.979.443.01175.688E − 4816.658100.483.19176.032E − 4315.801100.592.77576.021E − 4191.035110.594.01774.885E − 4535.889100.642.79375.969E − 4645.883100.652.79475.978E − 4396.663100.62FRD+0.974.847.24473.205E − 3816.658100.487.67473.411E − 3315.801100.596.66473.383E − 3191.035110.599.68672.760E − 3535.889100.646.70473.352E − 3645.883100.656.70773.370E − 3It could be seen from Table2 that at 39 and 81 the pore unit structure and the volume of hexahedral elements were similar (6.663E10 μm3 and 6.658E10 μm3), but due to the deformation rate of hexahedral elements (0.62 and 0.48), the specific surface areas of pore unit were different (P unit was 5.688E − 4 μm−1 and 6.032E − 4 μm−1; FRD unit was 3.205E − 3 μm−1 and 3.411E − 3 μm−1), and the specific surface area of the pore unit increases with the decrease of the deformation rate of hexahedral elements, while at 31 and 19 hexahedral elements, the volume deformation rate of hexahedral elements was the same (0.59), but due to the volume of hexahedral elements (5.801E10 μm3 and 1.035E11 μm3), the specific surface area of pore unit was also different (P unit was 6.021E − 4 μm−1 and 4.885E − 4 μm−1; FRD unit was 3.383E − 3 μm−1 and 2.760E − 3 μm−1), and the specific surface area of the pore unit decreases with the increase of the volume of hexahedral elements. However, at 53 and 64 hexahedral elements, due to the similarity of the hexahedral mesh volumes and deformation rates, the specific surface areas of the pore unit were closer (P unit was 5.969E − 4 μm−1 and 5.978E − 4 μm−1; FRD unit was 3.352E − 3 μm−1 and 3.370E − 3 μm−1).In conclusion, the porous scaffold constructed by TPMS unit could not only construct the internal pores with various shapes, but also adjust pore size and shape, porosity, specific surface area, and their spatial distributions through TPMS type,k value, and spatial configuration of hexahedral elements.


## 6.1.1. Porosity Distribution

The porosity distribution of P and FRD unit of porous scaffold underk = +0.9, 0, −0.9 is shown in Figure 14.Figure 14
Porosity distribution. It could be seen from Figure14 that the porosity of each pore unit of the porous scaffolds constructed by P and FRD unit showed different distributions. However, for the same type of TPMS, the porosity of each pore unit was a certain value under the same k values.It can be seen from Figure14 that the porosity of each pore unit of the porous scaffolds constructed by P and FRD unit showed different distributions. However, for the same type of TPMS, the porosity of each pore unit was a certain value under the same k values. Under k = −0.9, 0, and +0.9, the porosity of P unit was 79.44%, 56.46%, and 30.95%, and the FRD unit was 33.87%, 55.68%, and 74.84%. Above examples show that the porosity is related to the TPMS pore unit structure (TPMS type and k value), while the spatial configuration of hexahedral elements has no effect on the porosity.To further illustrate the influence ofk value on porosity, the porosity distribution of porous scaffolds under k values of +0.9, +0.7, +0.5, +0.3, +0.1, 0, −0.1, −0.3, −0.5, −0.7, and −0.9 is shown in Figure 15.Figure 15
Porosity distribution under differentk values. For the porous scaffolds constructed by P unit, the porosity decreases with the increase of k value, while for the FRD unit, the porosity increases with the increase of k value.It could be seen from Figure15 that, for the porous scaffolds constructed by P unit, the porosity decreases with the increase of k value, while for the FRD unit, the porosity increases with the increase of k value, and the change of porosity tends to be linear.The fitted linear equation is as follows:(16)ηP=−0.27462∗k+0.55794,ηFRD=0.22618∗k+0.55171,where ηP is the porosity of the pore unit constructed by P unit and ηFRD is the porosity of the pore unit constructed by FRD unit.


## 6.1.2. Pore Size

The pore size underk = +0.9, 0, and −0.9 is shown in Figure 16.Figure 16
Pore size distribution. (a) Pore size distribution withP unit. (b) Pore size distribution with FRD unit. The pore size of porous scaffold constructed by P and FRD unit showed different distribution rules with different TPMS types and k value.
(a)(b)It could be seen from Figure16 that the pore size of porous scaffold constructed by P and FRD unit showed different distribution rules with a different TPMS type and k value. However, compared with the pore size distribution of the same type of pore unit under different k values, it could be seen that although the pore size distribution varies with the k value, the change trend of the pore size distribution was similar, because the final shape of the pore unit was also controlled by the volume, deformation rate, and spatial distribution of the mapping hexahedral elements.


## 6.1.3. Specific Surface Area

Figure17 shows the specific surface area distribution of P and FRD unit that under k = +0.9, k = 0, and k = −0.9.Figure 17
Specific surface area distribution. (a) Specific surface area distribution with P unit. (b) Specific surface area distribution with FRD unit. Similar to the porosity distribution, the specific surface area of each pore unit also presents different distribution trends with the different pore unit structures (TPMS type andk values).
(a)(b)It could be seen from Figure17 that, similar to the porosity distribution, the specific surface area of each pore unit also presents different distribution trends with the different pore unit structures (TPMS type and k values). In order to explain the change rule of specific surface area, several pore units were selected for the explanation, as shown in Table 2.Table 2
Structural parameter of pore unit.
The serial number of pore unitVolume of hexahedral (μm3)Deformation rateTPMSkPorosity (%)Surface area (μm2)Specific surface area (μm−1)396.663100.62P−0.979.443.01175.688E − 4816.658100.483.19176.032E − 4315.801100.592.77576.021E − 4191.035110.594.01774.885E − 4535.889100.642.79375.969E − 4645.883100.652.79475.978E − 4396.663100.62FRD+0.974.847.24473.205E − 3816.658100.487.67473.411E − 3315.801100.596.66473.383E − 3191.035110.599.68672.760E − 3535.889100.646.70473.352E − 3645.883100.656.70773.370E − 3It could be seen from Table2 that at 39 and 81 the pore unit structure and the volume of hexahedral elements were similar (6.663E10 μm3 and 6.658E10 μm3), but due to the deformation rate of hexahedral elements (0.62 and 0.48), the specific surface areas of pore unit were different (P unit was 5.688E − 4 μm−1 and 6.032E − 4 μm−1; FRD unit was 3.205E − 3 μm−1 and 3.411E − 3 μm−1), and the specific surface area of the pore unit increases with the decrease of the deformation rate of hexahedral elements, while at 31 and 19 hexahedral elements, the volume deformation rate of hexahedral elements was the same (0.59), but due to the volume of hexahedral elements (5.801E10 μm3 and 1.035E11 μm3), the specific surface area of pore unit was also different (P unit was 6.021E − 4 μm−1 and 4.885E − 4 μm−1; FRD unit was 3.383E − 3 μm−1 and 2.760E − 3 μm−1), and the specific surface area of the pore unit decreases with the increase of the volume of hexahedral elements. However, at 53 and 64 hexahedral elements, due to the similarity of the hexahedral mesh volumes and deformation rates, the specific surface areas of the pore unit were closer (P unit was 5.969E − 4 μm−1 and 5.978E − 4 μm−1; FRD unit was 3.352E − 3 μm−1 and 3.370E − 3 μm−1).In conclusion, the porous scaffold constructed by TPMS unit could not only construct the internal pores with various shapes, but also adjust pore size and shape, porosity, specific surface area, and their spatial distributions through TPMS type,k value, and spatial configuration of hexahedral elements.


## 6.2. Gradient Pore Structure Design Based on Hexahedral Mesh Encryption

In this paper, the gradient structure was designed by changing the spatial configuration of the hexahedral elements. The initial hexahedral mesh was encrypted using fully encrypted template. The hexahedral mesh after encryption is shown in Figure18.Figure 18
Schematic diagram of hexahedral mesh. There are 648 hexahedral elements in the encrypted region, 312 hexahedral elements in the transition region, and 36 hexahedral elements in the sparse region.There are 648 hexahedral elements in the encrypted region, 312 hexahedral elements in the transition region, and 36 hexahedral elements in the sparse region. The volume distribution is shown in Figure19(a), and the deformation rate distribution is shown in Figure 19(b).Figure 19
Distribution of the volume and deformation rate of encrypted hexahedral elements. (a) The distribution of the volume. (b) The distribution of the deformation rate.
(a)(b)The design results are shown in Figure20 with the P unit as the basic pore unit for constructing micro-pores, and k value is 0.Figure 20
Design example two. (a) Gradient porous scaffold withP unit, k = 0. (b), (c) local amplification. The pores of the scaffold are gradient distribution, and the pores in different regions are excessively smooth.It could be seen from Figure20 that the pores of the scaffold are gradient distribution, and the pores in different regions are excessively smooth.


### 6.2.1. Porosity Distribution

The porosity distribution is shown in Figure21.Figure 21
Porosity distribution. The porosity of the hexahedral elements in transition region, encryption region, and sparse region is all 56.46%.It could be seen from Figure21 that the porosity of the hexahedral elements in transition region, encryption region, and sparse region is all 56.46%. This was because the porosity was determined by the structure of the pore unit, and changing the spatial configuration of the hexahedral elements had no effect on the porosity.


### 6.2.2. Pore Size and Distribution

The pore size and distribution are shown in Figure22.Figure 22
Pore size distribution. The minimum pore size was reduced by nearly half compared with that before hexahedral mesh encryption, and the distribution of pore size also changed.It could be seen from Figure22 that the minimum pore size was reduced by nearly half compared with that before hexahedral mesh encryption, and the distribution of pore size also changed. This was because the final shape of the pore unit was controlled not only by TPMS type and k value but also by the spatial configuration of the hexahedral elements. The volume, deformation rate, and spatial distribution of hexahedral elements were changed due to the hexahedral mesh encryption, thus indirectly changing the pore size and distribution of the porous scaffold.From Figure23, it could be seen that the specific surface area distribution is more concentrated in the encrypted region and sparse region, while in the transition region the specific surface area distribution is more scattered.Figure 23
The distribution of specific surface area. (a) Specific surface area in encrypted region. (b) Specific surface area in transition region. (c) Specific surface area in sparse region. The specific surface area distribution more concentrated in the encrypted region and sparse region, while in the transition region the specific surface area distribution more scattered.
(a)(b)(c)This is due to the similar spatial configuration of the hexahedral elements in encrypted and sparse region, where the mean volume of the hexahedral elements in encrypted and sparse region was 3.764E9 μm3 and 3.290E10 μm3, and the standard deviation 9.386E8 μm3 and 5.896E9 μm3. The average deformation rate was 0.776 and 0.576, and the standard deviation was 0.067 and 0.063. However, in transition region, the mean volume of the hexahedral elements was 4.828E9 μm3, the standard deviation was 2.604E9 μm3, the mean deformation rate was 0.389, and the standard deviation was 0.131. The spatial configuration of the hexahedral elements in transition region was more variable compared with the encryption and the sparse region, so the specific surface area distribution of the pore unit in this region was more dispersed, mainly concentrated in 1.97E − 3 μm−1 ∼ 7.12E − 3 μm−1, and in encrypted and sparse region was mainly concentrated in 2.5E − 3 μm−1 ∼  3.78E − 3 μm−1 and 9.8E − 4 μm−1 ∼ 1.34E − 3 μm−1.In conclusion, the porous scaffold constructed by TPMS unit, without changing the TPMS unit type andk value, realized the gradient pore structure design only by mesh encryption.


## 6.2.1. Porosity Distribution

The porosity distribution is shown in Figure21.Figure 21
Porosity distribution. The porosity of the hexahedral elements in transition region, encryption region, and sparse region is all 56.46%.It could be seen from Figure21 that the porosity of the hexahedral elements in transition region, encryption region, and sparse region is all 56.46%. This was because the porosity was determined by the structure of the pore unit, and changing the spatial configuration of the hexahedral elements had no effect on the porosity.


## 6.2.2. Pore Size and Distribution

The pore size and distribution are shown in Figure22.Figure 22
Pore size distribution. The minimum pore size was reduced by nearly half compared with that before hexahedral mesh encryption, and the distribution of pore size also changed.It could be seen from Figure22 that the minimum pore size was reduced by nearly half compared with that before hexahedral mesh encryption, and the distribution of pore size also changed. This was because the final shape of the pore unit was controlled not only by TPMS type and k value but also by the spatial configuration of the hexahedral elements. The volume, deformation rate, and spatial distribution of hexahedral elements were changed due to the hexahedral mesh encryption, thus indirectly changing the pore size and distribution of the porous scaffold.From Figure23, it could be seen that the specific surface area distribution is more concentrated in the encrypted region and sparse region, while in the transition region the specific surface area distribution is more scattered.Figure 23
The distribution of specific surface area. (a) Specific surface area in encrypted region. (b) Specific surface area in transition region. (c) Specific surface area in sparse region. The specific surface area distribution more concentrated in the encrypted region and sparse region, while in the transition region the specific surface area distribution more scattered.
(a)(b)(c)This is due to the similar spatial configuration of the hexahedral elements in encrypted and sparse region, where the mean volume of the hexahedral elements in encrypted and sparse region was 3.764E9 μm3 and 3.290E10 μm3, and the standard deviation 9.386E8 μm3 and 5.896E9 μm3. The average deformation rate was 0.776 and 0.576, and the standard deviation was 0.067 and 0.063. However, in transition region, the mean volume of the hexahedral elements was 4.828E9 μm3, the standard deviation was 2.604E9 μm3, the mean deformation rate was 0.389, and the standard deviation was 0.131. The spatial configuration of the hexahedral elements in transition region was more variable compared with the encryption and the sparse region, so the specific surface area distribution of the pore unit in this region was more dispersed, mainly concentrated in 1.97E − 3 μm−1 ∼ 7.12E − 3 μm−1, and in encrypted and sparse region was mainly concentrated in 2.5E − 3 μm−1 ∼  3.78E − 3 μm−1 and 9.8E − 4 μm−1 ∼ 1.34E − 3 μm−1.In conclusion, the porous scaffold constructed by TPMS unit, without changing the TPMS unit type andk value, realized the gradient pore structure design only by mesh encryption.


## 6.3. Gradient Pore Structure Design Based on Different k Values

In this section, gradient pore structure design was realized only by changing the distance functionk value. The P unit was still used as the basic pore unit for constructing microscopic pores. k1 of P unit was set as +0.9 in region A, k2 of P unit was set as 0 in region B, and region C was the transition region. The design results are shown in Figure 24.Figure 24
Design example three. (a) Schematic diagram of pore unit distribution. (b) Gradient porous scaffold model. The pore structure of designed scaffold was smooth, and the pore structure present was in a gradient distribution.
(a)(b)It could be seen from Figure24 that the pore structure of designed scaffold was smooth, and the pore structure present was in a gradient distribution. The distribution of the porosity, pore size, and specific surface area are shown in Figure 25.Figure 25
The distribution of porosity, pore size, and specific surface area. (a) The distribution of porosity. (b) The distribution of pore size. (c) The distribution of specific surface area.
(a)(b)(c)It could be seen from the porosity distribution in Figure25(a) that the porosity distribution varies in different regions because the porosity was controlled by the pore unit structure. For the gradient pore structure scaffold constructed by the same type of TPMS, the pore structure in local areas could be controlled by k value. In this case, the porosity in region A and region B was 30.95% and 56.46%, and the transition region C was 43.90%.For the pore size and distribution, it could be seen from Figure25(b) that the pore size and pore size distribution were different compared to the porous scaffold constructed from a single k, due to the different pore structures within each region.Similarly, the specific surface area distribution is shown in Figure25(c), the specific surface area showed different distributions in different regions, and due to the similar spatial configuration of the hexahedral elements, the specific surface area distribution was concentrated in each region. In region A, the specific surface area of each pore unit concentrated 1.38E − 3 μm−1 ∼ 1.82E − 3 μm−1, in region B the specific surface area of each pore unit concentrated 9.8E − 4 μm−1 ∼ 1.34E − 3 μm−1, and in the transition region C the specific surface area of each pore unit concentrated 1.25E − 3 μm−1 ∼ 1.65E − 3 μm−1.Therefore, for the same type of TPMS unit, when constructing gradient pore structure, the porosity, pore size, pore distribution, and specific surface area of bone scaffold could be effectively controlled by adjusting thek value in the local areas.


## 6.4. Gradient Pore Structure Design Based on Different TPMS Types

In this case, the gradient pore structure design was realized by TPMS type without changing the spatial shape of the hexahedron mesh element andk value. FRD and P unit were used as basic pore unit. FRD unit was filled in region A, and k1 is 0. P unit was filled in region B, and the k2 is 0. The design result is shown in Figure 26.Figure 26
Design example four. (a) Gradient porous scaffold model. (b), (c) Local amplification. The pore structure of designed scaffold was smooth, and the pore structure present was in a gradient distribution.It could be seen from Figure26 that the pore structure of designed scaffold was smooth, and the pore structure present was in a gradient distribution. The distribution of the porosity, pore size, and specific surface area is shown in Figure 27.Figure 27
The distribution of porosity, pore size, and specific surface area. (a) The distribution of porosity. (b) The distribution of pore size. (c) The distribution of specific surface area.
(a)(b)(c)The distribution of porosity is shown in Figure27(a). Due to the different types of pore unit in region A and B, the pore structure and porosity were also different. The porosity in region A and B was 55.68% and 56.46%, respectively, and the porosity in transition region C was 67.62% which was higher than the other two regions.As for the pore size and distribution, it could be seen from Figure27(b) that the pore size and pore size distribution were different compared to the porous scaffold constructed from a single type of TPMS, due to the different pore structures within each region.Similarly, the specific surface area distribution is shown in Figure27(c), the specific surface area showed different distributions in different regions, and due to the similar spatial configuration of the hexahedral elements, the specific surface area distribution was concentrated in each region. In region A, the specific surface area of each pore unit concentrated on the range of 1.65E − 3 μm−1 ∼ 2.16E − 3 μm−1, in region B the specific surface area of each pore unit concentrated in the range of 9.8E − 4 μm−1 ∼ 1.34E − 3 μm−1, and in the transition region C the specific surface area of each pore unit concentrated in the range of 1.38E − 3 μm−1 ∼  1.87E − 3 μm−1.It could be seen from the above experiments that the transition of different types of TPMS unit could not only enrich the geometry of pore unit in the scaffold, but also effectively control the porosity, pore size, and distribution of bone scaffolds.


## 6.5. SLM Preparation of Porous Scaffolds

The Ti6Al4V titanium alloy powder produced by EOS company was used as raw material, the porous support designed in example 4 was prepared by dimetal-100 h laser selective melting molding equipment, and the preparation process was carried out under the protection of argon gas. The parameters of the preparation were as follows: scanning speed was 1.2 m/s, the laser power was 280 W, the layer thickness was 35μm, the prepared porous gradient scaffold sample is shown in Figure 28(a), and Figures 28(b)–28(e) show the corresponding SEM photo.Figure 28
Preparation of porous bionic bone scaffold. (a) Porous bionic bone scaffold of example 4. (b), (c), (d), (e) SEM images of example 4. The molding quality of the parts was well, the transition between different types of TPMS unit was relatively smooth, and the prepared porous scaffold was consistent with the design model.It could be seen from Figure28 that the molding quality of the parts was well, the transition between different types of TPMS unit was relatively smooth, and the prepared porous scaffold was consistent with the design model, which verifies the feasibility of the method proposed in this paper.


## 7. Conclusion

(1)
Taking TPMS as the basic pore unit, the gradient porous scaffold modeling can be constructed by controlling the pore element structure and its spatial shape.(2)
Several examples show that the pore size, porosity, and specific surface area of porous scaffolds could be controlled by adjusting the deformation, distribution, volume of the hexahedral elements, and the structure of the TPMS pore units.(3)
In the next work, the pore types of TPMS unit will be further increased, and the structural model which is functionally equivalent to the real bone structure in terms of pore structure and biological properties will be studied.



---

*Source: 7129482-2022-09-05.xml*