# Basic Principles and Practical Applications of the Cahn–Hilliard Equation

Authors: Junseok Kim; Seunggyu Lee; Yongho Choi; Seok-Min Lee; Darae Jeong
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

The celebrated Cahn–Hilliard (CH) equation was proposed to model the process of phase separation in binary alloys by Cahn and Hilliard. Since then the equation has been extended to a variety of chemical, physical, biological, and other engineering fields such as spinodal decomposition, diblock copolymer, image inpainting, multiphase fluid flows, microstructures with elastic inhomogeneity, tumor growth simulation, and topology optimization. Therefore, it is important to understand the basic mechanism of the CH equation in each modeling type. In this paper, we review the applications of the CH equation and describe the basic mechanism of each modeling type with helpful references and computational simulation results.


---

## Body

## 1. Introduction

The Cahn–Hilliard (CH) equation is a mathematical model of the process of a phase separation in a binary alloy [1, 2]. Its physical applications have been extended to many areas of scientific fields such as spinodal decomposition [1, 3], diblock copolymer [4–6], image inpainting [7], multiphase fluid flows [8–10], microstructures with elastic inhomogeneity [11, 12], tumor growth simulation [13, 14], and topology optimization [15, 16]. In this paper, we review and describe the basic mechanism of each model. The basic building-block equation of the mathematical modeling and dynamical model is the CH equation: (1)∂ϕ∂tx,t=∇·Mϕx,t∇μx,t,x∈Ω,t>0,μx,t=F′ϕx,t-ϵ2Δϕx,t,∂ϕ∂nx,t=∂μ∂nx,t=0,x∈∂Ω,where Ω⊂Rd(d=1,2,3) is a bounded domain with a boundary ∂Ω. The quantity ϕ(x,t)=mα-mβ is a phase-field defined as the difference between the mole fractions of binary mixtures, where mα and mβ are the mole fractions of phases α and β. F(ϕ)=0.25(ϕ2-1)2 is a double-well potential of a homogeneous system of composition ϕ as shown in Figure 1. M(ϕ(x,t)) is a positive mobility, ϵ is a positive constant, and ∂/∂n is the outward normal derivative at the domain boundary. Originally, the concentration-dependent mobility has a form M(ϕ)=M0(1-ϕ)(1+ϕ), where M0 is a constant. This form has the same property in diseases transmission in that transmission happens by contact.Figure 1
Double-well potential functionF(ϕ)=0.25(ϕ2-1)2.The CH equation can be derived from the total free energy functional(2)Eϕ=∫ΩFϕ+ϵ22∇ϕ2dx.Now, we get the chemical potential from (2): μ=δE/δϕ=F′(ϕ)-ϵ2Δϕ. Then, by the mass conservation, we have the CH equation: ϕt=-∇·J, where the flux is defined as J=-M∇μ. The CH system (1) possesses the following two properties: (3)ddt∫Ωϕdx=0,ddtEt≤0. See [18] for more details about the physical, mathematical, and numerical derivations of the binary Cahn–Hilliard equation.The objective of this paper is to review the physical applications of the CH equation. This paper is organized as follows. In Section2, we present spinodal decomposition of a binary alloy, which is the original application of the CH equation. By adding a long-range nonlocal energy, we have a mathematical model for the diblock copolymers in Section 3. A slightly modified CH equation is applied to the inpainting of binary images in Section 4. In Section 5, we present a two-phase immiscible fluid flow model which uses the CH equation to represent two distinct fluids. Elastic strain energy takes an important role in a solid-state phase transformation. When the free energy functional includes this additional effect, the microstructure evolution can be modeled using the CH equation. We consider this in Section 6. In Section 7, we review tumor growth model in which the CH equation accounts for the cell species. In Section 8, we consider structure optimization problem in which we need to minimize mean compliance of the multimaterial structures. Finally, conclusions are drawn in Section 9.


## 2. Spinodal Decomposition

