Abstract
We study a linear stability analysis for the Cahn–Hilliard (CH) equation at critical and off-critical compositions. The CH equation is solved by the linearly stabilized splitting scheme and the Fourier-spectral method. We define the analytic and numerical growth rates and compare these growth rates with respect to the different average levels. In this study, the linear stability analysis is conducted by classifying three average levels such as zero average, spinodal average, and near critical point levels of free energy function, in the one-dimensional (1D) space. The numerical results provide insight for the dynamics of CH equation at critical and off-critical compositions.
1. Introduction
The Cahn–Hilliard (CH) equation describes the temporal evolution of the conserved phase-field by the following partial differential equation [1, 2]: where is a phase-field at position and time , is the free energy per unit volume and can be substituted for other free energy functions in numerical experiments, and constant determines the thickness of the interface. The CH equation (1) was originally introduced as the mathematical model which describes the spinodal decomposition in a binary mixture. This equation has been applied in diverse fields such as dealloying [3, 4], two-phase fluid flow [5–7], topology optimization [8–10], population dynamics [11], tumor growth [12–14], thin films [15–17], block copolymer [18, 19], and image processing [20, 21]. The CH equation is mass-constrained gradient flow in the dual space of zero average space of the Ginzburg–Landau free energy,
To confirm the total mass conservation and the energy dissipation, we differentiate the total mass of and equation (2) with respect to time and obtain the following results, respectively.
The results of equations (3) and (4), respectively, show that the conservation of total mass holds, and the total energy does not increase with respect to time. Figure 1 represents a double-well potential , its corresponding first and second derivatives. Here, the spinodal region is defined by the gray-colored region where the second derivative is negative, i.e., . In this region, is unstable and then phase separation occurs. With these various applications of the CH equation, many researchers have proposed new efficient numerical methods [22–24]. Usually, in the stage of verifying the numerical method, the linear stability analysis is used. In [25], the linear stability was used for local equilibrium of Boltzmann equation. In [26], the stability analysis of the advective CH equation was employed to define the perturbation to investigate the instability of wave packets. In this paper, we shall present a clear standard for this test by defining linear and nonlinear regimes.
The main purpose of this paper is to numerically investigate the dynamics of phase separation in binary mixtures with different average concentrations using the CH equation in the one-dimensional (1D) space. In addition, we propose criteria for linear regimes with small differences between numerical solutions and linear solutions.
This paper is organized as follows. The Fourier-spectral method is introduced in Section 2 for numerical solution. Section 3 provides numerical results and simulations with linear stability analysis. Conclusions are discussed in Section 4.
2. Numerical Solution
Now, we use the Fourier-spectral method [27] to find the numerical solution for the CH equation (1) in the 1D space ,
Let be the grid size in the -axis directions, where is positive even integer. At the edge points where , we simply denote the numerical representation into . Here, and is the temporal step size. Also, we set for the Fourier-spectral method. Let us consider the discrete Fourier transform and the inverse discrete Fourier transform
For simplicity of notation, we define . By applying the linearly stabilized splitting scheme [28] into the CH equation (5), we obtain the following discrete CH equation:
Letting , equation (7) comes into
Then, we can transform equation (8) into the discrete Fourier space as follows:
By rearranging the above one, we obtain
As a consequence, we can compute the updated numerical solution by applying the inverse discrete Fourier transform (6).
3. Numerical Experiments
In this section, we define growth rate of the CH equation and compare analytic growth rate with numerical growth rate. We also define an error between the analytic solution and the numerical solution. Through several numerical experiments about the growth rate and the error, we investigate the numerical behaviour at different spinodal points. Then, we propose linear and nonlinear regimes.
3.1. Analytic and Numerical Growth Rates
In this section, we shall find the analytic and numerical growth rates. We consider the following three different free energy functions [29], , , and at .
We linearize each of , , and by applying the Taylor expansion for different free energy functions; then, we obtain the following linearized terms. The linearized functions are explained in detail in Appendix.
Therefore, for the free energy functions , , and , the corresponding linear Cahn–Hilliard (LCH) equations are as follows:
To derive the ordinary differential equation for amplitude function , we take . For the positive even integer , it can satisfy the periodic boundary condition. Substituting above into LCH equations (17), (18), and (19), we have
Dividing on the both sides of equations (20), (21), and (22), we obtain
The analytic solutions of ordinary differential equations (23), (24), and (25) are as follows:
In this paper, we use unless the free energy function is mentioned; then, we write , . In the numerical approach, numerical growth rate is defined by
Here, denotes the maximum norm.
3.2. Convergence Test
In this section, we verify convergence of the numerical growth rate. For numerical test, we use , , and in .
Table 1 shows the numerical growth rate at time for various and . Here, the analytic growth rate is . In Table 1, we can see that the numerical growth rate converges to the analytic one as and decrease.
3.3. Linear Stability Analysis
From now on, we implement the linear stability analysis with the constant solution . In the every numerical experiment, we use the constant to adjust the interfacial thickness. As long as do not more mentioned, the numerical parameters use such as , , =e-, , , and is total time. In addition, we use three different value , and which are in spinodal region and near critical points of free energy. Note that these values are indicated points , and in Figure 1. For a positive even integer , the initial condition is selected as in . First, we investigate analytic and numerical growth rates for different mode numbers .
Before the next tests, we define linear solution as follows: where and is analytic growth rate, see equation (26). Figure 2 shows the change of the analytic and numerical growth rates with respect to each mode numbers for , , and in short-time evolution. As we can see, numerical growth rate estimates analytic growth rate well in short time (). The growth rates in spinodal region (, ) have positive values for some mode . On the other hand (), the growth rates are always negative regardless of mode number and also growth rates are monotonically decreased as mode number increases. Figures 3(a)–3(c) show the analytic growth rates according to mode numbers for , , and , respectively, where the analytic growth rates , , and are (26), (27), and (28). However, in the long time evolution, it is difficult to expect the numerical growth rate that fits well with the analytic one. In order to find this reason, we first observe dynamics of the numerical solution and linear solution . In Figure 4, we can check the behaviour of the numerical and linear solution in short and long time evolution for and , when and . Unlike the short-time cases, the long-time numerical solutions in all cases show a lot of difference from the linear solutions. This result means that over time, the CH solution has nonlinearity.
(a)
(b)
(a)
(b)
(c)
Let us consider Ginzburg–Landau free energy (2) by decomposing it as where
We compare short- and long-time numerical solution when and . In this test, we use mode and . As shown in Figure 5(a), we can see the nonlinearity in the long-time solution after a certain period of time, regardless of . Figures 5(b)–5(d) show results of double well potential , , and , respectively. Here, the solid and dashed lines mean the short- and long-time results of each function, respectively. Also, the area of the gray-colored region in Figures 5(b)–5(d) represent the energy , , and the Ginzburg–Landau free energy , respectively.
As shown in Figure 1, as moves away from , becomes smaller than when . For this result, in Figure 5(b), dashed line is located below solid line overall of . It means that becomes smaller over time in the spinodal region. Also, according to fluctuation of , we can see the property of in Figure 5(c). However, it can also be confirmed that the value has a slight effect on as shown in Figure 5(d). In conclusion, we know that the CH solution grows with time in the spinodal region, and it has nonlinearity in the existing solution in order to reduce the total free energy.
Figures 6(a)–6(c) show temporal evolution of , , and with different . All results of and decrease with respect to time . However, it can be seen that the larger the epsilon, that is, the thicker the interface of numerical solution, the slower the result is expressed.
From now on, we verify that the CH solution has nonlinearity over time. Therefore, when performing a linear stability test, it is accurate to examine short time evolution. Also, based on this fact, we propose the criteria for defining the following linear and nonlinear regimes. To find linear regime where numerical solution estimates linear solution well, we define error between numerical and linear solution as follows:
We can distinguish linear regime and nonlinear regime based on the user-criteria as follows: if ; then, linear regime; otherwise, non-linear regime. We propose a criterion using the initial condition as follows:
Figure 7 shows proposed criterion and temporal evolution of error for the time with dynamic of linear and numerical solutions.
4. Conclusions
We studied a linear stability analysis for the CH equation in spinodal region as different average levels. To solve the CH equation, we used linearly stabilized splitting scheme and Fourier-spectral method. Through the numerical simulations, we observed various dynamics of the CH equation and confirmed the numerical solution compared linear solution over time. We defined growth rate and also compared numerical growth rate and analytic growth rate. Using difference between numerical and linear solution, we defined error. By means of defined error, we distinguish linear regime where numerical solution and linear solution match because of the small error and nonlinear regime where do not.
Appendix
A. Appendix: Linearlization
By the Taylor theorem, for real valued-function , if has derivatives of all orders at , then for each positive integer , where .
Now, we approximate the linearized functions of various free energies at .
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The first author (S. Ham) was supported by the Brain Korea 21 FOUR through the National Research Foundation of Korea funded by the Ministry of Education of Korea. This study was supported by 2020 Research Grant from Kangwon National University. The author (D. Jeong) was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (NRF-2020R1F1A1A01075937). The corresponding author (J.S. Kim) was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF).