The nonlocal Cahn±Hilliard model [1,2] has found widespread use in various scientific fields, such as materials science, soft matter physics, and biophysics, to study complex systems exhibiting phase separations. Compared with the classical Cahn±Hilliard model [3], which is primarily a phenomenological model, the nonlocal Cahn±Hilliard model takes into account the effects of long-range forces such as van der Waals forces or electrostatic interactions, which are critical in nonuniform systems and significantly influence the dynamics of phase separation. In fact, a complete microscopic derivation of phase separation using the statistical mechanics approach leads to the nonlocal Cahn±Hilliard model [4] as opposed to the classical version. Therefore, from a physical perspective, when modeling two-phase flows with interfacial effects, coupling the nonlocal Cahn±Hilliard equation with the fluid dynamical equations to simulate a two-phase flow system is a more reasonable model than using the extensively employed local Cahn±Hilliard model [5±13] in the aspect of accounting for the long-range forces present in the system, thereby providing a more precise and physically consistent description of the fluid dynamics.
The primary objective of this article is to develop some efficient numerical schemes to solve the two-phase flow model resulting from the coupling of the flow field with the nonlocal Cahn±Hilliard equation, where the flow field comprises of two commonly-used regimes, the free flow (Navier±Stokes) and the Hele-Shaw flow (Darcy). Although theoretical studies related to topics such as global attractor, existence and uniqueness of strong solutions, regularities, and applications of these models (e.g., tumor models) exist, cf. [14±21], there has been a paucity of discussion regarding the development of numerical algorithms for them. The current focus of numerical method design has been primarily on the nonlocal Cahn±Hilliard equation itself, e.g., the convex-splitting method [22±24], the Invariant Energy Quadratization approach [25±29] or its various version of Scalar Auxiliary variable approach [30,31], linear stabilization method [32,33], etc. However, integrating the nonlocal Cahn±Hilliard equation with the flow field leads to a model with multiple nonlinear coupling terms, posing a daunting challenge in designing a trustworthy numerical method that meets critical requirements such as unconditional energy stability, linearity, high order temporal accuracy, and full decoupling. Notably, preserving the discrete energy dissipation law is of paramount importance when dealing with nonlinear dissipative systems. Such a pursuit eliminates spurious numerical solutions, while assuring the dependability and accuracy of the results. Furthermore, properties like linearity and full decoupling offer a host of advantages for enhanced ease and efficiency of computations. Consequently, maintaining secondorder accuracy in time, while also upholding these pivotal properties, is imperative for ensuring the precision of the numerical solution.
Hence, it is pertinent to pose the question of what obstacles are encountered in formulating numerical methodologies, as expounded earlier, for the coupled models comprising the nonlocal Cahn±Hilliard equation and either the Darcy or Navier±Stokes equations. A comprehensive investigation of the mathematical model structure uncovers numerous intricate aspects that demand meticulous attention during the discretization process, including: (i) the linear coupling between the velocity and pressure; (ii) the nonlinear coupling among the velocity and phasefield variable in the advection and surface tension; and (iii) the nonlinear convolutional term with nonlocal effects and the double-well potential with stiffness. Remarkably, numerous numerical methods have achieved success and established credibility for each of the three aspects. For instance, the projection-type approach is a well-known and efficient method for handling the linear coupling between velocity and pressure in (i) by solving two independent momentum equation and the pressure projection Poisson-type equation. A range of numerical methods exist for (iii), as mentioned earlier. In addition, various efficient approaches have been developed for (ii), such as fullyimplicit [9], semi-implicit [11], and stabilized explicit [6,7] methods. However, combining all three aspects in one model presents a formidable challenge, as the authors have noticed that no successful attempts have been undertaken to acquire numerical schemes that exhibit second-order accuracy, linearized, full decoupled, and unconditional energy stable. Consequently, there is currently a lack of fully discrete schemes available for the flow-coupled nonlocal Cahn±Hilliard model.
The objective of this research paper is to develop numerical schemes for solving the flow-coupled nonlocal Cahn±Hilliard model that exhibits the numerical advantages described above to effectively capture the dynamics and enable accurate simulations. To accomplish this goal, we propose to combine existing numerical algorithms, including the projection method for solving the fluid equations, the SAV method for linear discretization of the double-well potential, and the ªZero-Energy-Contributionº (ZEC) method for handling those coupled nonlinear advective and surface tension terms. The ZEC method has been widely used for constructing decoupled numerical schemes that involve the local Cahn±Hilliard equations and other equations, such as hydrodynamics, magnetic fields, or others (see [34±38]). Both the SAV and ZEC methods involve introducing additional nonlocal variables and constructing corresponding time-oriented ordinary differential equations (ODEs). In the SAV method, the nonlocal variable is obtained by forcing the nonlinear convolutional term quadratically, and its corresponding ODE is derived by subjecting the variable to a time derivative. The ZEC method, on the other hand, constructs the nonlocal variable by combining the inner products of coupling terms with certain variables. By utilizing these methods and discretizing the space using the Fourier-Spectral method, we can obtain a fully discrete numerical scheme that can effectively solve the coupled system by solving only several decoupled, linear elliptic equations with constant coefficients, thus allowing us to simulate the model efficiently, stably, and accurately.
The paper is organized as follows. Section 2 provides a brief introduction to the two coupled models consisting of the nonlocal Cahn±Hilliard equation coupled with Navier±Stokes/Darcy flow field. In Section 3, we develop two fully discrete Fourier-Spectral numerical schemes for solving these two models and establish their unconditional energy stability. We also provide a detailed implementation process. Section 4 presents several numerical examples to validate the effectiveness of the proposed schemes. Finally, we give some concluding remarks in Section 5.
In this section, we provide a brief introduction to the two-phase flow models that consist of the nonlocal Cahn± Hilliard equation coupled with the free flow (Navier±Stokes) model and the Hele-Shaw flow (Darcy) model. We refer to these models as the nCHNS model and the nCHD model, respectively. Both models have been extensively used to simulate the two-phase flow systems, with the nCHD model specifically designed to capture fluid dynamics in porous media.
We assume that Ω ∈ R d (d = 2, 3) denotes a smooth, bounded, and connected domain. For any functions u(x), v(x) ∈ L 2 (Ω ), we denote the L 2 inner product as (u, v) = ∫ Ω u(x)v(x)d x, and the L 2 norm as ∥u∥ = √ (u, u).
The phase-field function φ(x) denotes the mass (or volume) fraction of two different fluid components in the fluid mixture, i.e.,
The two values are connected by a smooth thin layer. The nonlocal operator L reads as (cf. [14,15]):
where J (x) is the interaction kernel that is Ω -periodic and satisfies the following properties [14,15], is described by the following equations:
where M > 0 denotes the mobility constant, λ > 0 denotes the surface tension parameter, ν > 0 denotes the fluid viscosity, p denotes the pressure, u denotes the fluid velocity field, f (φ) = F ′ (φ) = 1 ϵ 2 (φ 3φ), ϵ is related to the width of the diffusive interface, µ denotes the nonlocal chemical potential. We assume that all of the above variables satisfy the periodic boundary conditions. The initial conditions are set as follows:
The system (2.1a)±(2.1d) holds an energy dissipation property. To obtain it, we need two important identities, which also play an important role in deriving our analysis, read as
Then, by taking the L 2 inner products of (2.1a)±(2.1c) with µ, -φ t and u, respectively, and then using integration by parts and (2.2a)±(2.2b), we obtain
where
(2.4)
The nCHD model shares almost the same structure as the nCHNS model expect that the fluid momentum equation is replaced with the Darcy's law (see [20,21]):
)
where u denotes the dimensionless seepage velocity, τ is a positive parameter, α is the dimensionless hydraulic conductivity. We assume that the boundary conditions of all the above variables are periodic too. The initial conditions are set as u| t=0 = u 0 (x), p| t=0 = p 0 (x), φ| t=0 = φ 0 (x).
Using the similar way, it can be derived that the nCHD system (2.5a)±(2.5d) follows the energy dissipation property as
where
(2.7)
The Eqs. (2.2a)±(2.2b), (2.4), and (2.7) reveal an important feature of the nonlocal Cahn±Hilliard Navier± Stokes/Darcy system. Specifically, the two coupled terms φ∇µ, ∇ • (uφ) and the term (u • ∇)u do not contribute to the energy dissipation of the system, which is the so-called ªzero-energy-contributionº property as described in [34±38]. This unique property will be exploited in the construction of our numerical scheme.
To solve these two models, we discretize the space by using the Fourier-Spectral method [32,39], and integrate several efficient temporal discretization methods, including the projection method for the flow field, the SAV method [40,41] for linearizing the nonlinear term f (φ), and the ZEC approach to obtain the full decoupling structure [34±38, 42,43], to form fully discrete schemes.
We first introduce some basic concepts and notations for the Fourier spectral method. Let δt > 0 be a time step size and set t n = nδt for 0 ≤ n ≤ N T with T = N T δt. We assume Ω = (-L , L) 3 and N denotes an even number. Let h = 2L/N be the spatial step size. Define that
All periodic grid functions are denoted by M h , namely,
We then define the discrete inner product of two functions f and g in M h as
and define the associated L 2 norm of f as
where f i jk , g i jk represent the values of f, g taken at the grid point (x i , y j , z k ). For a function f ∈ M h , the discrete Fourier transform f = P f is defined by
The function f can be reconstructed by using the inverse Fourier transform f = P -1 f given by
flmn . The Fourier spectral approximations to the partial derivatives can be expressed as [39] D x = P -1
Dx P, D y = P -1 Dy P, D z = P -1 Dz P,
the discrete gradient, divergence and Laplace operators are computed by
The discrete nonlocal diffusion operator L h on any grid function f ∈ M h through the use of L h = P -1 Lh P by (see [32]):
where λ lmn are the eigenvalues of L h given by
.
Furthermore, we have the following discrete integration by parts formulas (cf. [32]) that read as:
To apply the SAV method to linearize the nonlinear term f (φ), we need introduce a nonlocal variable U (t) as
where the constant C 0 > 0. Then, we can use the variable U to reformulate the chemical potential µ as:
where
. In fact, by integrating the second equation in (3.5) over t and given the initial
we can obtain (3.4). This allows us to recast the equation of the chemical potential µ in (2.1b) and (2.5b) into the equivalent reformulation in (3.5). Therefore, it can be concluded that (3.5) is simply an alternative representation of the chemical potential µ in (2.1b) and (2.5b).
To apply the ZEC approach to get the decoupled structure, we shall introduce another nonlocal variable Q which satisfies the following ODE system as
From (2.2a)±(2.2b), we can see that the time derivative of Q is zero, i.e., Q t = 0. Furthermore, from the initial condition Q| t=0 = 1, we can conclude that Q remains constant at 1 for all time.
Thanks to (3.5) and (3.6), we can rewrite the nCHNS system (2.1a)±(2.1d) as
)
)
with the following initial conditions
Remark 3.1. Note that we multiply Q to the terms in Eq. (3.7d). However, by using the fact that Q ≡ 1, we know this modification will not affect the whole system. It is worth noting that the rationale behind the persistent modifications to the PDE lies in the realization that although these alterations may be equivalent at the PDE level ± that is, in the continuous setting ± they are essential in facilitating the development of the algorithms when dealing with discretizations. As such, one could assert that the modified PDE is algorithmically-favorable.
We can show that the equivalent system (3.7a)±(3.7f) holds the energy dissipation law.
Lemma 3.1. The system (3.7a)±(3.7f) holds the following energy dissipation property
where
Proof. By taking the L 2 inner product of (3.7a)±(3.7b) with µ and -φ t , respectively, we obtain
By multiplying (3.7c), (3.7f) with 2λU and Q, respectively, we have
By taking the L 2 inner product of (3.7d) with u and using (3.7e), we arrive at
By summing up the above five equations together, we obtain (3.8). □
Now we are able to develop the fully discrete numerical scheme as follows.
For n ≥ 1, assuming that (
∈ M h by solving the following discrete scheme:
where S > 0 is a stabilization parameter,
We show the fully discrete scheme (3.9a)±(3.9g) holds the energy stability unconditionally, as follows. Note the following two identities will be used in the proof:
Theorem 3.1. The scheme (3.9a)±(3.9g) satisfies the discrete energy stability as follows
where
Proof. Taking the L 2 inner product of (3.9e) with 2δt ũn+1 h , we obtain
By taking the L 2 inner product of (3.9f) with u n+1 h , we have
By using (3.13) and (3.10b), we deduce
We rewrite (3.9f) as 3 2δt
Square both sides of (3.15), and multiply the obtained equation with 2δt 2 3 , we derive
Further, by taking the L 2 inner product of (3.9f) with 2δtu n+1 h , we get 3 2
Taking the L 2 inner product of (3.9a) and (3.9b) with 2δtµ n+1 h and -(3φ
, respectively, we derive
and
By multiplying (3.9c) and (3.9d) by 2λU n+1 h and 2δt Q n+1 h , and using (3.10b), we obtain
and 1 2
Hence, by combining (3.18)±(3.22), and using (3.10a), we arrive at 1 2
, which completes the proof. □ Remark 3.2. From Theorem 3.1, we know that the fully discrete numerical scheme follows the modified discrete energy dissipation property. However, in view of numerical approximation, we find the proposed scheme preserve energy dissipation property of the original system. Remark 3.3. Despite its appearance as a coupled scheme rather than a full decoupling scheme, the fully discrete scheme outlined in Eqs. (3.9a)±(3.9g) can be effectively implemented with decoupled type computations. This can be achieved by leveraging the nonlocal nature of the variable Q, which enables the implementation of decoupled computations in the overall computation process, shown as follows.
In the subsequent five steps, we will elaborate on the fully-decoupled implementation process in greater detail.
Step 1:
(3.23)
By substituting (3.23) into (3.9a)±(3.9c), we obtain
Hence, according to Q n+1 h , we decompose the above system into the following two systems
We use the nonlocal variable U n+1 h1 and U n+1 h2 to continue to split the variables as:
Here, we assume that the boundary conditions for φ n+1 h11 , φ n+1 h12 , φ n+1 h21 , φ n+1 h22 , µ n+1 h11 , µ n+1 h12 , µ n+1 h21 , and µ n+1 h22 are all periodic.
As a result, the values of φ n+1 hi j and µ n+1 hi j for i, j = 1, 2 can be easily obtained, given that Eqs. (3.28)±(3.31) represent linear systems. Furthermore, the values of (φ n+1 h12 , µ n+1 h12 ) and (φ n+1 h22 , µ n+1 h22 ) satisfy the same system, which implies that φ n+1 h12 = φ n+1 h22 and µ n+1 h12 = µ n+1 h22 .
Step 2: We compute U n+1 h1 , U n+1 h2 . By using (3.9c) and (3.23), we arrive at
Together with (3.27), we have
Assuming that 1 -1 2 ⟨H n, * , φ n+1 h12 ⟩ ̸ = 0 and 1 -
⟨H n, * , φ n+1 h22 ⟩ ̸ = 0 (shown below), we get
.
We now prove the assumption. By taking the L 2 inner product of (3.29) with -2Mδt 3 µ n+1 h12 and φ n+1 h12 , respectively, we have
Noticed that φ n+1 h12 = φ n+1 h22 . Hence, we have 1 -
Here, we assume that the boundary conditions of ( ũn+1 h1 , ũn+1 h2 ) are periodic. One can easily obtain ( ũn+1 h1 , ũn+1 h2 ), due to the fact that (3.34) are constant-coefficient elliptic systems.
Step 4: We compute Q n+1 h . By using (3.9d), (3.23) and (3.33), we derive
Assuming that 3 2δt
η 2 ̸ = 0 (shown below), we have
.
We now prove the above assumption. By taking the L 2 inner product of (3.34) with (0, ũn+1 h2 ), we have
By taking the L 2 inner product of (3.26) with µ n+1 h2 and -3 2Mδt φ n+1 h2 , respectively, we derive
where the last inequality applies the second equation in (3.32). According to the above two equations, we know that -η 2 ≥ 0, i.e., 3 2δt
-
h by using (3.23) and (3.33).
Step 5: We compute p n+1 h and u n+1 h . By using (3.9f) and (3.9g), we obtain
h . Finally, we can get u n+1 h by using
). As evident from the aforementioned five steps, it is apparent that all the equations requiring a solution contain an independent variable. Hence, our numerical scheme is deemed fully decoupled, thus rendering it highly efficient.
We now attempt to construct an efficient fully discrete scheme for the nCHD system. Drawing upon a similar approach employed in the nCHNS model, we define a nonlocal variable Q and an associated ODE system, presented below:
Through the application of identity (2.2a), it can be deduced that (3.35) is equivalent to Q t = 0, Q| t=0 = 1 with the solution of Q(t) = 1. By combining (3.35) with the reformulation of the chemical potential µ, as shown in (3.5), and by multiplying the corresponding terms by the variable Q, we can rewrite the nCHD system (2.5a)±(2.5d) as
)
where
, and the system is equipped with the periodic boundary conditions and the following initial conditions u| t=0 = u 0 (x), p| t=0 = p 0 (x),
The system adheres to the energy dissipation law, as its process closely resembles that of Lemma 3.1. Therefore, we will skip over the process of deriving the energy law of the PDE and proceed directly to the fully discrete Fourier-Spectral scheme for the system (3.36a)±(3.36f), that reads as follows.
For n ≥ 1, assuming that (
)
where S > 0 is a stabilization parameter, and
). We now show that the fully discrete scheme (3.37a)±(3.37g) is unconditionally energy stable.
Theorem 3.2. The scheme (3.37a)±(3.37g) satisfies the discrete energy stability as follows 1
δt
where
Proof. By taking the L 2 inner product of (3.37e) with 2δt ũn+1 h , we obtain
By taking the L 2 inner product of (3.37f) with u n+1 h , we obtain
By taking the square of (3.37f) and using (3.37g), (3.3), we have
Further, by taking the L 2 inner product of (3.37f) with 2δtu n+1 h , we get 3τ 2
By taking the L 2 inner product of (3.37a)
Multiplying (3.37c), (3.37d) by 2λU n+1 h and 2δt Q n+1 h and then using identity (3.10b), we obtain
Hence, by combining (3.43)±(3.47), we arrive at
It complete the proof of the theorem. □ Next, we present a detailed implementation process for the scheme (3.37a)±(3.37g). It should be noted that the equations in this scheme are nearly identical to (3.9a)±(3.9g), with the exception of (3.37d) and (3.37e). Therefore, the computation of φ n+1 h , µ n+1 h , and U n+1 h can be computed in the same way as described in Steps 1 and 2 of the previous section. Now we turn to the computed scheme for the variables ũn+1
By substituting (3.48) into (3.37e) and decomposing the obtained equation according to
h , with the periodic boundary conditions. The above equations are very straight forward to be solved.
Step 4: We compute Q n+1 h . By using (3.23), (3.48) and (3.37d), we arrive at
where
We can obtain 3 2δtθ 2 ̸ = 0 by using the similar way as the scheme for the nCHNS model. By combining the above two steps 3 and 4 with the steps 1, 2, 5 for the nCHNS model, the solution to the scheme (3.37a)±(3.37g) is obtained using full decoupling type computations.
In this section, we aim to demonstrate the accuracy and energy stability of the suggested schemes through various numerical examples, including accuracy/stability tests, as well as simulations of spinodal decomposition, fingering instability, and bubble rising. These simulations enable us to evaluate the effectiveness of the proposed schemes in accurately and consistently capturing physical phenomena.
In this example, we set the computed domain as Ω = [-π, π] 2 , and set the initial conditions as φ 0 (x, y) = cos(2x) cos(2y), u 0 (x, y) = (0, 0) T , p 0 (x, y) = 0. To evaluate the convergence order in time, we employ a spatial grid with sufficiently fine resolution (consisting of 128 2 Fourier modes) to ensure that the spatial error is insignificant compared to the temporal error. We continue to verify the energy stability of the developed schemes. The computational domain is set as Ω = [-π, π] 2 and 128 2 Fourier modes are used to discretize the space. The initial conditions are set as follows:
), u 0 (x, y) = (0, 0) T , p 0 (x, y) = 0, where r 1 = 1.4, r 2 = 0.5, x 1 = -0.8, x 2 = 1.7. The kernel J (x, y) is defined as J (x, y) = 1 ϵ 4 exp(-x 2 +y 2 ϵ 2 ). For the nCHNS model, we set the parameters as ν = 100, λ = 0.01, ϵ = 0.05, M = 1, C 0 = 10, and S = 2, while for the nCHD model, we set the parameters as α = 100, τ = 1, λ = 0.01, ϵ = 0.05, M = 1, C 0 = 10, and S = 2. We plot the computed modified discrete energy (3.11) by varying time steps for the nCHNS model in Fig. 4.2(a). The obtained energy curves exhibit a monotonic decay, providing evidence for the energy decay property. To further illustrate the accuracy of the modified discrete energy, we compare the temporal evolution curves of the original energy (2.4) and the discrete energy (3.11) in Fig. 4.2(b), using δt = 0.01/2 3 . The overlap of these two energy curves supports the consistency between these two energies. Similarly, for the nCHD model, we carry out similar numerical tests, and the results are shown in Fig. 4.3(a) and (b). As the phenomena observed are exactly similar to those of the nCHNS model, we refrain from going into further detail.
In this numerical test, we perform a spinodal decomposition, which is a benchmark phenomenon of the Cahn± Hilliard equation. This phenomenon refers to the spontaneous separation of a homogeneous mixture of two or more substances into distinct regions or domains with different compositions. Spinodal decomposition is a critical research area in materials science, as it is a common mechanism in the formation of microstructures in various systems, including alloys, polymers, and liquid crystals, and the resulting microstructure can significantly affect the material's physical and chemical properties. 2 Fourier modes for 2D case, and 128 3 Fourier modes for 3D case. For the 2D case, we plot the profiles of φ at various times for the initial condition with φ0 = 0.5 in Fig. 4.4. The final equilibrium solution exhibits a circular phase due to the coarsening effects. Similarly, Fig. 4.5 shows the profiles of 3D isosurface of {φ = 0}, which also display a final circular equilibrium state. In Fig. 4.6, we plot the time evolution of the total free energy functional (3.38) for both 2D and 3D simulations, confirming the energy decaying property of our scheme.
When a liquid droplet is placed in a rotating Hele-Shaw cell, the resulting centrifugal force induces instability and deformation in the fluid interface. This instability leads to the development of finger-like protrusions, a common phenomenon observed in various fluid systems. The competition between surface tension, which tends to smoothen the interface, and the interfacial energy associated with the contact line, which drives the formation of fingers, characterizes this instability. The centrifugal force generated by the cell's rotation perturbs the fluid interface, leading to the formation of fingers. Over time, these fingers undergo growth in size and number, and their dynamics are governed by a combination of viscous and capillary forces.
In this numerical example, our focus is to examine the fingering pattern instability that arises in this scenario. This is achieved by incorporating the effect of the rotational force into the Darcy equation. Specifically, we replace the momentum equation (2.5c) with a modified version that accounts for the rotational force, namely, where e z = (0, 0, 1), and r = (x, y) for 2D and r = (x, y, 0) for 3D. The kernel J (x) is taken to be J
We also carry out 2D and 3D simulations respectively, where space is discretized by using 512 2 Fourier modes for a 2D domain of Ω = [-π, π] 2 , and 128 3 Fourier modes for 3D domain of Ω = [-π, π] 3 . The initial conditions of 2D simulation read as
), u 0 (x, y) = (0, 0) T , p 0 (x, y) = 0, where r = 1.3 + 0.01rand(x, y). The initial conditions of 3D simulation read as
), u 0 (x, y, z) = (0, 0, 0) T , p 0 (x, y, z) = 0, where r = 1 + 0.01rand(x, y, z). The model parameters are set as τ = 1, α = 400, M = 0.01, ϵ = 0.015, C 0 = 10, S = 2, ω = 5, g = 5.
In Fig. 4.7, we plot snapshots of the phase-field variable φ at different times where the surface tension λ = 1e-5, depicting the formation and evolution of finger-like structures on the fluid interface in a rotating Hele-Shaw cell. Initially, the circular droplet of fluid owns a smooth and continuous interface. As the rotation of the Hele-Shaw cell is initiated, the centrifugal forces generate perturbations on the fluid interface, leading to the formation of initial finger-like protrusions. These initial protrusions undergo significant growth in size and number, eventually evolving into a highly complex pattern of finger-like structures that exhibit highly non-uniform shapes and sizes, with some fingers growing faster than others. In Fig. 4.8, we conduct a series of numerical simulations with varying surface tension parameters λ, which demonstrate that the surface tension plays a crucial role in determining the number and shape of the finger-like structures formed on the fluid interface. Specifically, as surface tension decreases, we observe an increase in the number of finger-like structures formed. The 3D numerical simulations, presented in Fig. 4.9, also demonstrate the formation of finger-like structures over time, with the isosurfaces of {φ = 0} showing the evolution of these structures. Varying the surface tension parameter λ in Fig. 4.10 directly affects the number of fingers formed, with smaller values of λ leading to the formation of a larger number of fingers, consistent with the 2D simulations. These findings are also qualitatively consistent with the results obtained from the numerical simulations in [43] and the experimental results in [44,45].
The Rayleigh instability is a well-known hydrodynamic instability that occurs when a liquid thread or column is subjected to external perturbations such as gravity, surface tension, and viscosity. When the magnitude of the perturbation is small, the instability initially manifests as small ripples on the surface of the liquid thread. However, as time progresses, the instability causes the formation of bulges or beads along the length of the thread. Eventually, the bulges become so large that they break up into smaller droplets, which results in the fragmentation of the liquid thread.
In this example, we simulate the Rayleigh instability using the developed scheme (3.37a)±(3.37g) to solve the nCHNS model. The computational domain is set as Ω = [-π, π]×[-π/2, π/2]×[-π/2, π/2]. The initial condition is set as a liquid thread with a small perturbation on the interface, which is given by φ 0 (x, y, z) = tanh( 0.4 -√ y 2 + z 2 + 0.05 cos(2x) √ 2ϵ ), u 0 (x, y, z) = (0, 0, 0) T , p 0 (x, y, z) = 0.
In this example, we simulate the dynamics of a rising bubble under the influence of a gravity field using the nCHNS model and the numerical scheme (3.9a)±(3.9g). The gravity force is approximated using the Boussinesq approximation, which involves equipping the Navier±Stokes equation with a gravity field given by
where ρ 1 is the fluid density inside the bubble and ρ 2 is the density of the surrounding fluid medium, g is the gravitational acceleration with g = (0, -g 0 ) T for 2D and g = (0, 0, -g 0 ) T for 3D. The convolution kernel J (x) is taken as J (x) = 1 ϵ d+2 exp(- ), u 0 (x, y, z) = (0, 0, 0) T , p 0 (x, y, z) = 0.
In Fig. 4.13, we plot the isosurfaces of {φ = 0} at different times. Upon observation, it is evident that the dynamic movements of a single bubble rising in 3D are comparable to those observed in 2D environment (Fig. 4.12(a)).
As the droplet rises due to the force of gravity, it undergoes significant deformation, transforming from a spherical shape into a cap shape.
This paper presents efficient numerical schemes for a coupled system consisting of the nonlocal Cahn±Hilliard model with two flow regimes, the Navier±Stokes equation's free flow, and the Darcy flow's Hele-Shaw cell. The proposed scheme offers several desirable numerical features, such as second-order accuracy in time, linearity, unconditional energy stability, and a fully decoupled structure. The scheme combines the projection method, the SAV method, the Fourier spectral method, and the ZEC decoupling technique, leveraging the inherent properties of the coupled terms. The linearity and fully decoupled structure of the proposed scheme enable practical implementation efficiency. Rigorous proofs are provided for the scheme's unconditional energy stability and practical implementation process, and several 2D and 3D numerical examples are performed to demonstrate the scheme's accuracy and effectiveness.