A Review of the Numerical Methods for Diblock Copolymer Melts

1  Introduction

In this review, we present an overview of the numerical methodologies for diblock copolymer melts. The Ohta–Kawasaki equation, introduced by Ohta et al. [1], was originally developed to explain the microphase separation patterns in diblock copolymer melts. A diblock copolymer is defined as a polymer composed of two distinct subchains linked by a covalent bond, with each subchain composed of monomers that repel each other. When a large number of such molecules aggregate, they form a structure known as a melt. These diblock copolymer melts are of great interest in both scientific and engineering fields due to their unique phase separation behavior. Their distinct constituent monomers separate to form nanoscale structures with desirable properties [2]. This phase separation is a process known as spinodal decomposition, which enables the formation of various microstructures such as lamellar, cubic, hexagonal, and gyroid configurations depending on the synthetic conditions used [3]. These diverse structures can significantly influence the mechanical properties of the material [4], and it is therefore essential to anticipate the structure that will be formed experimentally.

Due to these unique characteristics, diblock copolymer melts have found applications across multiple fields, such as materials science [5,6], nanotechnology [7–9], and biomedical engineering [10]. Controlled pattern formation at the nanoscale can improve the mechanical, optical, and electrical properties of materials, and enable the creation of advanced materials tailored for specific functions [11]. Several computational algorithms have been proposed to analyze the equilibrium state of diblock copolymers. Techniques such as Monte Carlo methods [12–14], molecular dynamics [15], self-consistent field theory [16,17], cell dynamics simulations [18], dissipative particle dynamics simulations [19–21] and the phase-field model [22–24] are widely used for this purpose. The phase-field model, in particular, is a continuum model that describes the interface using the order parameter. This model is advantageous because it bypasses the need to solve complex boundary value problems, as it formulates the governing equations through the variation of an energy functional. Consequently, the phase-field model has become a widely adopted approach for reproducing various material structures, and it enables researchers to derive the governing equations in a relatively straightforward manner [25–27]. Researchers also constructed new or modified phase field equations for the diblock copolymer solutions for specific conditions and cases. Shen et al. [28] proposed a thermodynamically consistent phase field equation for the simulation of diblock copolymers. The electric and magnetic field is coupled in the proposed model. The authors constructed the phase field model based on the Onsager principle. In order to solve the proposed model, the energy quadratization technique is applied, while the thermodynamical consistency is preserved. Numerical simulations showed convergence rate of the proposed algorithm and the effects of electric fields and magnetic field on the diblock copolymer system.

The Cahn–Hilliard (CH) equation is a fundamental governing equation in phase-field models and is widely used to describe phase separation phenomena in binary systems. It effectively models the evolution of the compositional concentration field and reveals the role of interface thickness in structure dynamics. The Navier–Stokes equations can model the flow of phases, and various numerical methods have been developed for their solutions [29]. When coupled with the CH equation, the Navier–Stokes equations can be used to model the flow of polymer melts. Song et al. [30] introduced an innovative, unconditionally energy-stable data assimilation algorithm for solving the Navier–Stokes–Cahn–Hilliard equations. Their approach demonstrated robust performance in handling locally discrete observational data. Song et al. [31] proposed an unconditionally energy-stable numerical method for the CH equation and validated its efficiency and stability through numerical experiments involving complex initial conditions. The CH equation is frequently coupled with other equations to capture additional characteristics. Martínez–Agustín et al. [32] proposed coupling the CH equation with the Swift–Hohenberg equation in 3D to model phase transitions in diblock copolymers. The proposed method is solved using the fast Fourier transform and a pseudo-spectral implicit method. Since the model dynamics can produce various morphologies, the porous polymeric materials obtained from numerical simulations have potential applications in 3D printing. As the Ohta–Kawasaki model explains the microphase separation patterns in diblock copolymer melts, Barua et al. [33] rescaled the Ohta–Kawasaki model, free energy, and boundary conditions for more efficient calculations of sharp-interfacial symmetric diblock copolymer problem. They applied boundary integral formulation and small-scale decomposition to solve the governing equations numerically. Their simulations illustrate the evolution of the system from the initial to intermediate states and eventually to the steady. Meng et al. [34] proposed a solution for the nonlocal CH equation with a nonlocal diffusion operator by using a scalar auxiliary variable method. Their proposed methods are unconditionally stable and showed high-order accuracy. Numerical solutions were efficiently obtained using the fast Fourier transform and gradient approach to achieve low storage cost and calculation time. The nonlocal CH equation with degenerate mobility was studied by Elbar et al. [35], where the additional nonlocal term accounts for the surface tension to model long-range interactions.

The Ohta–Kawasaki model, defined on the domain Ω∈Rd for d=1,2,3, is given by

∂ϕ(x,t)∂t=Δμ(x,t)−α(ϕ(x,t)−ϕ¯),(1)

μ(x,t)=F′(ϕ(x,t))−ε2Δϕ(x,t),(2)

where ϕ=ρA−ρB is a phase-field variable representing the local volume fraction of two monomers and ϕ¯=∫Ωϕ(x,0)/∫Ωdx. The parameter α is the molecular chain length and ε is the suitably rescaled interfacial thickness at the different two monomer intersections. The term F(ϕ)=0.25(ϕ2−1)2 denotes the Helmholtz free energy functional. If α>0, the nonlocal Cahn–Hilliard (NCH) equation becomes the Ohta–Kawasaki model. The total free energy ℰ(ϕ) of the Ohta–Kawasaki model can divided into the short-range energy ℰs(ϕ) and long-range energy ℰl(ϕ), defined as follows:

