# Finite Element Analysis of Schwarz P Surface Pore Geometries for Tissue-Engineered Scaffolds

Authors: Jaemin Shin; Sungki Kim; Darae Jeong; Hyun Geun Lee; Dongsun Lee; Joong Yeon Lim; Junseok Kim
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

Tissue engineering scaffolds provide temporary mechanical support for tissue regeneration. To regenerate tissues more efficiently, an ideal structure of scaffolds should have appropriate porosity and pore structure. In this paper, we generate the Schwarz primitive (P) surface with various volume fractions using a phase-field model. The phase-field model enables us to design various surface-to-volume ratio structures with high porosity and mechanical properties. Comparing the Schwarz P surface's von Mises stress with that of triply periodic cylinders and cubes, we draw conclusions about the mechanical properties of the Schwarz P surface.


---

## Body

## 1. Introduction

In tissue engineering scaffolds, we should consider the sufficient mechanical strength and stiffness for the scaffold to support to the growing tissue [1]. In addition to mechanical properties, the permeability of the scaffold should be high enough to provide superior diffusion which would facilitate the inflow of nutrients and the disposal of metabolic waste [2]. And the sufficient permeability of scaffold pore architectures increases degree of seeding with cell or growth factors prior to implantation. A large surface area enables cell attachment and growth, whereas a large pore volume is required to accommodate and subsequently deliver a cell mass sufficient for tissue repair. Therefore, it is important for the scaffold to have the suitable pore size, high surface-to-volume ratio, sufficient porosity, and high pore interconnectivity [1].Fabricating solid freeform scaffolds [3, 4] with cubic lattice pores has shown that sharp corners pose machining difficulties. However, fabricating scaffolds with curved pores does not have these manufacturing hurdles. And triply periodic minimal surfaces (TPMS) in general and Schwarz P surface in particular have circular cross-sections throughout their extent making them the ideal pore morphologies for scaffold fabrication [5]. And Schwarz P porous medium has the largest fluid permeability among all of the other triply periodic porous media [6]. From among various pore architectures which have been studied, this study simulated pore microstructure by TPMS for the construction of tissue engineering scaffolds. The objective of this paper is to investigate the mechanical performance of TPMS compared to other configurations.Scaffolds with a high porosity provide more space for cells to move into and begin to thrive. Scaffold porosity greater than 75% is to ensure cell proliferation [7, 8]; however, the maximum porosity is limited by the required mechanical strength of the scaffold. Furthermore, pores that are too small will become clogged and cause a decrease in nutrient diffusion, cellular penetration, and possibly cell death; pores that are too large form gaps that the cells cannot cross. This balance often presents a tradeoff between a denser scaffold providing better function and a more porous scaffold providing better biofactor delivery. This paper reviews what is the best scaffold porosity under mechanical forces. We also accept for the optimal fluid permeability at a porosity volume ratio 1/2 [6] and pore size 1000 micrometers for bone cell proliferation [9–12]. This leads to the conjecture that the maximal stability for a triply periodic porous medium at a porosity volume ratio 1/2 is achieved by the structure that globally minimizes the specific surface.We also note that it is possible that different size pores and/or channels may be needed in different regions of the implant to facilitate the formation of different types of tissue and/or for development of the tissue’s blood supply. But we observe that orthogonally oriented pores with thick walls may be strong, while may not open towards the portion of the host tissue from which new tissue must invade. If that is the case then, the scaffold will act as a barrier to host tissue invasion. Thus, an optimal model of triply periodic scaffold should be chosen which has a smooth surface such as constant mean curvature surfaces.The contents of this paper are as follows. In Section2, we describe how constant mean curvature surface can be obtained through a phase-field algorithm. We introduce a new finite element mesh generation procedure. Geometric models, loading, and materials condition are stated in Section 3. Then Section 4 presents test results and discussion of numerical simulation. Finally, conclusions are given in Section 5.


## 2. Constant Mean Curvature Surfaces with Volume Constraint

### 2.1. The Schwarz P Surface with a Given Volume Fraction

