Porous Three-Dimensional Scaffold Generation for 3D Printing

: In this paper, we present an efficient numerical method for arbitrary shaped porous structure generation for 3D printing. A phase-field model is employed for modeling phase separation phenomena of diblock copolymers based on the three-dimensional nonlocal Cahn–Hilliard (CH) equation. The nonlocal CH equation is a fourth-order nonlinear partial differential equation. To efficiently solve the governing equation, an unconditionally gradient stable convex splitting method for temporal discretization with a Fourier spectral method for the spatial discretization is adopted. The standard fast Fourier transform is used to speed up the computation. A new local average concentration function is introduced to the original nonlocal CH equation so that we can locally control the morphology of the structure. The proposed algorithm is simple to implement and complex shaped structures can also be implemented with corresponding signed distance fields. Various numerical tests are performed on simple and complex structures. The computational results demonstrate that the proposed method is efficient to generate irregular porous structures for 3D printing.


Introduction
Porosity has many advantages in a variety of fields, including biomaterials, tissue engineering, and clinical medicine. Especially, 3D porous scaffolds have attracted considerable attention in tissue engineering. For instance, bones may be damaged due to accidents or aging in humans; one of the most common methods for severe bone injury is surgical treatment, in which the damaged area is removed and then transplanted for recovery. Because bone is a porous structure, scaffolds with similar porosity distribution and mechanical loading structures need to be studied [1,2].
Porous objects usually have relatively large internal surface areas and are resistant to physical force. Moreover, in the case of a model in which the pores are connected to each other, fluids can flow along them, thereby providing nutrition and cell growth [3]. Therefore, many researches related to generating porous structures have been conducted so far on materials [4], porous sizes [5], architecture [6], topology optimization [7], and 3D printing [8,9] for optimized porous scaffolds. Hollister [10] presented both hierarchical computational topology design and solid free-form fabrication to make 3D anatomic scaffolds with porous structure satisfying the balance between permeability and diffusion. Because creating a porous scaffold that balances mechanical function and material transport is an important issue, how to adjust it efficiently is constantly being studied so far. The authors in [11] used numerical and experimental methods to investigate the influence of the porous structure of triply periodic minimal surface (TPMS) scaffolds on the macroscopic permeability behavior. In addition, advances in 3D printing technology have led to active research into patient-specific structures. As such, the porous scaffolds are closely related to 3D printing. Mohammed et al. [8] examined the creation of unit cells from TPMS renderings as a basis for creating tailored porous structures. In [12], the authors presented a multi-scale porous structure-based lightweight method to make lightweight and strong three-dimensional printed models.
On the one hand, the development of various additive manufacturing technologies [13,14] has made it possible to use various infill technologies for generating 3D objects. These methods are effective for filling the interior of surfaces with porosity. Various methods are used to generate the porosity on surface, such as the Voronoi diagram [15], the Poisson disk sampling [16], and the example-based method [17], etc. Therefore, various combinations of the above methods can be applied to a variety of applications such as designing porosity architecture of bones, for example, by obtaining an internal structure suitable for porous surfaces. On the other hand, there are methods of creating porous objects smoothly using partial differential equations (PDE) to the 3D object itself directly. Creating porous structures using PDE begins with forming patterns on surface. In the previous study, using the reaction-diffusion type PDE, the numerical simulation results with generated patterns on surface were presented in [18]. In this study, the mesh is generated but the interior is not filled, therefore filling inside of object is needed to generate porosity. Then, one can create the object with the similar effect as above. Certain aspects of this approach have been examined implicitly by Jeong et al. [19]. The authors studied numerically a phase-field model for diblock copolymers based on the nonlocal Cahn-Hilliard (CH) equation. The reason why porosity can be obtained is that isosurfaces can be generated using phase-field models. A different average concentration yields a different pattern, i.e., a different porosity by adjusting the phase values.
Therefore, the main purpose of this paper is to use the nonlocal CH equation [19] with a new local average concentration function to generate arbitrary shaped porous structure for 3D printing with space-dependent porous patterns. The CH equation was originally derived from the spinodal decomposition and coarse dynamics of binary alloys. It has been applied to various fields such as image inpainting [20], image segmentation [21], multi-phase flow [22]. Once a signed distance field is given, our method solves the equation and can generate a porous object quickly by using the well-known method, the fast Fourier transform (FFT). Compared with traditional engineering algorithms for generating regular porous architectures, the proposed method allows one to easily adjust parameters to obtain irregular structures with specific volume, density, etc. Regular porous objects are easy to design and fabricate but have shortages such as controlling local pore sizes or distribution of pores. These limitations can be solved through the generation of irregular porous objects [9]. Likewise, these difficulties can be solved by using the modified space-dependent average concentration function in our model.
The paper is organized in the following manner. In Section 2, the governing equation and its computational method are described. In Section 3, the computational experiments are performed. In Section 4, conclusions are given.