ℰs(ϕ)=∫Ω(F(ϕ)+ε22|∇ϕ|2)dx,(3)

ℰl(ϕ)=α2∫Ω∫ΩG(x−y)(ϕ(x)−ϕ¯)(ϕ(y)−ϕ¯)dydx,(4)

where G is a Green’s function [36]. We can rewrite the long-range energy by considering the periodic or zero Neumann boundary conditions.

ℰl(ϕ)=α2∫ΩΔxω(x)[∫ΩΔyG(x−y)ω(y)]dx=α2∫Ω|∇ω(x)|2dx,(5)

where Δω=−(ϕ−ϕ¯) and ω is derived from Δω. The total free energy is obtained by summing the short-range and long-range energies.

ℰ(ϕ)=∫Ω(F(ϕ)+ε22|∇ϕ|2+α2|∇ω|2)dx.(6)

We also obtain the following energy dissipation law and mass conservation property [37].

ddtℰ(ϕ)=∫Ω[(F′(ϕ)−ε2Δϕ)ϕt−αωΔωt)dx=∫Ω(μ+αω)ϕtdx=−∫Ω|∇(μ+αω)|2dx≤0,(7)

and

ddt∫Ωϕdx=∫Ω[Δμ−α(ϕ−ϕ¯)]dx=∫∂Ω∂μ∂nds−α∫Ω(ϕ−ϕ¯)dx=0,(8)

where the integral by parts and periodic or zero Neumann boundary conditions are applied. Therefore, the Ohta–Kawasaki model satisfies the energy dissipation law and mass conservation by Eqs. (7) and (8), respectively. Because of these two fundamental properties of the Ohta–Kawasaki model, it is natural and important to construct proper numerical schemes that still satisfy these two properties in discrete manners. Singh et al. [9] developed a thermodynamically consistent nonlocal model for phase transformation and heat transfer in block copolymer directed self-assembly using the continuum theory of mixtures. Their model incorporates mass and energy balances, microforce balances, and constitutive relations to couple the phase field and temperature.They developed a finite element solution incorporating stable time-stepping schemes and performed computational tests to validate the equation’s effectiveness in simulating the self-assembly behavior of block copolymers. Chen et al. [38] proposed an efficient leapfrog time-marching method for the phase-field diblock copolymer equation, which significantly reduces the computational cost by requiring only the solution of a linear algebra system at each time step, ensures unconditional energy stability for large time steps, and guarantees the existence and uniqueness of the computational solution at each step, ensuring both reliability and accuracy. This method consists of the leapfrog integration scheme for time discretization and the Fourier spectral method for spatial discretization. The authors analyze the energy dissipation and mass conservation of the proposed numerical method for the Ohta–Kawasaki model and prove that the numerical scheme satisfies both properties. Li et al. [39] developed a highly efficient computational method for a flow-coupled phase-field equation for diblock copolymer melts. This model represents a complex nonlinear system consisting of the Navier–Stokes equations coupled with the CH type equation with an Ohta–Kawasaki potential. By combining decoupling techniques with the projection method, the authors developed a fully decoupled, energy-stable, second-order time-accurate computational method for the model. The decoupling technique is based on designing an auxiliary ordinary differential equation (ODE) that plays a critical role in achieving a fully decoupled structure while maintaining energy stability. The authors rigorously proved that the method satisfies unconditional energy stability. Diblock copolymers under various conditions are studied by Wang et al. [40], authors simulated a hydrodynamically coupled diblock copolymer in complex domains using a diffusion domain method which efficiently overcomes the problems derived from complex boundaries. Therefore, the computational efficiency is highly increased. In addition, 2D and 3D models are also considered using the second-order dimension splitting technique, which allows multi dimension problems to be decomposed into a 1D problem. Various numerical simulation tests in 1D, 2D and 3D verified the stability, accuracy and efficiency of the proposed method. The Hele–Shaw cell simulates a flow in a narrow gap between two flat parallel plates, which can model interfacial pattern formation in nonequilibrium systems [41]. Modeling a diblock copolymer melt in the Hele-Shaw cell is a complicated nonlinear system. Cao et al. [42] studied this system using a combination of explicit-invariant energy quadratization and projection methods. The proposed method is second order accurate in time and energy stable. Due to the introduction of two auxiliary variables, high efficiency is also achieved. Numerical simulations with random initial conditions confirmed various characteristics of the proposed method such as accuracy, stability. The spinodal decomposition in the rotating Hele–Shaw cell is also simulated using the proposed scheme. Xu et al. [43] studied the evolution of diblock copolymer in 2D domain where the governing equation is the extended Ohta–Kawasaki model. The projection operator is introduced to achieve mass conservation property of the proposed method. The projection method is constructed with a small computational cost, therefore simplifying the numerical algorithm. Zhang et al. [44] studied the microphase separation of diblock copolymers in three-dimensional space using a modified NCH equation, which is a variant of the CH equation with a viscous term. The viscous term introduces viscous effects to the governing equation. The integrating factor Runge–Kutta technique is adopted for large tie-stepping. Numerical simulation of diblock copolymers in 3D shows that the proposed scheme is capable of distinguishing different phase geometries. The maximum principle, mass conservation and energy stability is also satisfied in this numerical test. Luo et al. [45] studied the 2D and 3D nonlocal Ohta–Kawasaki model using the Fourier spectral method, where the nonlocal Ohta–Kawasaki model studies the pattern formation for the diblock copolymer system. The second-order backward finite difference method (FDM) is applied for Fourier collocation discretization. The asymptotic compatibility of the method is mainly discussed in this paper, proving in both analytic and numerical approaches. Numerical simulations also confirmed characteristics of the proposed method such as upper boundedness, effect of parameters, convergence rate and energy stability. Iqbal et al. [46] introduced a cell dynamic simulation model that predicts simulation outcomes through an examination of flow, deformation, and phase transitions in diblock copolymer systems under curvilinear coordinate systems.

In this review, we focus on various numerical methods and applications for solving the NCH equation for microphase separation patterns in diblock copolymer melts. The overall structure of the paper is as follows. Section 2 presents detailed explanations of the various numerical schemes for solving the NCH equation. Section 3 introduces the applications of the phase-field model for diblock copolymer melts. Section 4 presents the conclusion.

2  Numerical Methods

We describe the various numerical methods used to solve the NCH equation. We consider the computational domain Ω=(Lx,Rx)×(Ly,Ry) in 2D space. Let Δt be time step. We define h=(Rx−Lx)/Nx=(Ry−Ly)/Ny. The discrete domain is defined as Ωd={(xi=Lx+ih,yj=Ly+jh) | i=0,1,…,Nx, j=0,1,…,Ny} or Ωd={(xi=Lx+(i−0.5)h,yj=Ly+(j−0.5)h) | i=1,2,…,Nx,j=1,2,…,Ny} depending on whether we use cell corner points or cell center points. Fig. 1 shows the schematic illustrations of the cell corner points and cell center points on the computational domain. In this review paper, when boundary conditions are not explicitly specified, a periodic boundary condition is applied for a discrete domain with cell corner points, as shown in Fig. 1a, while a Neumann boundary condition is applied for a discrete domain with cell center points, as shown in Fig. 1b. Let ϕijn be the approximation of the ϕ(xi,yj,nΔt). The discrete l2-norm in is defined by

‖ϕ‖2=1NxNy∑i=1Nx∑j=1Ny(ϕij)2.

Figure 1: Schematic illustrations of (a) cell corner points and (b) cell center points on the computational domain

We can easily reduce or extend to one-dimensional and three-dimensional spaces.

When reviewing the numerical methods, we mainly focus on the accuracy and stability of the reviewed methods. High accuracy allows the user to predict phase separation dynamics in diblock copolymer systems with relatively low computational cost, such as coarse spatial grid or large time step. Therefore, high accuracy can lead to an effective numerical method. Stability is also concerned with time step and accuracy of the numerical method. High stability, and sometimes unconditional stability allows a larger time step to be used for the numerical simulation without the risk of blowing up. Therefore, we can obtain more freedom when modifying the time step size for a desired accuracy.

2.1 Finite Difference Method

The FDM is a widely employed numerical technique for solving partial differential equations, known for its simplicity and broad applicability. A key advantage of this method is its ease of implementation for a range of problems with regular grid structures, making it particularly effective for time-dependent simulations. Although the accuracy of the method is influenced by grid resolution and step size, it can be refined to achieve higher precision. Despite some limitations in handling complex boundaries or irregular domains, the FDM remains a valuable tool due to its computational efficiency and straightforward implementation. Xiao et al. [47] developed a space-time fourth-order method for two- and three-dimensional CH type equations. The authors used the operator splitting method with auxiliary variables for spatial differentiation terms to enable multi-thread computation. In addition, we extended the numerical scheme for the phase-field diblock copolymer model. Numerical methods for solving systems of equations with nonlinear terms, such as the NCH equation, are primal and difficult problems [48]. Jeong et al. [22] conduct a numerical investigation on controlling local defectivity in self-assembled diblock copolymer patterns by designing suitable substrates. The numerical solution algorithm is described as follows. Considering trench domain Ω with boundaries ∂Ω1 and ∂Ω2 shown in Fig. 2. For ϕ, we use the Dirichlet and periodic boundary conditions in ∂Ω1 and ∂Ω2, respectively. For μ, we apply the Neumann and periodic boundary conditions in ∂Ω1 and ∂Ω2, respectively.

Figure 2: Schematic diagram of the computational domain Ω with boundaries ∂Ω1 and ∂Ω2

We discretize Ω using cell corner points. The NCH Eqs. (1) and (2) with non-linearly stabilized splitting scheme [49] is rewritten as

ϕijn+1−ϕijnΔt=Δμijn+1−α(ϕijn+1−ϕ¯),(9)

μijn+1=(ϕijn+1)3−ϕijn−ε2Δϕijn+1,(10)

where ϕ¯=∑xij∈Ωϕij0/∑xij∈Ω1. To solve the NCH Eqs. (9) and (10), we apply the Gauss–Seidel iterative scheme. Given ϕijn, let ϕijn+1,0=ϕijn be an initial guess. We calculate ϕijn+1,m+1 and μijn+1,m+1 from ϕijn+1,m and μijn+1,m for m=0,1,… by