The maximum and minimum values of the normal curvature at a point on a regular surface are called the principal curvaturesκ1 and κ2, respectively (see Figure 1). The mean curvature H is defined as the arithmetic mean of the principal curvatures: H=(κ1+κ2)/2.Figure 1
The maximum and minimum curvatures (κ1 and κ2, resp.) at a point on a surface.Next, we introduce a remarkable class of triply periodic constant mean curvature surfaces (i.e., periodic in three directions). These constant mean curvature surfaces offer great attractions to physical scientists, biologists, and mathematicians. Examples of surfaces are the Schwarz primitive (P), the Schwarz diamond (D), and the Schoen gyroid (G) surfaces as shown in Figure2.Examples of triply periodic constant mean curvature surfaces. (a) Schwarz P, (b) Schwarz D, and (c) Schoen G surfaces.
(a)(b)(c)The triply periodic embedded Schwarz P porous medium is more permeable than other Schwarz models. It also has a smooth surface to open the portion of the host tissue from which new tissue must invade due to its constant mean curvature.To construct the Schwarz P triply periodic constant mean curvature surface with volume constraint, we use a phase-field model. In a phase-field model, the quantityϕ(x,t) is defined to be the mass concentration (with a volume mass) of one of the components. We take the periodic nodal surface (PNS) expansion [13] of the Schwarz P surface as an initial configuration with a desired symmetry, topology, and volume fraction:(2.1)P(x,y,z)=cos⁡2πx+cos⁡2πy+cos⁡2πz+0.5.Now, we describe our algorithm for constructing the Schwarz P surface with a given volume fraction. For clarity of exposition, let us consider a one-dimensional version of the Schwarz P surface and a desired volume fractionα which is greater than 0.5. Please refer to [14] for more details.Step 1.
Initialize a phase-fieldϕ as ϕistep1=cos2πxi+0.5 for i=1,…,N.Step 2.
Truncate nonphysical values such as negative and greater-than-one values:ϕistep2=1 if ϕistep1>1, ϕistep2=0 if ϕistep1<0, and ϕistep2=ϕistep1 otherwise.Step 3.
Relaxϕstep2 by solving the following equation: (ϕistep3-ϕistep2)/Δt = Δd(f(ϕistep3)-(5ϵ)2Δdϕistep3)-(1/4)Δdϕistep2, where the nonlinear function f(ϕ) = ϕ3-(3/2)ϕ2+(3/4)ϕ and we assume that the periodic boundary condition holds for ϕ. The above discrete system is solved by a nonlinear multigrid method. Numerical solution is described in [14].Step 4.
Find a parameterβ which makes a volume fraction of ϕstep3 to α. We define the average volume fraction Vave(ϕstep3,α,β) as Vave(ϕstep3,α,β)=∑i=1Nϕ̃i/N, where ϕ̃i=1 if ϕi>β and ϕ̃i=ϕi otherwise. Taking different β values, we can generate various average volume fractions within a prescribed tolerance. Find an approximate solution to Vave(ϕstep3,α,β)=α by taking the bisection algorithm [15]. Then we define the phase-field ϕstep4 with β: ϕstep4=1 if ϕistep3>β and ϕstep4=ϕistep3 otherwise.Step 5.
Adjust a volume fraction by shiftingϕistep4:ϕistep5=ϕistep4+α-(1/N)∑k=1Nϕkstep4 for i=1,…,N.Step 6.
Relaxϕstep5 by solving the following equation: (ϕistep6-ϕistep5)/Δt = Δd(f(ϕistep6)-ϵ2Δdϕistep6)-(1/4)Δdϕistep5.
Figure3 shows the phase-field ϕ at each step.Phase-fieldϕ at each step.
(a)(b)


### 2.2. Triply Periodic Cylinder and Cube

The volume of the triply periodic cylinder in a unit cubic cell is calculated as follows. LetA, B, and C be the cylinders which are parallel to x, y, and z axis, respectively, and V(A) denote the volume of A. Then,
(2.2)V(A∪B∪C)=V(A)+V(B)+V(C)-V(A∩B)-V(B∩C)-V(C∩A)+V(A∩B∩C),
where A∪B∪C, C,A∩C, and A∩B∩C are shown in Figures 4(a), 4(b), 4(c), and 4(d), respectively. Let a be the radius of the cylinder. Then, volume of each part is given as
(2.3)V(C)=πa2,V(A∩C)=8∫0a∫01-x21-x2dydx=163a3,V(A∩B∩C)=8∫-π/4π/4∫0ar1-r2cos⁡2(θ)drdθ=8(2-2)a3.Figure 4
Triply periodic cylinder: (a)A∪B∪C, (b) C, (c) A∩C, and (d) A∩B∩C.Therefore,V(A∪B∪C)=3πa2-82a3. Similarly, we can find the formula for the triply periodic cube with an edge length b (see Figure 5) as V(A∪B∪C)=3b2-2b3.Triply periodic cube: (a)A∪B∪C and (b) C.
(a)(b)


