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

Authors: Jaemin Shin; Yongho Choi; Junseok Kim
Journal: Mathematical Problems in Engineering
 (2019)

---

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


---

## Body

## 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 is(1)∂c∂tx,t=∇⋅Mx∇μx,t,x∈Ω,t>0,(2)μ=F′cx,t−ϵ2Δcx,t,where Ω is a domain, c is a mass concentration of a binary mixture, M is a mobility function, Fc=c2c−12/4 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 (Mx,y=0.01+0.99y) on phase separation through a numerical example, which is shown in Figure 1.Figure 1
Temporal evolution (a, b) of morphologies during spinodal phase separation of a ternary system. Reprinted from [12] with permission from Elsevier.
(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 Section2, 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 functionG 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Ωin and its boundary Γ=∂Ωin, we take a rectangular domain Ω=0,Lx×0,Ly embedding Ωin. Let Ωout=Ω/Ωin be outside of Ωin (see Figure 2(a)).Figure 2
(a)Ωin is a target geometry, (b) Ωinh is represented by the closed circles, and (c) boundary control function G is defined as G=1 at closed circles and G=0 at open circles.
(a)
(b)
(c)LetΔx=Lx/Nx and Δy=Ly/Ny be the space step sizes with even integers Nx and Ny, respectively; we consider a uniform mesh h=Δx=Δy. We define a set of cell centers as Ωh=xi,yj:1≤i≤Nx,1≤j≤Ny, where xi=i−0.5h and yj=j−0.5h. Let cijn and μijn be approximations of cxi,yj,tn and μxi,yj,tn, respectively, where tn=nΔt and Δt is a time step size. We define the inner and outer grid domains as Ωinh=Ωh∩Ωin and Ωouth=Ωh∩Ωout, respectively, see Figure 2(b). We denote Γh by the numerical interface, which is a staggered line between Ωinh and Ωouth. At the cell center, the boundary control function G is defined as(3)Gij=1,ifxij∈Ωinh,0,otherwise.At the edge,G is defined as Gi+1/2,j=GijGi+1,j and Gi,j+1/2=GijGi,j+1, 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Γh:(4)ni+1/2,j⋅∇dci+1/2,j=0,ni,j+1/2⋅∇dci,j+1/2=0,where n 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 Γh is defined to be zero, G has the information for the boundary condition. For example, if xi+1/2,j or xi,j+1/2 is the edge point of the boundary, then(5)Gi+1/2,j∇dci+1/2,j=0orGi,j+1/2∇dci,j+1/2=0.Figure3 shows that Ωinh, which is an approximation of the circular disk domain Ωin, converges to Ωin as we refine the mesh size.Figure 3
Convergence ofΩinh to Ωin as we refine the mesh. Mesh sizes are shown below each figure. (a) 16 × 16, (b) 32 × 32, (c) 64 × 64, and (d) 128 × 128.
(a)
(b)
(c)
(d)We define the discrete differentiation ofc as(6)Dxci+1/2,j=ci+1,j−cijh,Dyci,j+1/2=ci,j+1−cijh.The discrete divergence operator is denoted as(7)∇d⋅M∇dcij=Mi+1/2,jci+1,j−cij−Mi−1/2,jcij−ci−1,jh2+Mi,j+1/2ci,j+1−cij−Mi,j−1/2cij−ci,j−1h2.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:(8)cijn+1−cijnΔt=∇d⋅GMn∇dμn+1ij,(9)μijn+1=Fc′cijn+1−ϵ2∇d⋅G∇dcn+1ij+Fc′cijn,where Fcc and Fec are convex functions. The boundary conditions are included in the function G.


### 2.2. Inner Products and Mass Conservation

The discretel2 inner products are defined as(10)ϕ,ψh=h2∑i=1Nx∑j=1NyGijϕijψij,∇dϕ,∇dψe=h2∑i=0Nx∑j=1NyGi+1/2,jDxϕi+1/2,jDxψi+1/2,j+∑i=1Nx∑j=0NyGi,j+1/2Dyϕi,j+1/2Dyψi,j+1/2.Because of the wayG is defined, we can represent the summation as(11)∑i=1Nx∑j=1NyGij=∑Gij=1.Furthermore, the summations in⋅,⋅e can be rewritten in a similar manner:(12)∑i=0Nx∑j=1NyGi+1/2,j=∑Gi+1/2,j=1,∑i=1Nx∑j=0NyGi,j+1/2=∑Gi,j+1/2=1.From the definitions of the inner products and the boundary control functionG, the summation by parts is satisfied:(13)ϕ,∇d⋅G∇dψh=−∇dϕ,∇dψe.The detailed proof for the equality in equation (13) is as follows:(14)ϕ,∇d⋅G∇dψh=∑Gij=1ϕij∇d⋅G∇dψij=∑Gij=1ϕijGi+1/2,jψi+1,j−ψijh2−ϕijGi−1/2,jψij−ψi−1,jh2+∑Gij=1ϕijGi,j+1/2ψi,j+1−ψijh2−ϕijGi,j−1/2ψij−ψi,j−1h2=∑Gi+1/2,j=1ϕijDxψi+1/2,jh−∑Gi−1/2,j=1ϕijDxψi−1/2,jh+∑Gi,j+1/2=1ϕijDyψi,j+1/2h−∑Gi,j−1/2=1ϕijDyψi,j−1/2h=∑Gi+1/2,j=1ϕijDxψi+1/2,jh−∑Gi+1/2,j=1ϕi+1,jDxψi+1/2,jh+∑Gi,j+1/2=1ϕijDyψi,j+1/2h−∑Gi,j+1/2=1ϕi,j+1Dyψi,j+1/2h=−∑Gi+1/2,j=1Dxϕi+1/2,jDxψi+1/2,j−∑Gi,j+1/2=1Dyϕi,j+1/2Dyψi,j+1/2=−∇dϕ,∇dψe.By using the equality in equation (13), the discrete mass conservation is proved by(15)cn+1,1h=cn,1h+Δt∇d⋅GMn∇dμn+1,1h=cn,1h−ΔtGMn∇dμn+1,∇d1e=cn,1h,where 1=1,1,⋯,1.


### 2.3. Energy Dissipation

Next, we prove that the proposed scheme is unconditionally gradient stable. Let us define a discrete energy functional as(16)ℰhc=Fc,1h+ϵ22∇dc,∇dce.It is sufficient to show that the energy dissipation propertyℰhcn+1≤ℰhcn. First, we consider the property of convexity. Because of the convexity of Fcc and Fec, we have(17)Fcn+1−Fcn,1h≤Fc′cn+1−Fe′cn,cn+1−cnh.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 M>0, and then by summation by parts of equations (13) and (9), we have(18)Fc′cn+1−Fe′cn,cn+1−cnh=μn+1+ϵ2∇d⋅G∇dcn+1,cn+1−cnh=Δtμn+1,∇d⋅GMn∇dμn+1h−ϵ2∇dcn+1e2+ϵ2∇dcn+1,∇dcne=−Δt∇dμn+1,Mn∇dμn+1e−ϵ2∇dcn+1e2+ϵ2∇dcn+1,∇dcne≤−ϵ2∇dcn+1e2+ϵ2∇dcn+1,∇dcne.Next, using the inequality in equation (18), we observe that(19)ℰhcn+1−ℰhcn=Fcn+1−Fcn,1h+ϵ22∇dcn+1e2−ϵ22∇dcne2≤Fc′cn+1−Fe′cn,cn+1−cnh+ϵ22∇dcn+1e2−ϵ22∇dcne2=−ϵ22∇dcn+1e2+ϵ2∇dcn+1,∇dcne−ϵ22∇dcne2=−ϵ22∇dcn+1−∇dcne2≤0.Therefore, we can conclude that the schemes of equations (8) and (9) are unconditionally gradient stable.


## 2.1. Boundary Control Function

Given an arbitrary domainΩin and its boundary Γ=∂Ωin, we take a rectangular domain Ω=0,Lx×0,Ly embedding Ωin. Let Ωout=Ω/Ωin be outside of Ωin (see Figure 2(a)).Figure 2
(a)Ωin is a target geometry, (b) Ωinh is represented by the closed circles, and (c) boundary control function G is defined as G=1 at closed circles and G=0 at open circles.
(a)
(b)
(c)LetΔx=Lx/Nx and Δy=Ly/Ny be the space step sizes with even integers Nx and Ny, respectively; we consider a uniform mesh h=Δx=Δy. We define a set of cell centers as Ωh=xi,yj:1≤i≤Nx,1≤j≤Ny, where xi=i−0.5h and yj=j−0.5h. Let cijn and μijn be approximations of cxi,yj,tn and μxi,yj,tn, respectively, where tn=nΔt and Δt is a time step size. We define the inner and outer grid domains as Ωinh=Ωh∩Ωin and Ωouth=Ωh∩Ωout, respectively, see Figure 2(b). We denote Γh by the numerical interface, which is a staggered line between Ωinh and Ωouth. At the cell center, the boundary control function G is defined as(3)Gij=1,ifxij∈Ωinh,0,otherwise.At the edge,G is defined as Gi+1/2,j=GijGi+1,j and Gi,j+1/2=GijGi,j+1, 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Γh:(4)ni+1/2,j⋅∇dci+1/2,j=0,ni,j+1/2⋅∇dci,j+1/2=0,where n 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 Γh is defined to be zero, G has the information for the boundary condition. For example, if xi+1/2,j or xi,j+1/2 is the edge point of the boundary, then(5)Gi+1/2,j∇dci+1/2,j=0orGi,j+1/2∇dci,j+1/2=0.Figure3 shows that Ωinh, which is an approximation of the circular disk domain Ωin, converges to Ωin as we refine the mesh size.Figure 3
Convergence ofΩinh to Ωin as we refine the mesh. Mesh sizes are shown below each figure. (a) 16 × 16, (b) 32 × 32, (c) 64 × 64, and (d) 128 × 128.
(a)
(b)
(c)
(d)We define the discrete differentiation ofc as(6)Dxci+1/2,j=ci+1,j−cijh,Dyci,j+1/2=ci,j+1−cijh.The discrete divergence operator is denoted as(7)∇d⋅M∇dcij=Mi+1/2,jci+1,j−cij−Mi−1/2,jcij−ci−1,jh2+Mi,j+1/2ci,j+1−cij−Mi,j−1/2cij−ci,j−1h2.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:(8)cijn+1−cijnΔt=∇d⋅GMn∇dμn+1ij,(9)μijn+1=Fc′cijn+1−ϵ2∇d⋅G∇dcn+1ij+Fc′cijn,where Fcc and Fec are convex functions. The boundary conditions are included in the function G.


## 2.2. Inner Products and Mass Conservation

The discretel2 inner products are defined as(10)ϕ,ψh=h2∑i=1Nx∑j=1NyGijϕijψij,∇dϕ,∇dψe=h2∑i=0Nx∑j=1NyGi+1/2,jDxϕi+1/2,jDxψi+1/2,j+∑i=1Nx∑j=0NyGi,j+1/2Dyϕi,j+1/2Dyψi,j+1/2.Because of the wayG is defined, we can represent the summation as(11)∑i=1Nx∑j=1NyGij=∑Gij=1.Furthermore, the summations in⋅,⋅e can be rewritten in a similar manner:(12)∑i=0Nx∑j=1NyGi+1/2,j=∑Gi+1/2,j=1,∑i=1Nx∑j=0NyGi,j+1/2=∑Gi,j+1/2=1.From the definitions of the inner products and the boundary control functionG, the summation by parts is satisfied:(13)ϕ,∇d⋅G∇dψh=−∇dϕ,∇dψe.The detailed proof for the equality in equation (13) is as follows:(14)ϕ,∇d⋅G∇dψh=∑Gij=1ϕij∇d⋅G∇dψij=∑Gij=1ϕijGi+1/2,jψi+1,j−ψijh2−ϕijGi−1/2,jψij−ψi−1,jh2+∑Gij=1ϕijGi,j+1/2ψi,j+1−ψijh2−ϕijGi,j−1/2ψij−ψi,j−1h2=∑Gi+1/2,j=1ϕijDxψi+1/2,jh−∑Gi−1/2,j=1ϕijDxψi−1/2,jh+∑Gi,j+1/2=1ϕijDyψi,j+1/2h−∑Gi,j−1/2=1ϕijDyψi,j−1/2h=∑Gi+1/2,j=1ϕijDxψi+1/2,jh−∑Gi+1/2,j=1ϕi+1,jDxψi+1/2,jh+∑Gi,j+1/2=1ϕijDyψi,j+1/2h−∑Gi,j+1/2=1ϕi,j+1Dyψi,j+1/2h=−∑Gi+1/2,j=1Dxϕi+1/2,jDxψi+1/2,j−∑Gi,j+1/2=1Dyϕi,j+1/2Dyψi,j+1/2=−∇dϕ,∇dψe.By using the equality in equation (13), the discrete mass conservation is proved by(15)cn+1,1h=cn,1h+Δt∇d⋅GMn∇dμn+1,1h=cn,1h−ΔtGMn∇dμn+1,∇d1e=cn,1h,where 1=1,1,⋯,1.


## 2.3. Energy Dissipation

Next, we prove that the proposed scheme is unconditionally gradient stable. Let us define a discrete energy functional as(16)ℰhc=Fc,1h+ϵ22∇dc,∇dce.It is sufficient to show that the energy dissipation propertyℰhcn+1≤ℰhcn. First, we consider the property of convexity. Because of the convexity of Fcc and Fec, we have(17)Fcn+1−Fcn,1h≤Fc′cn+1−Fe′cn,cn+1−cnh.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 M>0, and then by summation by parts of equations (13) and (9), we have(18)Fc′cn+1−Fe′cn,cn+1−cnh=μn+1+ϵ2∇d⋅G∇dcn+1,cn+1−cnh=Δtμn+1,∇d⋅GMn∇dμn+1h−ϵ2∇dcn+1e2+ϵ2∇dcn+1,∇dcne=−Δt∇dμn+1,Mn∇dμn+1e−ϵ2∇dcn+1e2+ϵ2∇dcn+1,∇dcne≤−ϵ2∇dcn+1e2+ϵ2∇dcn+1,∇dcne.Next, using the inequality in equation (18), we observe that(19)ℰhcn+1−ℰhcn=Fcn+1−Fcn,1h+ϵ22∇dcn+1e2−ϵ22∇dcne2≤Fc′cn+1−Fe′cn,cn+1−cnh+ϵ22∇dcn+1e2−ϵ22∇dcne2=−ϵ22∇dcn+1e2+ϵ2∇dcn+1,∇dcne−ϵ22∇dcne2=−ϵ22∇dcn+1−∇dcne2≤0.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 Fcc=1/4c2 and Fec=1/4c4−1/2c3, we can rewrite equations (8) and (9) as(20)cijn+1−cijnΔt=∇d⋅GMn∇dμn+1ij,μijn+1=12cijn+1−ϵ2∇d⋅G∇dcn+1ij+cijn3−32cijn2.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 asMc=c1−c. The initial state is taken to be(21)cx,y,0=0.5+0.01⋅randx,y,in a domain Ω=0,1×0,1. The interface Γ is a disk whose radius is 0.45 and center is located at 0.5,0.5. Other numerical parameters are taken as h=1/256, ϵ=0.002, and Δt=0.05h.Figure4 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.Figure 4
Scaled discrete total energyℰhcn/ℰhc0 and average concentration of cx,y,t.


### 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. Figure5 shows the temporal evolution of the spinodal decomposition of a binary mixture with a constant mobility Mx,y=1 and space-dependent mobility Mx,y=0.25x2. In these simulations, the initial condition is(22)cx,y,0=0.5+0.01⋅randx,y,where randx,y is a random number selected from a random distribution between −1 and 1. A 512×128 mesh grid is used on the rectangular domain Ω=0,4×0,1. For the numerical parameters, Δt=10−5 and ϵ=0.0025 are employed. For the phase separation, Figures 5(a) and 5(b) show the solutions at t=0.01 and t=0.5, respectively, with constant mobility M=1. Figure 5(c) shows the different evolution after converting from constant mobility M=1 to space-dependent mobility M=0.25x2 from t=0.01. 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.Figure 5
Phase coarsening and domain growth: (a) phase separation with constant mobility at an early staget=0.01 and snapshots at t=0.5 with (b) constant and (c) space-dependent mobilities.
(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:(23)cx,y,0=0.5+0.01⋅randx,y.The domainΩin is a disk whose radius is 0.45 and center is 0.5,0.5. Variable mobility Mx,y,cn=x2cn1−cn is employed in the computational domain Ω=0,1×0,1. Other parameters are taken as h=1/256, ϵ=0.002, and T=8.Figure6 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.Figure 6
Evolution ofℰhcn/ℰhc0 for different time steps Δt.Figure 7
Numerical solutionsϕx,y,t with different time step Δt at t=8. (a) Δt=1/32, (b) Δt=1/16, (c) Δt=1/8, and (d) Δt=1/4.
(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 0.5,0.5,0.5 on the computational domain Ω=0,13 and Ωin is the inside of Γ. The initial condition is taken to be(24)cx,y,z,0=0.5+0.01⋅randx,y,zon Ωin. The mobility function is defined as M=c1−c32d5+0.01, where dx,y,z=x−0.52+y−0.52+z−0.52. For the computation, we set h=1/128, ϵ=0.005, and Δt=h. 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.Figure 8
Isosurfaces ofc at various times. (a) t=1, (b) t=4, and (c) t=16.
(a)
(b)
(c)Figure9 shows the isosurfaces of internal domains with various radii at t=4. We can show that the phase domain is getting thinner; however, the local average concentration is almost completely preserved.Figure 9
Isosurfaces ofc at t=4 cut with different radii (a) r=0.4, (b) 0.35, and (c) 0.3.
(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.Figure 10
(a) Cylindrical scaffold in [13], (b) a scaffold such as PDLLA, and (c) our design by applying phase separation with space-dependent mobility. Reprinted from [13] with permission from Elsevier.
(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 of(25)ψx,y,z=cosxsiny+cosysinz+coszsinx−0.6+0.05z−4π,within the region of x2+y2≤3.24π2 and 0≤z≤8π.We now consider the design of a specific cylindrical scaffold by applying the variable mobilityM=c1−cMsx with(26)Msx=121+tanh2z−52,which is gradually increasing along the longitudinal direction. We examine the evolution of a random perturbation with a small magnitude about mean composition(27)cx,0=0.5+0.01⋅randx,yon the domain Ω=0,1×0,1×0,2. For the simulation parameters, we take h=1/128, ϵ=0.003, and Δt=h.To highlight the variation of pore size in Figure10(c), we also present cross-sectional images at various heights in Figure 11. In this figure, the red, green, and blue regions indicate c=1, 0.5, and 0, respectively.Figure 11
Cross sections at the indicated heights. (a)z = 1.5 h, (b) z = 15.5 h, (c) z = 63.5 h, and (d) z = 127.5 h.
(a)
(b)
(c)
(d)While Figure11 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).Figure 12
Average concentration ofc⋅,z with respect to the longitudinal direction z.


## 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 asMc=c1−c. The initial state is taken to be(21)cx,y,0=0.5+0.01⋅randx,y,in a domain Ω=0,1×0,1. The interface Γ is a disk whose radius is 0.45 and center is located at 0.5,0.5. Other numerical parameters are taken as h=1/256, ϵ=0.002, and Δt=0.05h.Figure4 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.Figure 4
Scaled discrete total energyℰhcn/ℰhc0 and average concentration of cx,y,t.


## 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. Figure5 shows the temporal evolution of the spinodal decomposition of a binary mixture with a constant mobility Mx,y=1 and space-dependent mobility Mx,y=0.25x2. In these simulations, the initial condition is(22)cx,y,0=0.5+0.01⋅randx,y,where randx,y is a random number selected from a random distribution between −1 and 1. A 512×128 mesh grid is used on the rectangular domain Ω=0,4×0,1. For the numerical parameters, Δt=10−5 and ϵ=0.0025 are employed. For the phase separation, Figures 5(a) and 5(b) show the solutions at t=0.01 and t=0.5, respectively, with constant mobility M=1. Figure 5(c) shows the different evolution after converting from constant mobility M=1 to space-dependent mobility M=0.25x2 from t=0.01. 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.Figure 5
Phase coarsening and domain growth: (a) phase separation with constant mobility at an early staget=0.01 and snapshots at t=0.5 with (b) constant and (c) space-dependent mobilities.
(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:(23)cx,y,0=0.5+0.01⋅randx,y.The domainΩin is a disk whose radius is 0.45 and center is 0.5,0.5. Variable mobility Mx,y,cn=x2cn1−cn is employed in the computational domain Ω=0,1×0,1. Other parameters are taken as h=1/256, ϵ=0.002, and T=8.Figure6 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.Figure 6
Evolution ofℰhcn/ℰhc0 for different time steps Δt.Figure 7
Numerical solutionsϕx,y,t with different time step Δt at t=8. (a) Δt=1/32, (b) Δt=1/16, (c) Δt=1/8, and (d) Δt=1/4.
(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 0.5,0.5,0.5 on the computational domain Ω=0,13 and Ωin is the inside of Γ. The initial condition is taken to be(24)cx,y,z,0=0.5+0.01⋅randx,y,zon Ωin. The mobility function is defined as M=c1−c32d5+0.01, where dx,y,z=x−0.52+y−0.52+z−0.52. For the computation, we set h=1/128, ϵ=0.005, and Δt=h. 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.Figure 8
Isosurfaces ofc at various times. (a) t=1, (b) t=4, and (c) t=16.
(a)
(b)
(c)Figure9 shows the isosurfaces of internal domains with various radii at t=4. We can show that the phase domain is getting thinner; however, the local average concentration is almost completely preserved.Figure 9
Isosurfaces ofc at t=4 cut with different radii (a) r=0.4, (b) 0.35, and (c) 0.3.
(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.Figure 10
(a) Cylindrical scaffold in [13], (b) a scaffold such as PDLLA, and (c) our design by applying phase separation with space-dependent mobility. Reprinted from [13] with permission from Elsevier.
(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 of(25)ψx,y,z=cosxsiny+cosysinz+coszsinx−0.6+0.05z−4π,within the region of x2+y2≤3.24π2 and 0≤z≤8π.We now consider the design of a specific cylindrical scaffold by applying the variable mobilityM=c1−cMsx with(26)Msx=121+tanh2z−52,which is gradually increasing along the longitudinal direction. We examine the evolution of a random perturbation with a small magnitude about mean composition(27)cx,0=0.5+0.01⋅randx,yon the domain Ω=0,1×0,1×0,2. For the simulation parameters, we take h=1/128, ϵ=0.003, and Δt=h.To highlight the variation of pore size in Figure10(c), we also present cross-sectional images at various heights in Figure 11. In this figure, the red, green, and blue regions indicate c=1, 0.5, and 0, respectively.Figure 11
Cross sections at the indicated heights. (a)z = 1.5 h, (b) z = 15.5 h, (c) z = 63.5 h, and (d) z = 127.5 h.
(a)
(b)
(c)
(d)While Figure11 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).Figure 12
Average concentration ofc⋅,z with respect to the longitudinal direction z.


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



---

*Source: 1710270-2019-12-28.xml*