(1Δt+α)ϕijn+1,m+1+4h2μijn+1,m+1=ϕijnΔt+αϕ¯+μi−1,jn+1,m+1+μi+1,jn+1,m+μi,j−1n+1,m+1+μi,j+1n+1,mh2,(11)

(−4ε2h2−3(ϕijn+1,m)2)ϕijn+1,m+1+μijn+1,m+1=−ϕijn−2(ϕijn+1,m)3−ε2ϕi−1,jn+1,m+1+ϕi+1,jn+1,m+ϕi,j−1n+1,m+1+ϕi,j+1n+1,mh2.(12)

We need to apply certain special formulas considering the boundary conditions. To facilitate understanding, we consider one example illustrated in Fig. 3. For the points ϕij near the boundary within the domain Ω, the values of ϕa and ϕb at ∂Ω1, as shown in Fig. 3a, are determined by the Dirichlet boundary condition and correspond to the maximum values of the numerical equilibrium solution, ‖ϕeq‖∞. Moreover, as shown in Fig. 3b, the values of μa and μb at ∂Ω1 are determined by the Neumann boundary condition and correspond to μp and μq, respectively, which are obtained by the linear interpolation of neighboring points. Thus, μa=μq and μb=μp. Here, μa and μp have the same values and are located at distinct positions. Similarly, μb and μq also share identical values and are situated at separate locations. Other points near the boundary are similar. We can obtain the values of ϕa and ϕb at the boundary ∂Ω1 using the Dirichlet boundary condition. We define the numerical equilibrium solution as ϕk if ‖ϕk−ϕk−1‖2≤10−6 for k≥1. The discrete second derivatives near boundary for the Dirichlet boundary condition are defined by

ΔxxDϕij=(ϕa−ϕijαh−ϕij−ϕi−1,jh)/(αh+h2),ΔyyDϕij=(ϕb−ϕijβh−ϕij−ϕi,j−1h)/(βh+h2),

where ϕa and ϕb are the maximum values of the numerical equilibrium solution and 0<α,β<1. Then, for the Dirichlet boundary condition, we define the discrete Laplace operator as ΔDϕijn+1=ΔxxDϕijn+1+ΔyyDϕijn+1. Next, we define the discrete second derivatives near boundary for the Neumann boundary condition as

ΔxxNμij=(μa−μijαh−μij−μi−1,jh)/(αh+h2),ΔyyNμij=(μb−μijβh−μij−μi,j−1h)/(βh+h2),

where μa=μp=pμij+(1−p)μi,j−1 and μb=μq=(1−q)μij+qμi−1,j by the linear interpolation. Therefore, for the Neumann boundary condition, the discrete Laplace operator is defined by ΔNμijn+1=ΔxxNμijn+1+ΔyyNμijn+1. We repeat the above iterations until the l2-norm error between numerical solutions at successive time steps is smaller than a given tolerance tol: ‖ϕn+1,m+1−ϕn+1,m‖2<tol. We perform the numerical simulations using the parameters L∗=0.375, h=L∗/10, Nx=500, Ny=300, Δt=0.1h, ε=1/(252), α=100, the maximum values of the numerical equilibrium solution ‖ϕeq‖∞=0.6134. The initial condition on the computational domain Ω=(−25L∗,25L∗)×(−15L∗,15L∗) is defined as ϕ(x,y,0)=0.01rand(x,y). Fig. 4 shows the numerical solution solved using the numerical method [22].

Figure 3: Schematic representation of (a) Dirichlet boundary condition and (b) Neumann boundary condition

Figure 4: Snapshots of the numerical solutions with (a) θ=30 and (b) θ=60

Jeong et al. [24] developed a numerical method to investigate microphase separation patterns in diblock copolymer melts on curved surfaces. This method employs a discrete narrow band grid adjacent to the curved surface and applies a pseudo-Neumann boundary condition for the near boundary using the closest point scheme. Therefore, the Laplace–Beltrami operator can be replaced with the standard Laplace operator. We define the σ-neighborhood band domain as Ωσ={(x,y,z)∣(a,b,c)∈S,(x,y,z)=(a,b,c)+θn(a,b,c) for |θ|<σ}, which includes the surface S. Here, n(a,b,c) is the unit normal vector at (a,b,c) on the surface S. The global domain that includes the Ωσ is defined by Ω=(Lx,Rx)×(Ly,Ry)×(Lz,Rz). Let Nx, Ny, and Nz be positive integers and let h=(Rz−Lz)/Nz=(Ry−Ly)/Ny=(Rx−Lx)/Nx be a size of spatial step. The discrete global domain is defined as Ωd={(xi=Lx+ih,yj=Ly+jh,zk=Lz+kh) | i=0,1,…,Nx, j=0,1,…,Ny,k=0,1,…,Nz}. Let ϕijkn and μijkn be approximations of ϕ(xi,yj,zk,nΔt) and μ(xi,yj,zk,nΔt), respectively. The discrete narrow band domain is defined by Ωσd={(xi,yj,zk):|ψijk|<σ}. To find the boundary points of discrete narrow band domain Ωσd, define Pijk=0 if (xi,yj,zk)∈Ωσd and Pijk=1 otherwise. Let ∂Ωσd={(xi,yj,zk):Pijk|∇Pijk|≠0} be the boundary points of discrete narrow band domain. Here, the gradient of P is defined by