## 2.1. The Schwarz P Surface with a Given Volume Fraction

The maximum and minimum values of the normal curvature at a point on a regular surface are called the principal curvaturesκ1 and κ2, respectively (see Figure 1). The mean curvature H is defined as the arithmetic mean of the principal curvatures: H=(κ1+κ2)/2.Figure 1
The maximum and minimum curvatures (κ1 and κ2, resp.) at a point on a surface.Next, we introduce a remarkable class of triply periodic constant mean curvature surfaces (i.e., periodic in three directions). These constant mean curvature surfaces offer great attractions to physical scientists, biologists, and mathematicians. Examples of surfaces are the Schwarz primitive (P), the Schwarz diamond (D), and the Schoen gyroid (G) surfaces as shown in Figure2.Examples of triply periodic constant mean curvature surfaces. (a) Schwarz P, (b) Schwarz D, and (c) Schoen G surfaces.
(a)(b)(c)The triply periodic embedded Schwarz P porous medium is more permeable than other Schwarz models. It also has a smooth surface to open the portion of the host tissue from which new tissue must invade due to its constant mean curvature.To construct the Schwarz P triply periodic constant mean curvature surface with volume constraint, we use a phase-field model. In a phase-field model, the quantityϕ(x,t) is defined to be the mass concentration (with a volume mass) of one of the components. We take the periodic nodal surface (PNS) expansion [13] of the Schwarz P surface as an initial configuration with a desired symmetry, topology, and volume fraction:(2.1)P(x,y,z)=cos⁡2πx+cos⁡2πy+cos⁡2πz+0.5.Now, we describe our algorithm for constructing the Schwarz P surface with a given volume fraction. For clarity of exposition, let us consider a one-dimensional version of the Schwarz P surface and a desired volume fractionα which is greater than 0.5. Please refer to [14] for more details.Step 1.
Initialize a phase-fieldϕ as ϕistep1=cos2πxi+0.5 for i=1,…,N.Step 2.
Truncate nonphysical values such as negative and greater-than-one values:ϕistep2=1 if ϕistep1>1, ϕistep2=0 if ϕistep1<0, and ϕistep2=ϕistep1 otherwise.Step 3.
Relaxϕstep2 by solving the following equation: (ϕistep3-ϕistep2)/Δt = Δd(f(ϕistep3)-(5ϵ)2Δdϕistep3)-(1/4)Δdϕistep2, where the nonlinear function f(ϕ) = ϕ3-(3/2)ϕ2+(3/4)ϕ and we assume that the periodic boundary condition holds for ϕ. The above discrete system is solved by a nonlinear multigrid method. Numerical solution is described in [14].Step 4.
Find a parameterβ which makes a volume fraction of ϕstep3 to α. We define the average volume fraction Vave(ϕstep3,α,β) as Vave(ϕstep3,α,β)=∑i=1Nϕ̃i/N, where ϕ̃i=1 if ϕi>β and ϕ̃i=ϕi otherwise. Taking different β values, we can generate various average volume fractions within a prescribed tolerance. Find an approximate solution to Vave(ϕstep3,α,β)=α by taking the bisection algorithm [15]. Then we define the phase-field ϕstep4 with β: ϕstep4=1 if ϕistep3>β and ϕstep4=ϕistep3 otherwise.Step 5.
Adjust a volume fraction by shiftingϕistep4:ϕistep5=ϕistep4+α-(1/N)∑k=1Nϕkstep4 for i=1,…,N.Step 6.
Relaxϕstep5 by solving the following equation: (ϕistep6-ϕistep5)/Δt = Δd(f(ϕistep6)-ϵ2Δdϕistep6)-(1/4)Δdϕistep5.
Figure3 shows the phase-field ϕ at each step.Phase-fieldϕ at each step.
(a)(b)


## 2.2. Triply Periodic Cylinder and Cube

