Abstract
In this study, we apply a finite difference scheme to solve the Cahn–Hilliard equation with generalized mobilities in complex geometries. This method is conservative and unconditionally gradient stable for all positive variable mobility functions and complex geometries. Herein, we present some numerical experiments to demonstrate the performance of this method. In particular, using the fact that variable mobility changes the growth rate of the phases, we employ space-dependent mobility to design a cylindrical biomedical scaffold with controlled porosity and pore size.
1. Introduction
The Cahn–Hilliard (CH) equation was derived for modeling the phase separation of a binary alloy system [1, 2] and has been utilized in applications in various fields [3]. The CH equation iswhere is a domain, c is a mass concentration of a binary mixture, M is a mobility function, is a free energy function, and ϵ is a positive constant. The zero Neumann or periodic boundary conditions are generally applied to complete the system. The reader can refer to [4] for more details on the CH equation. Furthermore, the CH equation has been widely researched in complex domains and on surfaces with complex geometries and topologies using meshfree methods [5] and isogeometric analysis [6, 7].
Many studies on the CH equation have assumed that the mobility is constant; however, the equation was originally formulated with degenerate mobility [1]. Numerical methods and simulations with variable mobility can be easily found (for example [8–11]). Kim [12] showed the effect of space-dependent mobility () on phase separation through a numerical example, which is shown in Figure 1.
(a)
(b)
Variable mobility changes the rate of phase separation and the coarsening process. This observation motivated us to employ space-dependent mobility to the manufacturing of a biomedical scaffold because many researchers have recently demonstrated that adaptive porosity and pore size are important requirements for a scaffold [13, 14]. A general scaffold has been constructed with a uniform pattern [15–17]. However, various types of biomedical scaffolds have also been proposed to provide reinforcing features [18, 19].
To numerically design a biomedical scaffold, we require an algorithm that can deal with complex geometry. Thus, we consider a simple and efficient method that was proposed in [20]. In this paper, we completely prove the energy stability of the proposed method with variable mobility. Applying space-dependent mobility, we have a concrete result for a cylindrical scaffold with control of both porosity and pore size.
The remainder of this paper is organized as follows. In Section 2, we describe our numerical method for the CH equation with variable mobility in arbitrarily shaped domains. In Section 3, several numerical experiments are reported. Finally, conclusions are presented in Section 4.
2. Numerical Methods
For simplicity, we describe the numerical algorithm in the two-dimensional space, and then its three-dimensional discretization is analogously defined. To solve the given equations in the complex domain, we use the boundary control function G introduced in [20], which we have described below. The summation by parts formula is shown with the proposed inner products. Owing to this formula, we can simply prove the mass conservation and the decreasing of the energy functional.
2.1. Boundary Control Function
Given an arbitrary domain and its boundary , we take a rectangular domain embedding . Let be outside of (see Figure 2(a)).
(a)
(b)
(c)
Let and be the space step sizes with even integers and , respectively; we consider a uniform mesh . We define a set of cell centers as , where and . Let and be approximations of and , respectively, where and is a time step size. We define the inner and outer grid domains as and , respectively, see Figure 2(b). We denote by the numerical interface, which is a staggered line between and . At the cell center, the boundary control function G is defined as
At the edge, G is defined as and , see Figure 2(c). By defining the boundary control function G, we can reuse the multigrid algorithm which is natural to the rectangular domain.
Here, we consider a zero Neumann boundary condition on the staggered boundary :where is the outward normal vector at the cell edge. In this paper, we focus on the method with guaranteed mass conservation and energy stability. Because the value of G at the boundary is defined to be zero, G has the information for the boundary condition. For example, if or is the edge point of the boundary, then
Figure 3 shows that , which is an approximation of the circular disk domain , converges to as we refine the mesh size.
(a)
(b)
(c)
(d)
We define the discrete differentiation of c as
The discrete divergence operator is denoted as
Now, we present the fully discrete scheme for the CH equation with variable mobility in a complex domain, which is based on the convex splitting method [21]. A semi-implicit time and centered difference space discretization of equations (1) and (2) are as follows:where and are convex functions. The boundary conditions are included in the function G.
2.2. Inner Products and Mass Conservation
The discrete inner products are defined as
Because of the way G is defined, we can represent the summation as
Furthermore, the summations in can be rewritten in a similar manner:
From the definitions of the inner products and the boundary control function G, the summation by parts is satisfied:
The detailed proof for the equality in equation (13) is as follows:
By using the equality in equation (13), the discrete mass conservation is proved bywhere .
2.3. Energy Dissipation
Next, we prove that the proposed scheme is unconditionally gradient stable. Let us define a discrete energy functional as
It is sufficient to show that the energy dissipation property . First, we consider the property of convexity. Because of the convexity of and , we have
For a detailed proof of the inequality in equation (17), please refer to [22]. We now proceed to expand the right-hand side of the inequality in equation (17). We assume that , and then by summation by parts of equations (13) and (9), we have
Next, using the inequality in equation (18), we observe that
Therefore, we can conclude that the schemes of equations (8) and (9) are unconditionally gradient stable.
3. Numerical Experiments
First, we report a two-dimensional simulation to numerically demonstrate the mass conservation and energy dissipation. Next, we present the numerical results which highlight different evolutions with constant and space-dependent mobilities in a rectangular domain. Finally, we report the numerical experiments on spinodal decomposition in the spherical and cylindrical domains.
Prior to describing the numerical experiments, the specific numerical method for the convex splitting in equations (8) and (9) is introduced. Applying a possible choice of the linear convex splitting as and , we can rewrite equations (8) and (9) as
The system is solved via the multigrid method [23]. For details on the restriction and prolongation operators in the complex domain with the multigrid method, please refer to [20].
3.1. Mass Conservation and Energy Dissipation
To demonstrate that the numerical scheme inherits the energy decreasing property in a complex domain, we display the evolution of the discrete total energy. The variable mobility is taken as . The initial state is taken to bein a domain . The interface is a disk whose radius is 0.45 and center is located at . Other numerical parameters are taken as , , and .
Figure 4 shows the evolution of the scaled energy functional and average concentration, which implies that the total discrete energy is nonincreasing and the mass is preserved. The insets show the evolution of the concentration field at the times marked by closed circles.
3.2. Rectangular Domain with a Space-Dependent Mobility
We show an example of the phase separation to compare the effect of space-dependent mobility. Figure 5 shows the temporal evolution of the spinodal decomposition of a binary mixture with a constant mobility and space-dependent mobility . In these simulations, the initial condition iswhere is a random number selected from a random distribution between and 1. A mesh grid is used on the rectangular domain . For the numerical parameters, and are employed. For the phase separation, Figures 5(a) and 5(b) show the solutions at and , respectively, with constant mobility . Figure 5(c) shows the different evolution after converting from constant mobility to space-dependent mobility from . For the constant mobility case, uniform coarsening and domain growth are observed. For the space-dependent mobility case, nonuniform coarsening and domain growth are observed, i.e., smaller- and larger-scale coarsening occurs for small and large values of mobility, respectively.
(a)
(b)
(c)
3.3. Stability Test
To demonstrate the energy stability of the proposed scheme, we calculate the numerical solution for large time steps with the initial condition:
The domain is a disk whose radius is 0.45 and center is . Variable mobility is employed in the computational domain . Other parameters are taken as , , and .
Figure 6 shows the evolution of the scaled discrete energy for various time steps. All the curves are nonincreasing and imply that the numerical solutions are energy stable. We note that a time delay occurs if we use a large time step, as shown in Figure 7.
(a)
(b)
(c)
(d)
3.4. Spherical Domain with a Space-Dependent Mobility
We now report a numerical experiment in a spherical domain. The surface is a sphere with radius 0.45 and center on the computational domain and is the inside of . The initial condition is taken to beon . The mobility function is defined as , where . For the computation, we set , , and . Figure 8 shows the isosurfaces of the evolution of c. For visibility, we plot the longitude and latitude of the boundary of the domain . The typical phase separation and coarsening procedure is observable; however, the growth rates of the inner and outer regions are quite different.
(a)
(b)
(c)
Figure 9 shows the isosurfaces of internal domains with various radii at . We can show that the phase domain is getting thinner; however, the local average concentration is almost completely preserved.
(a)
(b)
(c)
3.5. Application to Design the Scaffold
In this section, we consider manufacturing a cylindrically shaped biomedical scaffold for the application of space-dependent mobility. This is motivated by [13] which indicates that a proper spatial gradient of concentration is required for effective regeneration of tissues and organs. Figure 10(a) shows the PCL/F127 cylindrical scaffold, which has a gradient on the concentrations along the longitudinal direction. As shown in Figures 5 and 9, the space-dependent mobility in the CH equation gives the ability to control both the porosity and pore size of the scaffold.
(a)
(b)
(c)
Adaptively changing the pore size is not a new concept in the field of biomedical scaffold manufacturing [24–26]. For example, we can consider the design in Figure 10(b) as the PDLLA in [25]; however, they could not control the porosity of the scaffold. Figure 10(b) is obtained by the isosurface ofwithin the region of and .
We now consider the design of a specific cylindrical scaffold by applying the variable mobility withwhich is gradually increasing along the longitudinal direction. We examine the evolution of a random perturbation with a small magnitude about mean compositionon the domain . For the simulation parameters, we take , , and .
To highlight the variation of pore size in Figure 10(c), we also present cross-sectional images at various heights in Figure 11. In this figure, the red, green, and blue regions indicate , 0.5, and 0, respectively.
(a)
(b)
(c)
(d)
While Figure 11 shows the dynamical variation of the pore size, Figure 12 shows that the porosity of all cross sections is near 50%. Consequently, both the porosity and pore size can be simultaneously controlled by applying phase separation with space-dependent mobility, as shown in Figure 10(c).
4. Conclusions
We proposed a conservative and stable finite difference method to solve the CH equation with variable mobility in complex geometries. We proved that the numerical scheme is conservative and energy dissipative for any time step size. In addition, the CH equation with space-dependent mobility was considered. The two- and three-dimensional numerical results show that local length scale in phase separations can be controlled by using space-dependent mobilities. We expect that space-dependent mobility can be applied to the fabrication of scaffolds to control both local pore size scale and porosity.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
J. Shin was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) and funded by the Korean Government, Ministry of Education (2017R1D1A1B0-3032422). J. S. Kim expresses thanks for the support from the BK21 PLUS program.