∇Pijk=(Pi+1,jk−Pi−1,jk,Pi,j+1,k−Pi,j−1,k,Pij,k+1−Pij,k−1)2h.

We define the discrete Laplace operator as

Δϕijkn=ϕij,k+1n+ϕij,k−1n+ϕi,j+1,kn+ϕi,j−1,kn+ϕi+1,jkn+ϕi−1,jkn−6ϕijknh2.

Then, we discretization the NCH Eqs. (1) and (2) by applying the unconditionally stable scheme.

ϕijkn+1−ϕijknΔt=Δμijkn+1−α(ϕijkn+1−ϕ¯),(13)

μijkn+1=(ϕijkn+1)3−ϕijkn−ε2Δϕijkn+1.(14)

The numerical closest point of xijk=(xi,yj,zk) on S is defined by

cp(xijk)=xijk−∇|ψijk||∇|ψijk|||ψijk|.

We use pseudo-Neumann boundary condition on ∂Ωσd as

ϕijkn+1=ϕn+1(cp(xijk))  and  μijkn+1=μn+1(cp(xijk)).

Fig. 5a,b shows schematic illustrations of the narrow band domain for the cross section of the sphere in three- and two-dimensions, respectively. Since cp(xijk) is not included in the grid points above the discrete narrow band domain, we compute ϕn(cp(xijk)) using the trilinear interpolation. Then, we solve the discrete NCH equation using the Jacobi iterative method [50], which iteratively updates the solution until convergence. Given ϕijkn and μijkn on Ωσd, the initial guesses ϕijkn+1,0 and μijkn+1,0 are used as starting points for the iterative process and are set to ϕijkn and μijkn, respectively. Specifically, ϕijkn+1,0=ϕijkn and μijkn+1,0=μijkn. Then, we calculate the Eqs. (13) and (14) by the following Jacobi iteration for m≥0.

Figure 5: Schematic illustrations of (a) the narrow band domain for surface S and (b) the closest points for the boundary ∂Ωσ of narrow band domain

(1Δt+α)ϕijkn+1,m+1+6h2μijkn+1,m+1=ϕijknΔt+αϕ¯+1h2(μi+1,jkn+1,m+μi−1,jkn+1,m+μi,j+1,kn+1,m)+1h2(μi,j−1,kn+1,m+μij,k+1n+1,m+μij,k−1n+1,m),(15)

−(3(ϕijkn+1,m)2+6ε2h2)ϕijkn+1,m+1+μijkn+1,m+1=−ϕijkn+2(ϕijkn+1,m)3+ε2h2(ϕi+1,jkn+1,m+ϕi−1,jkn+1,m)+ε2h2(ϕi,j+1,kn+1,m+ϕi,j−1,kn+1,m+ϕij,k+1n+1,m+ϕij,k−1n+1,m),(16)

Here, we iteratively compute Eqs. (15) and (16) until ‖ϕn+1,m+1−ϕn+1,m‖2<tol on Ωσd is satisfied for given tolerance tol. Therefore, we obtain the ϕijkn+1 on the narrow band domain Ωσd. For more details of the numerical algorithm, the reader is referred to [24].

Fig. 6 displays the temporal evolution of the computational solution for the bunny surface with initial conditions ϕ(x,0)=−0.3+0.1rand(x), Nx=Ny=Nz=138, h=0.031, Δt=0.1, ε=1/(302), and α=150, where rand (x) is the random value from −1 to 1 at the point x.

Figure 6: Snapshot of the numerical solution using the numerical algorithm from [24] at time t=20Δt, 25Δt, and 600Δt on the bunny surface

We observe that numerical simulation results form appropriate patterns in complex bunny surfaces. Yang [23] presented a linear time-marching scheme for the Ohta–Kawasaki model coupled with incompressible fluid flow, describing the phase-field model for diblock copolymers in fluid environments. The numerical algorithm employs the scalar auxiliary variable (SAV) approach, which ensures energy stability, even with large time steps. The 2D and 3D spatial discretizations are conducted using the FDM, providing a practical and efficient framework for computation. Furthermore, the authors analytically demonstrated the existence of unique solutions and proved the method’s energy stability.

2.2 Fourier Spectral Method

The Fourier spectral method is a technique for solving partial differential equations, that approximates a solution using a sum of functions from a certain function space. To solve problems with periodic boundary conditions, use Fourier transformation composed of sine and cosine functions, and to solve problems with Neumann boundary conditions, use cosine transformation composed only of cosine functions. We consider only two dimensions, and for one or three dimensions it is easily derived from this. The Fourier transformation satisfies the periodic boundary condition, thus we define the discrete domain using the cell corner points. The discrete Fourier transform and its inverse transform are defined by

ϕ^pqn=∑j=1Ny∑i=1Nxϕijne−i(xiξp+yjηq),(17)

ϕijn=1NxNy∑q=1−Ny/2Ny/2∑p=1−Nx/2Nx/2ϕ^pqnei(xiξp+yjηq),(18)