The volume of the triply periodic cylinder in a unit cubic cell is calculated as follows. LetA, B, and C be the cylinders which are parallel to x, y, and z axis, respectively, and V(A) denote the volume of A. Then,
(2.2)V(A∪B∪C)=V(A)+V(B)+V(C)-V(A∩B)-V(B∩C)-V(C∩A)+V(A∩B∩C),
where A∪B∪C, C,A∩C, and A∩B∩C are shown in Figures 4(a), 4(b), 4(c), and 4(d), respectively. Let a be the radius of the cylinder. Then, volume of each part is given as
(2.3)V(C)=πa2,V(A∩C)=8∫0a∫01-x21-x2dydx=163a3,V(A∩B∩C)=8∫-π/4π/4∫0ar1-r2cos⁡2(θ)drdθ=8(2-2)a3.Figure 4
Triply periodic cylinder: (a)A∪B∪C, (b) C, (c) A∩C, and (d) A∩B∩C.Therefore,V(A∪B∪C)=3πa2-82a3. Similarly, we can find the formula for the triply periodic cube with an edge length b (see Figure 5) as V(A∪B∪C)=3b2-2b3.Triply periodic cube: (a)A∪B∪C and (b) C.
(a)(b)


## 3. Numerical Simulation

In this section, we investigate the effect of volume fraction of scaffolds in each shape.


### 3.1. Model and Material

The 3D scaffold models in three shapes such as triply periodic cube, triply periodic cylinder, and Schwarz P were developed by commercially available software ABAQUS 6.8.1 (Dassault Systems, Providence, RI) [16]. And the material of the cube, cylinder, and Schwarz P unit cells are assumed to be homogeneous, isotropic, and linearly elastic poly-DL-lactide (PDLLA). In Table 1, the values of the density (δ), Young’s modulus (E), and Poisson’s ratio (ν) of this material are listed. The values for the material were obtained from the literature [17].Table 1
Values of the density (δ), Young’s modulus (E), and Poisson’s ratio (ν) of the material used in the models.
MaterialDensity (δ)Young’s modulus (E)Poisson’s ratio (ν)PDLLA316 g/cm33.3 GPa0.3The finite element method was used to evaluate and simulate the distribution and magnitude of the von Mises stress on three unit scaffolds under different conditions such as loadings and volume fractions.


### 3.2. Mesh Design

First, we describe a mesh generation procedure for constant mean curvature surfaces in ABAQUS. For the Schwarz P surface model, the outline of mesh generation procedures is as follows.As the phase-fieldϕ is the value from 0 to 1, we divide one-half contour line into regular intervals at some heights. The heights are determined by the circumference of the side and the arc of the cutting line, and connect the neighboring points. The Schwarz P surface consists eight equivalent surfaces, generates the edges of one surface as shown in Figure 6(a), and then makes closed surface by using reflections, see Figure 6(b). Note that the fan-shaped edges in Figure 6(a) are adopted to make a closed surface. Increasing the division points at the circumference and the arc, we can obtain more smooth surface. We use the eight division points for both curves.Generation of the closed surface mesh.
(a)
Edges of one surface(b)
Closed surface with edgesA three-dimensional volume part of the Schwarz P surface was created using the closed edges in GAMBIT 2.2.30 (Ansys Inc., Canonsburg, PA). Using ABAQUS, a volume made by GAMBIT was partitioned in order to reduce the misshapen mesh appearance. Then, the solid model was meshed using hexahedral element to accommodate the irregular geometry features (see Figure7).Mesh generation in ABAQUS.
(a)
Partition(b)
MeshThe other models were partitioned in the same way as the Schwarz P surface model. The meshes were generated with hexahedral elements. A meshed model with a50% volume fraction is shown in Figure 8.The meshed models with hexahedral elements at a50% volume fraction.
(a)
Triply periodic cube(b)
Triply periodic cylinder(c)
Schwarz P


### 3.3. Loading Conditions

A quasistatic compression test was simulated by placing six rigid circle plates on each surface of the unit-cell. As shown in Figure9, each rigid plate moves down slowly by using three concentrated forces (10 N, 30 N, 50 N) to press each surface of the model, respectively. The contact between the rigid plate and the unit cell model was modeled using the penalty model with a friction coefficient of 0.3 due to the rubber-like alginate.Figure 9
Loading condition.


## 3.1. Model and Material

The 3D scaffold models in three shapes such as triply periodic cube, triply periodic cylinder, and Schwarz P were developed by commercially available software ABAQUS 6.8.1 (Dassault Systems, Providence, RI) [16]. And the material of the cube, cylinder, and Schwarz P unit cells are assumed to be homogeneous, isotropic, and linearly elastic poly-DL-lactide (PDLLA). In Table 1, the values of the density (δ), Young’s modulus (E), and Poisson’s ratio (ν) of this material are listed. The values for the material were obtained from the literature [17].Table 1
Values of the density (δ), Young’s modulus (E), and Poisson’s ratio (ν) of the material used in the models.
MaterialDensity (δ)Young’s modulus (E)Poisson’s ratio (ν)PDLLA316 g/cm33.3 GPa0.3The finite element method was used to evaluate and simulate the distribution and magnitude of the von Mises stress on three unit scaffolds under different conditions such as loadings and volume fractions.