A system of the CH equations (1) is the leading model of spinodal decomposition in binary alloys. Spinodal decomposition is a process by which a mixture of two materials can separate into distinct regions with different material concentrations [1]. Morphological patterns during spinodal decomposition and subsequent coarsening for interface-diffusion-controlled dynamics are shown in Figure 2.Figure 2
Morphological patterns during spinodal decomposition and subsequent coarsening. The upper and lower rows represent the results with (a)ϕ¯=-0.4 and (b) ϕ¯=0, respectively. Times are shown below each figure.
(a) (b)The simulation is performed on a two-dimensional domainΩ=(0,1)×(0,1) with 128×128 mesh grid; that is, the spatial step size h=1/128. The initial condition is given by ϕx,y=ϕ¯+0.05rand⁡x,y, where ϕ¯ is an average value of ϕ(x,y) and rand⁡x,y is a random number function distributed uniformly in [-1,1]. The temporal step size Δt=h2 and ϵ=0.0075 are used. As shown in Figure 2, we can see two different patterns with respect to ϕ¯. The upper and lower rows in Figure 2 represent the time evolution of numerical results with ϕ¯=-0.4 and ϕ¯=0, respectively. In these tests, we used the nonlinearly stabilized splitting scheme [19]. There are many other numerical methods [20–27] that can be applied to the CH equations.


## 3. Diblock Copolymer

A diblock copolymer is a linear-chain molecule consisting of two subchains joined covalently to each other [4]. Let ϕ(x,t)=ρA(x,t)-ρB(x,t) be the order parameter which is defined as the difference between the local volume fractions of A and B monomers at point x and time t. Then, the governing equation by the Ohta–Kawasaki model [28] is as follows: (4)∂ϕx,t∂t=ΔF′ϕx,t-ϵ2Δϕx,t-αϕx,t-ϕ¯,x∈Ω,where Ω⊂Rd(d=1,2,3) is a domain, ϕ¯=∫Ωϕdx/|Ω| is the spatial average of ϕ, and α is inversely proportional to the square of the total chain length of the copolymer [29]. Equation (4) can be derived from the following free energy functional: (5)Eϕ=∫ΩFϕx+ϵ22∇ϕx2dx+α2∬ΩGx,yϕx-ϕ¯ϕy-ϕ¯dxdy, where F(ϕ)=0.25(ϕ2-1)2 and G(x,y) denotes Green’s function of -Δ on Ω with periodic boundary conditions [4].As a test problem, we perform a numerical test with respect toϕ¯ on the domain Ω=(0,L)×(0,L). In this test, we set the initial condition as ϕx,y,0=ϕ¯+0.1rand⁡x,y, where rand(x,y)∈[-1,1] is a randomly generated number. The other parameters used are h=0.03, Δt=0.1, α=100, ϵ=0.035, L=128h, and 128×128 mesh grid. Figures 3(a) and 3(b) show the evolution of the numerical solutions with ϕ¯=-0.3 and ϕ¯=0, respectively. As shown in Figure 3(a), figures in the first row change to the hexagonal geometry as time increases. Unlike this, we can see the lamellar geometry in Figure 3(b) at the equilibrium state.Figure 3
(a) Hexagonal pattern whenϕ¯=-0.3 and (b) lamella pattern when ϕ¯=0 with h=0.03, Δt=0.1, α=100, ϵ=0.035, and 128×128 mesh grid. Times are shown below each column.
(a) (b)


## 4. Image Inpainting

Image inpainting is the filling in of damaged or missing regions of an image with the use of information from surrounding areas [7]. In [7], authors modified the CH equation to achieve fast inpainting of a binary image. Let f(x), where x=(x,y), be a given binary image in a domain Ω and D⊂Ω be the inpainting domain. The image is scaled so that 0≤f≤1. Let c(x,t) be a phase-field which is governed by the following modified CH equation: (6)ctx,t=Δμx,t+λxfx-cx,t,x∈Ω,t>0,μx,t=F′cx,t-ϵ2Δcx,t, where F(c)=0.25c2(1-c)2 and λ(x)=0 if x∈D; otherwise λ(x)=λ0.We perform a test of the inpainting of a double stripe as shown in Figure4(a). Here, gray region in the initial configuration denotes the inpainting region. In this test, we start the numerical test with the following parameters: λ0=10, ϵ=10, h=1/128, and Δt=0.1. And then we switch these values to λ0=0.1, ϵ=0.1, and Δt=1 after 3000 iterations and stop the simulation after 4000 iterations. Figures 4(b), 4(c), and 4(d) are numerical solutions at t=40, 500, and 1300, respectively. The final inpainting result represents the completed stripe. Numerical solution algorithms can be found in [7, 30].Figure 4
(a) Initial data with inpainting region in gray. (b), (c), and (d) are numerical solutions att=40, 500, and 1300, respectively.
(a) (b) (c) (d)