where ξp=2πp/(Rx−Lx) and ηq=2πq/(Ry−Ly). Let ϕ(x,y,nΔt) be sufficiently smooth continuous extension counterpart of ϕijn. The partial derivatives in Fourier space for the variables x and y are

∂ϕ(x,y,nΔt)∂x=1NxNy∑q=1−Ny/2Ny/2∑p=1−Nx/2Nx/2(iξp)ϕ^(ξp,ηq,nΔt)ei(ξpx+ηqy),(19)

∂ϕ(x,y,nΔt)∂y=1NxNy∑q=1−Ny/2Ny/2∑p=1−Nx/2Nx/2(iηq)ϕ^(ξp,ηq,nΔt)ei(ξpx+ηqy).(20)

Thus, Laplacian is defined using Eqs. (19) and (20) as

Δϕ(x,y,nΔt)=∂2ϕ(x,y,nΔt)∂x2+∂2ϕ(x,y,nΔt)∂y2=−1NxNy∑q=1Ny∑p=1Nx(ξp2+ηq2)ϕ^(ξp,ηq,nΔt)ei(ξpx+ηqy).

Then, the linearly stabilized splitting method is used to the NCH equation as

ϕijn+1−ϕijnΔt=Δ(2ϕijn+1+fijn−ε2(Δϕn+1)ij)ij−α(ϕijn+1−ϕ¯),(21)

where fijn=(ϕijn)3−3ϕijn. Thus, we applied the discrete Fourier transformation to Eq. (21).

ϕ^pqn+1−ϕ^pqnΔt=−(ξp2+ηq2)(2ϕ^pqn+1+ε2(ξp2+ηq2)ϕ^pqn+1+f^pqn)−α(ϕ^pqn+1−ϕ¯).(22)

Hence, we obtain the numerical solution in the Fourier space from Eq. (22) as

ϕ^pqn+1=ϕ^pqn+Δt(αϕ¯1^pq−(ξp2+ηq2)f^pqn)1+Δt(α+2(ξp2+ηq2)+ε2(ξp2+ηq2)2).

Then, we use the inverse discrete Fourier transform to obtain the computational solution ϕijn for i=1,2,…,Nx, j=1,2,…,Ny.

ϕijn+1=1NxNy∑q=1−Ny/2Ny/2∑p=1−Nx/2Nx/2ϕ^pqn+1ei(xiξp+yjηq).

The Fourier spectral method based on discrete cosine transformation can be used to solve phase-field models with Neumann boundary conditions [51]. We define the computational discrete domain using the cell center points to describe the Fourier spectral method with discrete cosine transform. The discrete cosine transform and its inverse transform are defined by

ϕ^pqn=upvq∑j=1Ny∑i=1Nxϕijncos⁡(ξpx^i)cos⁡(ηqy^j)ϕijn=∑q=1Ny∑p=1Nxupvqϕ^pqncos⁡(ξpx^i)cos⁡(ηqy^j),

where x^i=xi−Lx, y^j=yj−Ly, ξp=π(p−1)/(Rx−Lx), ηq=π(q−1)/(Ry−Ly), u1=1/Nx, v1=1/Ny, up=2/Nx for p=2,3,…,Nx, and vq=2/Ny for q=2,3,…,Ny. Therefore, the partial derivatives in the Fourier space for the variables x and y are defined by

∂ϕ(x,y,nΔt)∂x=−∑q=1Ny∑p=1Nxupvqξpϕ^(ξp,ηq,nΔt)sin⁡(ξpx^)cos⁡(ηqy^)∂ϕ(x,y,nΔt)∂y=−∑q=1Ny∑p=1Nxupvqηqϕ^(ξp,ηq,nΔt)cos⁡(ξpx^)sin⁡(ηqy^),

where x^=(x−Lx) and y^=(y−Ly). We get the Laplacian by using the partial derivatives in Fourier space:

Δϕ(x,y,nΔt)=∂2ϕ(x,y,nΔt)∂x2+∂2ϕ(x,y,nΔt)∂y2=−∑q=1Ny∑p=1Nx(ξp2+ηq2)upvqϕ^(ξp,ηq,nΔt)cos⁡(ξpx^)cos⁡(ηqy^).

Hence, similar to the Fourier spectral method with discrete Fourier transform, we can solve Eq. (21) using discrete cosine transformation. Xia et al. [52] employed a phase-field model within a Lagrange multiplier framework to investigate crystal phase transitions and nucleation processes, demonstrating the effectiveness of the Fourier spectral method in capturing complex phase dynamics. Li et al. [53] used the Fourier spectral method to solve biological transport networks in complex domains, providing insights into the optimization properties and adaptive mechanisms of network structures. Refer to [54] for detailed information. Fig. 7 shows the schematic illustrations of the discrete Fourier transform, the discrete cosine transform, and the discrete sine transform, from top to bottom. The dots on the x-axis represent the points of the discrete domain, while the dashed lines connecting the dots indicate the corresponding values of ϕ(x) and the solid line is the sum of the basis functions for each transformation that includes all the discrete points. We can observe that the curves represent the sums of the bases of each discrete transform, satisfying the respective boundary conditions.

Figure 7: Schematic illustrations of (a) the discrete Fourier, (b) the discrete cosine, and (c) the discrete sine transforms