Mathematical Model and Numerical Method
Let φ(x, t) be the order parameter. Then, the nonlocal CH equation for modeling a diblock copolymer is given as where Ω is the domain, and α are positive constants, andφ is the average concentration [19].
In this paper, the nonlocal CH equation with a local average concentration functionφ(x) is proposed to generate arbitrary shaped porous structure for 3D printing: With the local average concentration functionφ(x), arbitrary shaped porous structure for 3D printing can be easily generated. When we consider a structure represented by the signed distance function ψ(x, y, z), a local average concentration functionsφ(x) can be defined as follows: where f (x) is a space-dependent function. Figure 1 illustrates a structure in the domain cube. An unconditionally stable Fourier spectral method is used for the nonlocal CH Equation (1) Fourier spectral method has been applied in many studies [23].
k=1 φ s mnk e −i(ξ p x m +η q y n +ν r z k ) be the discrete Fourier transform with a given data {φ s mnr |m = 1, · · · , N x , n = 1, · · · , N y , and k = 1, · · · , N z }, where ξ p = π p/L x , η q = πq/L y , and ν r = πr/L z . The inverse discrete Fourier transform is s pqr e i(ξ p x m +η q y n +ν r z k ) . Let s pqr e i(ξ p x+η q y+ν r z) . Then, By using the linearly stabilized splitting scheme [24] to Equation (1), we have where F (φ) = φ 3 − 3φ. Thus, Equation (5) can be transformed into the discrete Fourier space as follows: Therefore, the following discrete Fourier transform is obtained as follows: Then, using Equation (4), we get the updated numerical solution φ s+1 mnk as s+1 pqr e i(ξ p x m +η q y n +ν r z k ) .
The standard FFT has been used to speed up the computations [25].

Numerical Experiments
In this section, we perform several numerical experiments on various structures and print the resulting 3D models. The interfacial length parameter m is defined by m = mh/[2 √ 2 tanh −1 (0.9)] ≈ 0.24015 mh. This means that there is approximately mh transition layer width when α = 0 [26].

Simple Structures
In this subsection, numerical tests are conducted to generate porous structures inside simple surfaces for which distance functions are well known. The followings are the signed distance functions ψ(x, y, z) for a sphere, a cube, a cylinder, and a torus [27]. The periodic nodal surface (PNS) approximation of the Schwarz primitive (P) surface [28,29] is also listed below.
First, let us consider the case that f (x, y, z) is a constant. The initial phase condition is given as is a random number between −1 and 1. Figure 2a shows the surfaces of the sphere, the cube, the cylinder, the torus, and the P surface with the signed distance functions (7)    Moreover, the following space-dependent average concentration function is used where (x, y, z) are inside the cylinder. Here, the parameters are used as follows: N x = N y = N z = 100, = 2 , ∆t = 10h 2 , and α = 200. Figure 5 shows the phase separation of the cylinder structure over time. Figure 5d is a 3D printed model of the structure in Figure 5c. As post-processing, the above cylinder result is covered on the side with 6h and on the bottom 2h thickness. Figure 6 shows the results in Figure 5c with closed cylinder. Here, the parameters N x = N y = 60, N z = 100, h = 0.01, = 2 , ∆t = 10h 2 , and α = 800 are used. In Figure 7, there are more porous upwards.