## 5. Two-Phase Fluid Flows

In the two-phase fluid flow problem, we use the CH equation for capturing the interface location between two immiscible fluids. The CH equation provides a good mass conservation property. We model the variable quantities such as viscosity and density by using the phase-field. Also, we model the surface tension effect with the phase-field. The velocity field is governed by the modified Navier–Stokes equation. The phase-field is advected by the bulk velocity. A typical coupled model is(7)ρϕ∂u∂t+u·∇u=-∇p+1ReΔu+1WeSF+ρϕFrg,∇·u=0,∂ϕ∂t+∇·ϕu=1Pe∇·Mϕ∇μ,μ=F′ϕ-CΔϕ,where u is the velocity, p is the pressure, ρ(ϕ) is the density, SF is the surface tension force, g is the gravitational force, ϕ is the phase-field function, M(ϕ) is the mobility function, and μ is the chemical potential. The dimensionless parameters are used as Reynolds number Re=ρ∗L∗V∗/η∗, Weber number We=ρ∗L∗V∗2/σ, Froude number Fr=V∗2(gL∗), Cahn number ϵ2/μ∗, and diffusional Peclet number Pe=L∗V∗/(M∗μ∗), where V∗ is the characteristic velocity, L∗ is the characteristic length, η is the viscosity, σ is the surface tension coefficient, g is the gravitational constant, ϵ is the interfacial parameter, and subscript ∗ denotes the corresponding characteristic values. See [8, 31, 32] for numerical solution algorithms.As a test problem, we consider a drop falling test under the gravity. A schematic diagram of a two-phase domain is shown in Figure5. Figure 6 shows the interfaces of the falling drop at time t=0, 1, 2, and 2.7 from left to right. The simulation is performed on the domain Ω=(0,1)×(0,2) with a 128×256 mesh grid. The temporal and spatial step sizes are Δt=0.0005 and h=1/128, respectively. We use ϵ=0.00151, ρ1/ρ2=3, Re=100, We=12.5, and Fr=1. Also, the radius of the initial circle is set to 0.25.Figure 5
Schematic diagram of a two-phase fluid domain.Figure 6
Evolution of falling drop shapes withρ1/ρ2=3, Re=100, We=12.5, and Fr=1. Times are shown below in each figure.


## 6. Microstructure with Elastic Inhomogeneity