Jeong et al. [55] investigated the energy-minimizing wavelength in the equilibrium state of diblock copolymers in the hex-cylinder phase by solving the NCH equation using the Fourier spectral method with the discrete Fourier transform. The boundary condition is a periodic boundary condition with the discrete domain using the cell center point. The authors performed computations in rectangular domains with an aspect ratio of 3 until a numerical equilibrium state was reached, and repeated the calculations in progressively larger domains. Fig. 8 displays the schematic diagram of the hexagonal pattern and the domain with aspect ratio 1:3 for the one-period hexagonal pattern.

Figure 8: Schematic diagram of (a) the hexagonal pattern and (b) the domain with aspect ratio 1:3 for the one-period hexagonal pattern

At time nΔt, the discrete total energy of the numerical solution for the NCH equation is defined by

ℰd(ϕn)=∑m=1Nx∑n=1Ny(h2F(ϕijn)+ε22[(ϕi+1,jn−ϕijn)2+(ϕi,j+1n−ϕijn)2]+α2[(ψi+1,jn−ψijn)2+(ψi,j+1n−ψijn)2]).

We perform the numerical test for the discrete total energy. For numerical simulation, Ω=(0,Rx=hxNx)×(0,Ry=3Ry) is used, where hx is spatial step size for x-axis. The parameters used hx=0.0015, Nx=512, α=100, Δt=0.125, ε=1/(202). The initial condition is ϕ(x,y,0)=−0.3+0.2rand(x,y) on Ω. We perform simulations until the computational solution becomes a numerical equilibrium solution. Fig. 9 shows the discrete total energy dissipation and the snapshot of the numerical solutions at the time indicated by the dot.

Figure 9: The discrete total energy with a snapshot of the numerical solutions

We observe the formation of a hexagonal pattern in the discrete equilibrium state and discrete total energy dissipation from numerical results. Chen et al. [56] presented a hydrodynamically-coupled phase-field model for diblock copolymer melts based on a conservative Allen–Cahn equation that preserves the volume fraction of the two monomers. In addition, the authors developed the linear and second-order time-marching method for the presented phase-field model. This method is easy to implement and can also be applied to a variety of phase-field models, such as CH equation for the diblock copolymer melts.

2.3 Alternating Direction Explicit

Yang et al. [57] developed an explicit FDM for the Ohta–Kawasaki model to describe microphase separation patterns in diblock copolymer melts. Their approach employs a Saul’yev-type scheme, which is grounded in a linearly stabilized convex splitting method, to achieve effective discretizations of the model equations. This method enhances the numerical stability of the simulations, allowing for more stable predictions of the complex behavior exhibited by diblock copolymer melts. The discrete domain is defined by using cell center points. The NCH Eqs. (1) and (2) are rewritten by applying the linear convex splitting-type scheme [58] as

ϕijn+1−ϕijnΔt=Δ[(ϕijn)3−3ϕijn]+2Δϕijn+1−ε2Δ2ϕijn+1−α(ϕijn+1−ϕ¯),(23)

where

Δϕij=ϕi−1,j+ϕi+1,j+ϕi,j−1+ϕi,j+1−4ϕijh2 and Δ2ϕij=Δϕi−1,j+Δϕi+1,j+Δϕi,j−1+Δϕi,j+1−4Δϕijh2.

Then, we use the Saul’yev-type method [59]. There are a total of 8 cases considering the order of i and j, and whether i is increasing or decreasing, as well as whether j is increasing or decreasing. One of these cases is as follows, and the remaining 7 cases are derived similarly.

For j=1,2,…,Ny,  for i=1,2,…,Nx,(24)

ϕijn+1−ϕijnΔt=Δ[(ϕijn)3−3ϕijn]+2h2(ϕi−1,jn+1+ϕi+1,jn−2ϕijn−2ϕijn+1+ϕi,j−1n+1+ϕi,j+1n)−ε2h4[ϕi−2,jn+1+ϕi+2,jn+ϕi,j−2n+1+ϕi,j+2n+2(ϕi−1,j−1n+1+ϕi−1,j+1n+ϕi+1,j−1n+1+ϕi+1,j+1n)−8(ϕi−1,jn+1+ϕi+1,jn+ϕi,j−1n+1+ϕi,j+1n)+10ϕijn+10ϕijn+1]−α(ϕijn+1−ϕ¯).(25)

We can simplify Eq. (25) as

ϕijn+1=1r{ϕijnΔt+Δ[(ϕijn)3−3ϕijn]+2h2(ϕi−1,jn+1+ϕi+1,jn−2ϕijn+ϕi,j−1n+1+ϕi,j+1n)−ε2h4[ϕi−2,jn+1+ϕi+2,jn+ϕi,j−2n+1+ϕi,j+2n+2(ϕi−1,j−1n+1+ϕi−1,j+1n+ϕi+1,j−1n+1+ϕi+1,j+1n)−8(ϕi−1,jn+1+ϕi+1,jn+ϕi,j−1n+1+ϕi,j+1n)+10ϕijn]+αϕ¯},(26)

where r=1/Δt+4/h2+10ε2/h4+α.

We consider the irregular domain Ω^, which is depicted as the gray region in Fig. 10. The initial condition on the full domain Ω=(−2,2)×(−2,2) is given by

