Research Article | Open Access

Jaemin Shin, Yongho Choi, Junseok Kim, "The Cahn–Hilliard Equation with Generalized Mobilities in Complex Geometries", *Mathematical Problems in Engineering*, vol. 2019, Article ID 1710270, 10 pages, 2019. https://doi.org/10.1155/2019/1710270

# The Cahn–Hilliard Equation with Generalized Mobilities in Complex Geometries

**Academic Editor:**Chris Goodrich

#### 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.

#### References

- J. W. Cahn and J. E. Hilliard, “Free energy of a nonuniform system. I. Interfacial free energy,”
*The Journal of Chemical Physics*, vol. 28, no. 2, pp. 258–267, 1958. View at: Publisher Site | Google Scholar - J. W. Cahn and J. E. Hilliard, “Spinodal decomposition: a reprise,”
*Acta Metallurgica*, vol. 19, no. 2, pp. 151–161, 1971. View at: Publisher Site | Google Scholar - J. Kim, S. Lee, Y. Choi, S.-M. Lee, and D. Jeong, “Basic principles and practical applications of the Cahn–Hilliard equation,”
*Mathematical Problems in Engineering*, vol. 2016, Article ID 9532608, 11 pages, 2016. View at: Publisher Site | Google Scholar - D. Lee, J.-Y. Huh, D. Jeong, J. Shin, A. Yun, and J. Kim, “Physical, mathematical, and numerical derivations of the Cahn-Hilliard equation,”
*Computational Materials Science*, vol. 81, pp. 216–225, 2014. View at: Publisher Site | Google Scholar - F. Amiri, S. Ziaei-Rad, N. Valizadeh, and T. Rabczuk, “On the use of local maximum entropy approximants for Cahn-Hilliard phase-field models in 2D domains and on surfaces,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 346, pp. 1–24, 2019. View at: Publisher Site | Google Scholar - N. Valizadeh and T. Rabczuk, “Isogeometric analysis for phase-field models of geometric PDEs and high-order PDEs on stationary and evolving surfaces,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 351, pp. 599–642, 2019. View at: Publisher Site | Google Scholar - C. Zimmermann, D. Toshniwal, C. M. Landis, T. J. R. Hughes, K. K. Mandadapu, and R. A. Sauer, “An isogeometric finite element formulation for phase transitions on deforming surfaces,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 351, pp. 441–477, 2019. View at: Publisher Site | Google Scholar - C. M. Elliott and D. A. French, “Numerical studies of the Cahn-Hilliard equation for phase separation,”
*IMA Journal of Applied Mathematics*, vol. 38, no. 2, pp. 97–128, 1987. View at: Publisher Site | Google Scholar - Y. He, Y. Liu, and T. Tang, “On large time-stepping methods for the Cahn–Hilliard equation,”
*Applied Numerical Mathematics*, vol. 57, no. 5–7, pp. 616–628, 2007. View at: Publisher Site | Google Scholar - J. Kim, “A numerical method for the Cahn-Hilliard equation with a variable mobility,”
*Communications in Nonlinear Science and Numerical Simulation*, vol. 12, no. 8, pp. 1560–1571, 2007. View at: Publisher Site | Google Scholar - J. Zhu, L.-Q. Chen, J. Shen, and V. Tikare, “Coarsening kinetics from a variable-mobility Cahn-Hilliard equation: application of a semi-implicit Fourier spectral method,”
*Physical Review E*, vol. 60, no. 4, pp. 3564–3572, 1999. View at: Publisher Site | Google Scholar - J. Kim, “Phase field computations for ternary fluid flows,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 196, no. 45–48, pp. 4779–4788, 2007. View at: Publisher Site | Google Scholar - S. H. Oh, T. H. Kim, and J. H. Lee, “Creating growth factor gradients in three dimensional porous matrix by centrifugation and surface immobilization,”
*Biomaterials*, vol. 32, no. 32, pp. 8254–8260, 2011. View at: Publisher Site | Google Scholar - B. Rentsch, R. Bernhardt, D. Scharnweber, W. Schneiders, S. Rammelt, and C. Rentsch, “Embroidered and surface coated polycaprolactone-co-lactide scaffolds,”
*Biomatter*, vol. 2, no. 3, pp. 158–165, 2012. View at: Publisher Site | Google Scholar - H. Shui, Q. Shi, N. M. Pugno, Q. Chen, and Z. Li, “Effect of mechanical stimulation on the degradation of poly(lactic acid) scaffolds with different designed structures,”
*Journal of the Mechanical Behavior of Biomedical Materials*, vol. 96, no. 96, pp. 324–333, 2019. View at: Publisher Site | Google Scholar - W. Wang, B. Huang, J. J. Byun, and P. Bártolo, “Assessment of PCL/carbon material scaffolds for bone regeneration,”
*Journal of the Mechanical Behavior of Biomedical Materials*, vol. 93, pp. 52–60, 2019. View at: Publisher Site | Google Scholar - Y. Zhang, Y. Xu, Y. Liu et al., “Porous decellularized trachea scaffold prepared by a laser micropore technique,”
*Journal of the Mechanical Behavior of Biomedical Materials*, vol. 90, pp. 96–103, 2019. View at: Publisher Site | Google Scholar - N. Iwasaki, Y. Kasahara, S. Yamane, T. Igarashi, A. Minami, and S.-I. Nisimura, “Chitosan-based hyaluronic acid hybrid polymer fibers as a scaffold biomaterial for cartilage tissue engineering,”
*Polymers*, vol. 3, no. 1, pp. 100–113, 2010. View at: Publisher Site | Google Scholar - S. Kim, S. Ha, S. Lee et al., “Scaffold-type microrobots for targeted cell delivery,” in
*Proceedings of the 2015 12th International Conference on IEEE Ubiquitous Robots and Ambient Intelligence (URAI)*, pp. 526-527, Goyang, Korea, October 2015. View at: Google Scholar - J. Shin, D. Jeong, and J. Kim, “A conservative numerical method for the Cahn-Hilliard equation in complex domains,”
*Journal of Computational Physics*, vol. 230, no. 19, pp. 7441–7455, 2011. View at: Publisher Site | Google Scholar - D. J. Eyre,
*An Unconditionally Stable One-Step Scheme for Gradient Systems*, Unpublished article, 1998. - S. M. Wise, C. Wang, and J. S. Lowengrub, “An energy-stable and convergent finite-difference scheme for the phase field crystal equation,”
*SIAM Journal on Numerical Analysis*, vol. 47, no. 3, pp. 2269–2288, 2009. View at: Publisher Site | Google Scholar - U. Trottenberg, C. W. Oosterlee, and A. Schuller,
*Multigrid*, Academic Press, Cambridge, MA, USA, 2000. - O. A. Abdelaal and S. M. Darwish, “Review of rapid prototyping techniques for tissue engineering scaffolds fabrication,” in
*Characterization and Development of Biosystems and Biomaterials*, Springer, Berlin, Heidelberg, Germany, 2013. View at: Google Scholar - F. P. W. Melchels, K. Bertoldi, R. Gabbrielli, A. H. Velders, J. Feijen, and D. W. Grijpma, “Mathematically defined tissue engineering scaffold architectures prepared by stereolithography,”
*Biomaterials*, vol. 31, no. 27, pp. 6909–6916, 2010. View at: Publisher Site | Google Scholar - D. Yoo, “Heterogeneous minimal surface porous scaffold design using the distance field and radial basis functions,”
*Medical Engineering & Physics*, vol. 34, no. 5, pp. 625–639, 2012. View at: Publisher Site | Google Scholar

#### Copyright

Copyright © 2019 Jaemin Shin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.