## 3.2. Mesh Design

First, we describe a mesh generation procedure for constant mean curvature surfaces in ABAQUS. For the Schwarz P surface model, the outline of mesh generation procedures is as follows.As the phase-fieldϕ is the value from 0 to 1, we divide one-half contour line into regular intervals at some heights. The heights are determined by the circumference of the side and the arc of the cutting line, and connect the neighboring points. The Schwarz P surface consists eight equivalent surfaces, generates the edges of one surface as shown in Figure 6(a), and then makes closed surface by using reflections, see Figure 6(b). Note that the fan-shaped edges in Figure 6(a) are adopted to make a closed surface. Increasing the division points at the circumference and the arc, we can obtain more smooth surface. We use the eight division points for both curves.Generation of the closed surface mesh.
(a)
Edges of one surface(b)
Closed surface with edgesA three-dimensional volume part of the Schwarz P surface was created using the closed edges in GAMBIT 2.2.30 (Ansys Inc., Canonsburg, PA). Using ABAQUS, a volume made by GAMBIT was partitioned in order to reduce the misshapen mesh appearance. Then, the solid model was meshed using hexahedral element to accommodate the irregular geometry features (see Figure7).Mesh generation in ABAQUS.
(a)
Partition(b)
MeshThe other models were partitioned in the same way as the Schwarz P surface model. The meshes were generated with hexahedral elements. A meshed model with a50% volume fraction is shown in Figure 8.The meshed models with hexahedral elements at a50% volume fraction.
(a)
Triply periodic cube(b)
Triply periodic cylinder(c)
Schwarz P


## 3.3. Loading Conditions

A quasistatic compression test was simulated by placing six rigid circle plates on each surface of the unit-cell. As shown in Figure9, each rigid plate moves down slowly by using three concentrated forces (10 N, 30 N, 50 N) to press each surface of the model, respectively. The contact between the rigid plate and the unit cell model was modeled using the penalty model with a friction coefficient of 0.3 due to the rubber-like alginate.Figure 9
Loading condition.


## 4. Results and Discussion

The von Mises stress is defined as(4.1)σe=(σx-σy)2+(σy-σz)2+(σz-σx)2+6(τxy2+τyx2+τzx2)2.Values of von Mises stress equivalent stress at the unit cell scaffold were computed for all variations such as shape of unit cell and different volume fraction.Figures10(a), 10(b), and 10(c) show the von Mises stress for three cases with different volume fractions with the 10 N, 30 N, and 50 N, respectively. The result of Figure 10 indicates that the von Mises stresses on the three models have similar distribution independent of volume fractions and loading forces. Furthermore, the von Mises stresses of the Schwarz P model decreased faster as the volume fraction increased as did the other models.The von Mises stress for three cases with different volume fractions with (a) 10 N, (b) 30 N, and (c) 50 N.
(a)
10 N(b)
30 N(c)
50 NFigure11 shows the magnitude and distribution of the von Mises stress on loading force 50 N for the unit cells with different volume fractions, respectively. It is clearly seen that the stress concentration along sharp edges of the cubic and cylinder unit cell has been significantly reduced in the unit cells of Schwarz P model. The simulation also reveals the stresses on the unit cells of Schwarz P model are more smoothly distributed than those on the other models.Von Mises stress distribution. Volume fractions are (a)30%, (b) 40%, (c) 50%, (d) 60%, (e) 70% with identical loading condition and material property for triply periodic cube (left), triply periodic cylinder (middle), and Schwarz P (right).
(a)(b)(c)(d)(e)It is obvious from these results that in general, as the volume fractions increase for given cells, the von Mises stress also increases. But results imply that Schwarz P scaffold has more stable structures than each of the triply periodic-structures.


## 5. Conclusions

This study focuses on the determination of the relationship between porosity and structural behavior of tissue scaffolds. Finite element analysis of scaffold compression has been conducted to determine the relationship between scaffold porosity and structural behavior under quasistatic loading. The analysis reveals the optimal stress distribution on the triply periodic constant mean surface unit cells. We have evaluated the mechanical property through three different periodic models. The key conclusion we draw in research is that the Schwarz P surface has good stability as well as fluid permeability [6]. These results provide insights into the reason behind the natural choice of constant mean curvature surface forms by biological systems.



---

*Source: 694194-2012-04-30.xml*