ϕ(x,y,0)={ϕ¯+0.15rand(x,y)if(x,y)∈Ω^,−0.75otherwise.

Figure 10: Schematic illustrations of the irregular domain Ω^

For the numerical test, we used parameters Nx=Ny=800, ε=0.01, α=800, ϕ¯=0, −0.3, and Δt=h2. Fig. 11 displays the temporal evolution of the numerical solution for the NCH equation on the complex domain using the alternating direction explicit FDM for the phase-field equation.

Figure 11: Pattern formations in irregular domain Ω^. From the top to bottom, ϕ¯=0 and −0.3

Next, we perform numerical simulations to investigate how different boundary conditions, such as Dirichlet and Neumann boundary conditions, influence the simulation results. The initial condition is given by ϕ(x,y,0)=rand(x,y) on Ω=(−1,1)×(−1,1). For zero Dirichlet boundary condition, the boundary values are set as ϕ−1,j=ϕ0,j=ϕi,−1=ϕi,0=ϕNx+1,j=ϕNx+2,j=ϕi,Ny+1=ϕi,Ny+2=0 for i=1,2,…,Nx, j=1,2,…,Ny. For non-zero Dirichlet boundary condition, the boundary values are assigned as ϕ−1,j=ϕ0,j=ϕi,−1=ϕi,0=ϕNx+1,j=ϕNx+2,j=ϕi,Ny+1=ϕi,Ny+2=−0.7 for i=1,2,…,Nx, j=1,2,…,Ny. For Neumann boundary condition, the boundary values are set as ϕ−1,j=ϕ2,j, ϕ0,j=ϕ1,j, ϕNx+1,j=ϕNx,j, and ϕNx+2,j=ϕNx−1,j for j=1,2,…,Ny, ϕi,−1=ϕi,2, ϕi,0=ϕi,1, ϕi,Nx+1=ϕi,Nx, and ϕi,Nx+2=ϕi,Nx−1 for i=1,2,…,Nx.

Fig. 12 shows the numerical solutions obtained under different boundary conditions. From top to bottom, the results correspond to ϕ¯=0 and ϕ¯=−0.3. The numerical solutions using Dirichlet boundary conditions with values of 0 and −0.7 are shown in Fig. 12a,b, respectively. The solution using the Neumann boundary condition is presented in Fig. 12c. We observed that the Dirichlet boundary condition enforces fixed values near the boundaries, influencing pattern formation. The non-zero Dirichlet boundary condition with a boundary value of −0.7 produces a different pattern near the boundary compared to the zero Dirichlet boundary condition. In contrast, the Neumann boundary condition maintains continuity at the boundaries, resulting in smoother and more regular patterns.

Figure 12: Numerical solutions at time t=1000Δt with different boundary conditions. (a) zero Dirichlet, (b) non-zero (−0.7) Dirichlet, and (c) Neumann boundary conditions. From the top to bottom, ϕ¯=0 and −0.3

2.4 Invariant Energy Quadratization Approach

The invariant energy quadratization (IEQ) method is a numerical method designed to ensure energy stability when solving gradient flow problems. Originally proposed by Yang [60], the IEQ method transforms nonlinear partial differential equations into a form that allows for constructing energy-stable time-marching schemes. We simply describe an IEQ approach for the NCH equation. The auxiliary variable is defined as q(ϕ)=F(ϕ)+C. Then, the total free energy is rewritten as

ℰ(q,ϕ)=∫Ω(ε22|∇ϕ|2+α2|∇ω|2+q2)dx−C.

We get

∂ϕ∂t=Δμ,(27)

μ=qH(ϕ)−ε2Δϕ+αω,(28)

dqdt=12H(ϕ)∂ϕ∂t,(29)

H(ϕ)=F′(ϕ)F(ϕ)+C.

Therefore, we can obtain the energy dissipation law of the system (27)–(29) by taking the L2 inner product of Eqs. (27)–(29) with −μ, ϕt, and −2q, respectively.

ddtℰ(q,ϕ)=−‖∇μ‖2≤0.

The first-order time discretization IEQ method for the system (27)–(29) is defined by

ϕn+1−ϕnΔt=Δμn+1μn+1=Sqn+1Hn−ε2Δϕn+1+ωn+1qn+1−qn==12Hn(ϕn+1−ϕn).

The solution algorithms for the CH equation using the IEQ scheme is studied by Chen et al. [61], which are unconditionally stable. Two numerical methods, each first-order and second-order are reported. Numerical simulations showed that a large time step can be adopted while maintaining the energy decrease, therefore verifying the applicability of the proposed method.

2.5 Scalar Auxiliary Variable Approach

The SAV scheme was originally proposed by Shen et al. [62] and based on the IEQ approach, constructs energy-stable and efficient time discretization methods for gradient flows. The SAV scheme is applicable to different gradient flows and can be extended to higher-order by applying the backward differentiation formula (BDF) and Adam–Bashforth methods. We describe an SAV scheme for the NCH equation. The scalar auxiliary variable r(t) is defined by

r(t)=∫ΩF(ϕ(x,t))dx+C,

where C is a non-negative constant, which guarantees the value beneath the square root is not zero. Thus, we have

ℰ(r,ϕ)=∫Ω(ε22|∇ϕ|2+α2|∇ω|2)dx+r2−C.