Strain energy stored by the elastic deformation of the object is called the elastic strain energy. This energy performs an important function in determining the transformation path and the corresponding microstructure evolution [12]. Chen and Khachaturyan [33] and Wang et al. [34] have first derived the phase-field model from the microscopic theory of Khachaturyan [35, 36]. The standard phase-field model is combined with the microelasticity theory that allows us to express the elastic energy as a function of composition and order parameters [37]. This model has been applied to a variety of solid-state phase transformations [37, 38] and used to investigate the effect of elastic strain energy on the morphological evolution of the microstructure. For more details, see [12].Now, we derive the phase-field model for microstructural evolution. Let us consider a two-phase microstructure with a composition fieldc(x) and elastic strain ϵijel(x). The total free energy in an inhomogeneous system is defined as (8)F=∫VFc+Eel+κ2∇c2dV.Here, F(c) is the local chemical free energy density defined by a quartic polynomial F(c)=2.5(c-cα)2(c-cβ)2 with two local minima cα and cβ. Assuming that materials are cubic anisotropy, the elastic strain energy density Eel(c,u) is defined as (9)Eelc,u=-C11+C12ecE11+E22-ec+2C44E122+0.5C11E112+E222+C12E11E22,where u=(u,v) is the displacement and C11, C12, and C44 are the cubic elastic parameters and are dependent on the order parameter c. Here, Ck=Ck1c+Ck0(1-c) for k=11, 12, and 44, where Ck0 and Ck1 are elastic constants of matrix and precipitate, respectively. The eigenstrain is e(c)=η(c-cs), which obeys Vegard’s law. Here, cs is the average composition and η is a constant. And E11=ux, E22=vy, and E12=0.5(uy+vx). κ is the gradient energy coefficient which can be related to interatomic interaction parameters. By the variational derivative δF/δc, we obtain the following equation: (10)δFδc=F′c-κΔc+∂Eel∂c.From the total free energy functional (8), we derive the modified CH equation for the two-dimensional microstructure evolution with strong elastic inhomogeneity: (11)∂c∂t=MΔμ,μ=F′c-ϵ2Δc+∂Eel∂c,where M is the mobility. From the definition, we get (12)
∂
E
e
l
∂
c
=
-
C
11
1
-
C
11
0
+
C
12
1
-
C
12
0
e
c
E
11
+
E
22
-
e
c
-
η
C
11
+
C
12
E
11
+
E
22
-
2
e
c
+
2
C
44
1
-
C
44
0
E
12
2
+
0.5
C
11
1
-
C
11
0
E
11
2
+
E
22
2
+
C
12
1
-
C
12
0
E
11
E
22
.Under the condition for plain strain and the quasi-static approximation, we have (13)C11ux+C12vyx+C44uy+C44vxy=C11+C12ecx,C11vy+C12uxy+C44vx+C44uyx=C11+C12ecy,where the subscripts x and y denote partial derivatives with respect to the corresponding arguments.Now, we perform a numerical test withκ=1.2247,cα=0.053,cβ=0.947, the composition expansion coefficient η=0.05, spatial step size h=1, temporal step size Δt=1, and the computational domain Ω=(0,256)×(0,256). With an initial condition, (14)cx,y,0=0.947,ifx-1282+y-1282≤400.060,otherwise,we investigate the effect of elastic inhomogeneity κ. We take the elastic constants C110=232,C120=153, and C440=117 in the same values in [39, 40], which have a cubically anisotropic system. In order to show the effect of elastic inhomogeneity, we keep the same bulk modulus B=C11+C12 and the same ratio of anisotropy δ=2C44/(C11-C12), while changing the ratio of shear modulus κ=C44P/C44M in the matrix and precipitate phases, respectively [39].Figures7(a)–7(f) show the numerical results at time t=1700 with the ratio of shear modulus κ=1.7, 1.35, 1, 0.6, 0.45, and 0.3, respectively. The decrease of κ means that C44P=κC44M; that is, the stiffness of the precipitate is less than the matrix [41–43]. Therefore, the precipitate phase is easy to deform to star shape when κ<1 as shown in Figures 7(e) and 7(f). For a comparison study of the effect of free energies (polynomial and logarithmic functions) on the dynamics, see [44].Figure 7
Morphologies of a circular precipitate with different shear modulusκ at time t=1700.
(a)
κ = 1.7 (b)
κ = 1.35 (c)
κ = 1 (d)
κ = 0.6 (e)
κ = 0.45 (f)
κ = 0.3


## 7. Tumor Growth Simulation

