Basic Principles and Practical Applications of the Cahn – Hilliard Equation

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.


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][5][6], image inpainting [7], multiphase fluid flows [8][9][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: where Ω ⊂ R  ( = 1, 2, 3) is a bounded domain with a boundary Ω.The quantity (x, ) =   −   is a phase-field defined as the difference between the mole fractions of binary mixtures, where   and   are the mole fractions of phases  and .() = 0.25( 2 − 1) 2 is a double-well potential of a homogeneous system of composition  as shown in Figure 1.((x, )) is a positive mobility,  is a positive constant, and / is the outward normal derivative at the domain boundary.Originally, the concentration-dependent mobility has a form () =  0 (1 − )(1 + ), where  0 is a constant.This form has the same property in diseases transmission in that transmission happens by contact.
The CH equation can be derived from the total free energy functional equation:   = −∇⋅J, where the flux is defined as J = −∇.The CH system (1) possesses the following two properties: 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 Section 2, 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.

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.

Diblock Copolymer
A diblock copolymer is a linear-chain molecule consisting of two subchains joined covalently to each other [4].Let (x, ) =   (x, ) −   (x, ) be the order parameter which is defined as the difference between the local volume fractions of  and  monomers at point x and time .Then, the governing equation by the Ohta-Kawasaki model [28] is as follows: where Ω ⊂ R  ( = 1, 2, 3) is a domain,  = ∫ Ω  x/|Ω| 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: where () = 0.25( 2 − 1) 2 and (x, y) denotes Green's function of −Δ on Ω with periodic boundary conditions [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 (x), where x = (,), be a given binary image in a domain Ω and  ⊂ Ω be the inpainting domain.The image is scaled so that 0 ≤  ≤ 1.Let (x, ) be a phase-field which is governed by the following modified CH equation: where () = 0.25 2 (1 − ) 2 and (x) = 0 if x ∈ ; otherwise (x) =  0 .We perform a test of the inpainting of a double stripe as shown in Figure 4(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, ℎ = 1/128, and Δ = 0.1.And then we switch these values to  0 = 0.1,  = 0.1, and Δ = 1 after 3000 iterations and stop the simulation after 4000 iterations.Figures 4(b), 4(c), and 4(d) are numerical solutions at  = 40, 500, and 1300, respectively.The final inpainting result represents the completed stripe.Numerical solution algorithms can be found in [7,30].

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 where u is the velocity,  is the pressure, () is the density, SF is the surface tension force, g is the gravitational force,  is the phase-field function, () is the mobility function, and  is the chemical potential.subscript * denotes the corresponding characteristic values.
As a test problem, we consider a drop falling test under the gravity.A schematic diagram of a two-phase domain is shown in Figure 5. Figure 6 shows the interfaces of the falling drop at time  = 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 Δ = 0.0005 and ℎ = 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.

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 field (x) and elastic strain  el  (x).The total free energy in an inhomogeneous system is defined as Here, () is the local chemical free energy density defined by a quartic polynomial () = 2.5( −   ) where u = (,V) is the displacement and  11 ,  12 , and  44 are the cubic elastic parameters and are dependent on the order parameter .Here,   =  1   +  0  (1 − ) for  = 11, 12, and 44, where  0  and  1  are elastic constants of matrix and precipitate, respectively.The eigenstrain is () = ( −   ), which obeys Vegard's law.Here,   is the average composition and  is a constant.And  11 =   ,  22 = V  , and  12 = 0.5(  + V  ). is the gradient energy coefficient which can be related to interatomic interaction parameters.By the variational derivative F/, we obtain the following equation: From the total free energy functional (8), we derive the modified CH equation for the two-dimensional microstructure evolution with strong elastic inhomogeneity: where  is the mobility.From the definition, we get  where the subscripts  and  denote partial derivatives with respect to the corresponding arguments.Now, we perform a numerical test with  = 1.2247,   = 0.053,   = 0.947, the composition expansion coefficient  = 0.05, spatial step size ℎ = 1, temporal step size Δ = 1, and the computational domain Ω = (0, 256) × (0, 256).With an initial condition, we investigate the effect of elastic inhomogeneity .We take the elastic constants  0 11 = 232,  0 12 = 153, and  0 44 = 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  =  11 +  12 and the same ratio of anisotropy  = 2 44 /( 11 −  12 ), while changing the ratio of shear modulus  =   44 /  44 in the matrix and precipitate phases, respectively [39].
Figures 7(a)-7(f) show the numerical results at time  = 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   44 =    44 ; that is, the stiffness of the precipitate is less than the matrix [41][42][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].

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 Ω  , Ω  , and Ω  be the healthy, viable, and dead tumor tissues, respectively.We also define Ω  = Ω  ∪ Ω  as the tumor cell.See Figure 8 for the schematic of these definitions.
Let   and   be the volume fractions of tumor and dead tumor tissues.Then, the nondimensional equations for the volume fraction of tumor cells are where (  ) is a continuous cut-off function.The nondimensional equations for the volume fraction of necrotic cells are where H is the Heaviside function.The nondimensional cell velocity is where  is a nondimensional measure of the adhesion force.The velocity satisfies  The nondimensional quasi-steady nutrient equation is where the diffusion coefficient and nutrient capillary source term are where   is the nutrient diffusivity in the healthy tissue and (  ) 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][46][47].

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    ) .
The bulk energy is defined as where ] = 1 −  −  −  is the fourth phase and  is a constant.
Here,  * (()) is defined as More details about the modeling and simulation can be found in [16].
Figure 10 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 11 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. Δ = 0.1ℎ,  = 0.05, and  = 0.005 are used.The FEM mesh is of 64 × 64 quadrilateral elements.

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.
) and3(b)  show the evolution of the numerical solutions with  = −0.3 and  = 0, respectively.As shown in Figure3(a), figures in the first row change to the hexagonal geometry as time increases.Unlike this, we can see the lamellar geometry in Figure3(b) at the equilibrium state. t

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.

Figure 5 :
Figure 5: Schematic diagram of a two-phase fluid domain.

Figure 7 :Figure 8 :
Figure 7: Morphologies of a circular precipitate with different shear modulus  at time  = 1700.

Figure 9 :
Figure 9: Evolution of the   = 0.5 isosurface during the growth of an asymmetrical 3D tumor (reprinted from Wise et al. [17], with permission from Elsevier).

Figure 10 :
Figure 10: Design domain of topology optimization.

Figure 11 :
Figure 11: Solution results for the four-material bridge structure at given time steps : material 1 in  (softest), material 2 (softer) in , material 3 (hardest) in , and the void in  (reprinted from Zhou and Wang [16], with permission from Springer).