Complex Structures
Further examples for complex structures with various average concentration functionsφ(x) are presented. Assume that there is point cloud data of surface for complex structure. First, the Stanford Bunny from the Stanford 3D repository is used as an example for the complex structure. Note that point2trimesh function in MATLAB is used to make a discrete signed distance function for the Stanford bunny. Once the discrete signed distance function is computed, then the procedure is the same as before. Unless otherwise noted, the following initial condition is used for phase separations of complex structure case: Figure 9a-c show the phase separation on the Stanford Bunny at t = 40∆t, 60∆t, and 80∆t, respectively. Here, N x = 150, N y = 150, N z = 130, h = 1/150, ∆t = h 2 , α = 1000, and = 4 are used. Furthermore, the following discrete average concentration function is used: Compared to the constant average function used in the above simulations, the usage of a space-dependent average function produces various patterns in complex structures. For the next step, a similar investigation for the Stanford Armadillo is performed. Figure 10a-c show the phase separation on the Stanford Armadillo at t = 40∆t, 60∆t, and 80∆t, respectively. Here, N x = 146, N y = 160, N z = 139, h = 1/N x , ∆t = h 2 , α = 1000, = 4 , and the same discrete average concentration function Equation (13) are used, where x ∈ (0, 1) × (0, 160h) × (0, 139h). Moreover, 3D printed models of Figures 9c and 10a can be printed as shown in Figure 11.
Next, the following examples show that the proposed method for making porous structures can be useful for bone tissue engineering studies because bone is a porous material (see Figure 12). Let us consider a bone and a spine as shown in Figures 13a and 14a. Similar to the Stanford Bunny, we create a discrete signed distance function using point2trimesh function and use the same initial condition (12). Figure 13b,c illustrate the phase separation on the given bone structure at t = 20∆t and 40∆t, respectively. Here, N x = N y = 86, N z = 2N x , h = 1/N x , ∆t = h 2 , α = 750, and = 4 are used. Moreover, we use the discrete average concentration function (13), where x ∈ (0, 1) × (0, 1) Figure 14b,c illustrate the compact bone made by giving 3h thickness on the given spine, and spongy bone represented by the phase separation on the given spine at t = 40∆t. Figure 14d describes the combination of these two structures as post-processing. Here, N x = N z = 203, N y = 123, h = 1/N x , ∆t = 4h 2 , α = 750, and = 8 are used. Moreover, we use the discrete average concentration function (13), where x ∈ (0, 1) × (0, hN y ) × (0, 1),
(a) P surface (b) D surface (c) G surface A porous scaffold structure can be created by configuring several TPMSs in the shape of cylinder. In [8], authors made 3D porous scaffolds in a tetrahedral, cubic, and cylindrical configuration using Gyroid and P surfaces. They used commercial programs to make Gyroid and Schwarz P unit cells and create scaffold structures. However, the nonlocal CH equation (2) can be used to make similar results. The following parameters on the computational domain Ω = (0, 1) × (0, 1) × (0, 1) are used: N x = N y = N z = 100, = 2 , ∆t = 10h 2 , φ ave = −0.3, and α = 1000. Equation (14) is used as the initial condition in the cylinder (9). Figure 16a-c illustrate the numerical results at 10∆t, 20∆t, and 50∆t, respectively.  Figure 17 shows cylindrical scaffolds using the D surface (15) and G surface (16) with the same parameters as in the previous test. Figures 16 and 17 show the shapes of the scaffolds soften over time.
To smooth the outside of the cylindrical scaffold structure, a post-processing is added after applying Equation (6). Let the whole domain be Ω, the cylindrical domain be Ω out , the slightly smaller cylinder be Ω in as shown in Figure 18.  First, we make a cylinder with the initial condition and a local average concentration, and save the initial data in the domain Ω in . At every iteration, the solution φ s+1 pqr is replaced by the initial value for each (x p , y q , z r ) ∈ Ω in . Figure 19 illustrates the numerical results at t = 50∆t. Both methods based on the simulation results are compared: On the one hand, Equation (6) is solved on the domain Ω. The results are shown in Figure 19a,b; On the other hand, φ s+1 pqr for (x p , y q , z r ) ∈ Ω in is solved and replaced by the initial values. The results are shown in Figure 19c,d. With the second method, we can obtain the cylindrical scaffolds which maintain the original shape of TPMSs inside the cylinder and smooth outside it. Figure 19. Comparison between phase separation of two methods at t = 50∆t: (a,b) Equation (6) is solved on the domain Ω; (c,d) φ s+1 pqr for (x p , y q , z r ) ∈ Ω in is updated with the initial condition.

TPMS Scaffolds Inside Complex Surfaces
Both methods used in the TPMS section above are applied to the Stanford Bunny and Armadillo. Figures 20 and 21 show snapshots of temporal evolutions with average concentration function Equation (13) for the Stanford Bunny and Armadillo, respectively. All the parameters are the same listed in Section 3.2 except for α = 500 and T = 160∆t. Figure 20. (a-c) are snapshots of phase separation on the Stanford Bunny at t = 40∆t, t = 120∆t, and t = 160∆t, respectively. (d-f) are snapshots of phase separation with replacing initial values inside the object at t = 40∆t, t = 120∆t, and t = 160∆t, respectively. Here, the local average concentration function Equation (13) is used. Figure 21. (a-c) are snapshots of phase separation on the Stanford Armadillo at t = 40∆t, 120∆t, and 160∆t, respectively. (d-f) are snapshots of phase separation with replacing initial values inside the object at t = 40∆t, 120∆t, and 160∆t, respectively. Here, the local average concentration function Equation (13) is used.

Conclusions
In this paper, the simple numerical method was proposed to generate irregular porous objects using the phase-field model for diblock copolymers based on the nonlocal CH equation. Using the space-dependent local average concentration function, the computational simulation results of generating various typed porous structures were presented. The most important findings are: • The proposed method can generate porous structures from a solid volume using a distance function. • The porous structure generating algorithm is simple and easy to implement. • The proposed algorithm can control pore shapes by the space-dependent average concentration function.
In this article, we focused on the porous structure generation from solid volumes. In future work, we will investigate the generation of porous structure on arbitrary surfaces in three-dimensional space.