To provide optimal strategies for treatments, a mathematical modeling is very useful since it gives systematic investigation. LetΩ be a computational domain. Let ΩH, ΩV, and ΩD be the healthy, viable, and dead tumor tissues, respectively. We also define ΩT=ΩV∪ΩD as the tumor cell. See Figure 8 for the schematic of these definitions.Figure 8
Schematic of the domain for tumor simulation.Ω, ΩH, ΩV, ΩD, and ΩT are the computational domain, healthy cell, viable cell, dead cell, and tumor, respectively.LetϕT and ϕD be the volume fractions of tumor and dead tumor tissues. Then, the nondimensional equations for the volume fraction of tumor cells are (15)∂ϕT∂t=M∇·ϕT∇μ+ST-∇·ϕTuS,μ=F′ϕT-ϵ2ΔϕT,ST=nGϕTϕT-ϕD-λLϕD,where G(ϕT) is a continuous cut-off function. The nondimensional equations for the volume fraction of necrotic cells are (16)∂ϕD∂t=M∇·ϕD∇μ+SD-∇·ϕDuS,SD=λA+λNHnN-nϕT-ϕD-λLϕD, where H is the Heaviside function. The nondimensional cell velocity is (17)u=-κϕT,ϕD∇p-γϵμ∇ϕT, where γ is a nondimensional measure of the adhesion force. The velocity satisfies (18)∇·u=ST. The nondimensional quasi-steady nutrient equation is (19)0=∇·DϕT∇n+TcϕT,n-nϕT-ϕD, where the diffusion coefficient and nutrient capillary source term are (20)DϕT=DH1-QϕT+QϕT,TCϕT,n=νPH1-QϕT+νPTQϕTnC-n,where DH is the nutrient diffusivity in the healthy tissue and Q(ϕT) is an interpolation function. For more details about the functions and parameters, see [17]. Also, numerical simulation using phase-field method can be found in [45–47].Figures9(a)–9(c) show the numerical results at times t=0, t=50, and t=100, respectively [17]. The parameters are νPH=0.5, nC=1, DH=1, λA=0, λL=1, M=10, νPT=0, κ=1, nN=0.4. λN=3, γ=0, and ϵ=0.1.Figure 9
Evolution of theϕT=0.5 isosurface during the growth of an asymmetrical 3D tumor (reprinted from Wise et al. [17], with permission from Elsevier).
(a)
t = 0 (b)
t = 50 (c)
t = 100


## 8. Topology Optimization

Zhou and Wang [16] extended the phase-field approach to the problem of minimizing the mean compliance of a multimaterial structure. For a four-phase system, we have (21)∂c∂t=∇·Mc∇μ,μ=F′c-ΓϵΔc+wc, where (22)
c
=
α
β
γ
,
μ
=
ξ
ζ
ψ
,
Γ
ϵ
=
2
ϵ
2
ϵ
2
ϵ
2
ϵ
2
2
ϵ
2
ϵ
2
ϵ
2
ϵ
2
2
ϵ
2
,
f
c
=
∂
F
∂
α
∂
F
∂
β
∂
F
∂
γ
,
w
c
=
∂
W
∗
∂
α
∂
W
∗
∂
β
∂
W
∗
∂
γ
.The bulk energy is defined as (23)Fc=κα2β2+γ2+ν2+β2γ2+ν2+γ2ν2+α2β2γ2+ν2+β2γ2ν2,where ν=1-α-β-γ is the fourth phase and κ is a constant. Here, W∗(u(c)) is defined as (24)W∗uc=12∫Ωϵiju:EijklcϵkludΩ.More details about the modeling and simulation can be found in [16].Figure10 shows a schematic of a bridge structure. The design domain has a length to height ratio of 2 : 1. The left bottom corner is fixed while the right bottom corner is simply supported.Figure 10
Design domain of topology optimization.Figure11 shows the system changes from a random initial mix to a final optimal bridge structure. The volume ratio for the three hard material phases is 0.1 : 0.1 : 0.2 and the ratio of Young’s modulus is 1 : 3 : 9. Δt=0.1h, η=0.05, and ϵ=0.005 are used. The FEM mesh is of 64×64 quadrilateral elements.Figure 11
Solution results for the four-material bridge structure at given time stepsm: material 1 in green (softest), material 2 (softer) in blue, material 3 (hardest) in red, and the void in black (reprinted from Zhou and Wang [16], with permission from Springer).
(a)
m = 0 (b)
m = 100 (c)
m = 1000 (d)
m = 2000 (e)
m = 5000 (f)
m = 10000


## 9. Conclusions

In this article, we reviewed and described the basic mechanism of the phase-field modeling with the CH equation. It is important to understand the fundamental dynamics of the CH equation in various models since we can use it to model physical phenomena involving deformable interfaces. As a future work, it would be interesting to review the multicomponent CH equation.



---

*Source: 9532608-2016-10-27.xml*