Surfactants are introduced into heterogeneous solutions to modify or decrease the surface tension, thereby stabilizing the mixture and preventing phase separation. By lowering surface tension, surfactants enhance the miscibility of otherwise immiscible fluids [1], allowing for more uniform mixing. This is particularly useful in industrial processes such as emulsification, where stable interfaces between fluid phases, like oil and water, are essential for producing stable products such as detergents, paints, and pharmaceuticals. Since Laradji's seminal work decades ago [2,3], the phase-field method, also known as the diffusive interface approach, has become scheme, providing rigorous proof of unconditional energy stability, demonstrating solvability, and detailing the fully-decoupled type implementations. In Section 4, we present several numerical examples to demonstrate the accuracy and efficiency of the proposed method, with particular emphasis on the simulation of Saffman-Taylor fingering patterns on various complex topological surfaces. Finally, concluding remarks are provided in Section 5.
In this section, we present the binary phase-field surfactant model coupled with the Darcy flow system in a Hele-Shaw cell on a closed surface. To accurately simulate the underlying physical mechanisms of this fluid system, it is essential to implement an effective parameterization method on the surface. Therefore, we first provide a brief overview of the relevant differential geometry operators, followed by a detailed presentation of the governing system.
This subsection introduces some fundamental concepts related to surface parameterization. For any given surface domain, an appropriate parameterization is essential for mathematical modeling. Let us consider a smooth, orientable surface ā ∶= {Ė( 1 , 2 ) * R 3 ∶ ( 1 , 2 ) * D ⊂ R 2 }, where 1 and 2 represent the surface parameters. Using this parameterization, the surface derivatives can be expressed through the following general formulas:
The coefficients of the first fundamental form for a surface are given by = (Ė , Ė ) with , = 1, 2, where (ç, ç) represents the inner product. The important surface normal is given by Č = Ė 1 ×Ė 2 ‖Ė 1 ×Ė 2 ‖ , where ç × ç and ‖ ç ‖ stand for the cross product and the Euclidean norm in R 3 , respectively.
We also denote [ ] = [ ] -1 , which represents the inverse of the matrix formed by the coefficients of the first fundamental form, and = det [ ], the determinant of the first fundamental form. These notations are essential for establishing the following geometric framework.
• Tangential gradient operator. On the function space 1 (ā), for any * 1 (ā), the tangential gradient operator ∇ is defined as
• Divergence operator. On the smooth vector field Ĕ * [ 1 (ā)] 3 , the divergence operator ∇ ç is defined as
(2.3)
• Laplace-Beltrami operator. For any * 2 (ā), the Laplace-Beltrami operator (LBO) is defined as
• Sobolev space on surface. For a given constant and a function * @ (ā), we denote ∇ the th order covariant derivative of function , with the convention ∇ 0 = . Let
and ā be a compact surface with at least th order smoothness Sobolev space (ā), which is the completion of (ā) in the sense of norm
.
(2.5)
For any functions , * 1 (ā), two classical inner products on the surface ā are described as ( , ) = + ā dĖ, and (∇ , ∇ ) = + ā ∇ ç ∇ dĖ.
In this subsection, we describe the Darcy flow-coupled binary phase-field surfactant model on a closed surface. The phasefield model utilizes two scalar variables, (Ė, ) and (Ė, ), to capture the dynamic behavior of the system. The scalar variable represents the local density variation between the two fluids, where = 1 corresponds regions dominated by fluid I, and = -1 denotes areas where fluid II (the ambient medium) prevails. To avoid discontinuity of , we introduce a smooth transition layer Computer Methods in Applied Mechanics and Engineering 436 (2025) 117733 with width ( ), which connects these distinct regions. Consequently, the fluid interface is characterized through the zero-level set = {Ė ∶ (Ė, ) = 0}. The second scalar variable, , denotes the local concentration of the surfactant within the binary system. For the binary surfactant system, the free energy ( , ) on the surface domain ā, coupling the two variables and , is given as (cf. [13,19,32,86,87]):
in which, the functions ( ), ( ) and ( , ) read as
2 2 + 4 |∇ | 4 , (2.7) with parameters ( 1 , 2 , 3 , 4 ) = (3.62, -7.25, 7.30, -3.68). The parameters of 1 , 2 , , , , 1 , 2 , are all positive, controlling various aspects of the free energy and model dynamics.
We provide a brief explanation about the total free energy ( , ) as follows. The total free energy is composed of four primary components: the bulk potentials for both and , the -coupling potential (the first two terms in ( , )), and the quartic gradient potential (the last term in ( , )). The bulk energy of consists of two terms: the gradient potential |∇ | 2 , which accounts for the hydrophilic tendency of the two fluid components (represented by = ±1), and the double-well potential ( ), which models their hydrophobic behavior. The bulk energy for the concentration variable has a similar structure, but with the quartic polynomial potential ( ) in another different formation (see Remark 2.1 for further discussion on ( )), where acts as the penalty parameter. The coupling potential ( , ) governs the interaction between and , with the term 1 controlling the surfactant's effect on the fluid interface, and the 2 term guaranteeing a low concentration of surfactant in the bulk phase. The quartic gradient potential term is included to keep the total free energy bounded from below, with the parameter chosen to be small ( ≪ 1), as seen in [32,87].
Remark 2.1. As demonstrated in numerous existing works, the choice of the hydrophobic functional, ( ), is quite flexible. Several formulations have been proposed for it thus far. The earliest version utilizes a quadratic polynomial, as seen in [13,32,86], while a fourth-order double-well potential is introduced in [6]. Additionally, the Flory-Huggins logarithmic form, expressed as ln + (1 -) ln(1 -), has been adopted in [12,13]. In this paper, we adopt the quartic polynomial potential for ( ) given in (2.7), as provided in [87]. This formulation offers a suitable approximation of the Flory-Huggins logarithmic potential while mitigating numerical instabilities, such as overflow issues, commonly encountered with logarithmic functions.
Assuming the fluid motion adheres to the mechanical principles of porous media and follows the generalized Fick's law, where the mass flux is proportional to the chemical potential gradient, the binary surfactant system coupled with Darcy flow on a surface given by:
)
where
where 1 , 2 are two mobility parameters for and , respectively, ē represents the dimensionless seepage velocity field, is the dimensionless hydraulic conductivity, and is the pressure. The fluid viscosity ( ) is modeled as a linear interpolation between the viscosities of the two fluids, 1 and 2 , with ( ) = 1 2 1 (1 + ) + 1 2 2 (1 -) where 1 and 2 correspond to the viscosities of fluid I ( = 1) and fluid II ( = -1), respectively. It should be noted that the temporal derivative of the seepage velocity ē is retained for flows in porous media, as also discussed in [61,[88][89][90].
We show the energy dissipation law followed by the above system (2.8)-(2.13) on the surface, as follows.
Theorem 2.1. The system (2.8)-(2.13) on the closed surface holds the law of energy dissipation:
) Computer Methods in Applied Mechanics and Engineering 436 (2025) 117733 where ( , , ē) = + ā ( 2 |ē| 2 + 1 ( 1 2 |∇ | 2 + 1 2 ( )) + 2 ( 2 |∇ | 2 + 1 2 ( )) + ( , ) ) dĖ. (2.16)
Proof. By taking the 2 -inner products of (2.8) with and applying integration by parts, we get
(2.17)
By taking the 2 -inner product (2.9) withand applying integration by parts, we obtain
By taking the 2 -inner products of (2.10) with and applying integration by parts, we get
By taking the 2 -inner product (2.11) withand applying integration by parts, we obtain
By taking the 2 -inner products of (2.12) with ē, and applying the divergence free condition (2.13) and integration by parts, we obtain
By combining the above five equalities (2.17)-(2.21), the energy dissipation law (2.15) is obtained, thus completing the proof. ¦ Note that, in deriving the energy law (2.15), we make use of the following two identities:
and
which are obtained after applying integration by parts. These identities reveal that both the surface tension and nonlinear advection terms do not contribute to the total energy of the system, thereby satisfying the ''zero-energy-contribution'' property. When coupled nonlinear terms exhibit this property, the ZEC decoupling method [20,83,84] can be effectively employed to construct decoupled-type algorithms, which will be described in the following section.
In this section, we construct the spatiotemporal full-discrete scheme for solving the surfactant-Darcy coupled system (2.8)-(2.13). The domain of interest consists of a surface with arbitrary topologies and geometric shapes, making the parameterization of the associated differential geometric operators highly challenging. To address this challenge, we adopt the IGA framework, which integrates CAD shape parameterization techniques with FEM discretization methods. Specifically, we utilize the recently developed subdivision-based IGA approach [17,18] for spatial discretization, thereby providing a geometrically accurate foundation for our scheme. Building upon this framework, a decoupled, full-discrete, second-order accurate in time, energy-stable algorithm is derived.
A surface domain is discretized into a finite element mesh, referred to as a subdivision element mesh in our method. Integral calculations must be performed within each subdivision element. Due to the complexity of the surface we are working with, its discretized form often involves unstructured grids with arbitrary topological connectivity. Loop subdivision is a generalization of quartic box splines defined on triangular meshes, offering the advantage of compatibility with any topological connectivity. To address spatial discretization, we perform numerical integration for the integral terms arising from the weak form of the PDE system. This integration is achieved by connecting a non-overlapping tiling from the parameter element to the physical domain of interest. It is important to note that integral computations on the subdivision element mesh are nontrivial, as they are not strictly elements in the classical finite element sense.
We first recall the computation of classical quartic box spline element. Consider the local topology of a regular patch, where each of its three control vertices has six incident edges. This regular patch can be exactly calculated using the following analytical expression Ė( , ) = 12 1 =1 ( , )Ė , (3.1)
where ( , , 1 --) are the barycentric coordinates of the unit parameter triangle, and are the 12 quartic box-spline basis functions associated with their corresponding 2-ring neighboring control vertices Ė (see Fig. 3.1(a) and details in [43]).
When addressing the local topology of an irregular triangle patch, a direct exact expression is not feasible because at least one of its three control vertices lacks six incident edges. Fortunately, we can apply the hierarchical strategy of quartic box spline expressions developed by Stam [43]. It requires that the considered triangle patch has only one control vertex with a valence other than six, which can be achieved with one-time initial refinement. In each refinement step, a triangle is divided into four sub-triangles, with the subdivision stencils providing topological refinement for both control vertices and control edges. The global mesh can then be refined repeatedly, allowing parameter values of interest to be restricted within a regular patch topologically. Ultimately, this repeated subdivision of the irregular patch produces a sequence of regular sub-patches spanned over the following three local domains:
where the necessary subdivision times = +1 -log( + )∕ log 2,. The system of hierarchical sub-triangles can be projected onto the unit triangle using the shift
Hence, the physical domain is defined by its hierarchical restriction to every triangular element:
where three sets of twelve control vertices Ė , are chosen from the Next, we use the stencils to describe the subdivision rules in Fig. 3.2, which essentially represent topological refinement. Loop subdivision is an approximation algorithm, where both of old and new vertex positions must be calculated at each refinement marching. After one subdivision, two types of control vertices are presented on the newly generated mesh. Among them, the control points inherited from the previous mesh need to have their positions updated according to the formula:
- [ 5 8 -( 3 8 + 1 4 cos 2
) 2 ] and is its valence; and (b):
subdivision stencil of an edge.
The IGA framework in our method is based on the Loop subdivision, which owns the hierarchical property and is the generalization of quartic box splines. In the implementation, we use adaptive refinement to achieve Gaussian integration of elements on the coarse mesh. The integrations for computing the matrix elements are computed by a 12-point Gaussian quadrature rule. That is, each domain triangle is subdivided into four sub-triangles and a four-point Gaussian quadrature rule is employed on each of the sub-triangles. Evidently, it only needs a few subdivision steps in order to bring these Gaussian quadrature knots into a quartic box-spline patch. In this paper, On each sub-triangle, we need to estimate its minimum subdivision times in advance for the Gaussian knot with the parameter values ( 215 , 2 15 ), ( 1115 , 2 15 ), ( 215 , 11 15 ) and ( 1 3 , 1 3 ), which are 2 times, 1 times, 1 times and 1 times respectively, then we can adaptively carry out the fast evaluation according to the box-spline formulation according to [43] . First, this approach is efficient as we do not substantially refine the given coarse mesh during the calculation process, meaning the number of unknown variables to be solved remains unchanged. Second, the hierarchical subdivision bases can maintain their linear independence according to subdivision stencils, which are only related to the local topological structure around them. Most importantly, this subdivision relationship ensures that the subdivision bases are geometrically refined without altering the object's geometric shape, i.e., the hierarchical subdivision bases maintain geometric exactness. For traditional linear basis functions, their support set is within one neighboring ring. In contrast, subdivision basis functions have the support sets within two neighboring rings, while retaining three desirable properties: unit partitioning, positivity, and linear independence. This subdivision function space belongs to the global 1 continuity class, and is 2 continuous everywhere except at extraordinary control vertices. We use subdivision surfaces to represent the physical domain and simultaneously apply the same set of hierarchical subdivision bases as the discretization method for the PDE weak form. In addition to being a versatile tool for geometric representation, smooth subdivision surfaces are also well-suited for discretizing high-order PDEs.
Given our objective to develop a linear-type scheme, linearization approaches naturally emerge as the primary candidates. Recognizing that the IEQ method [68][69][70][71] offers the dual benefits of effectively linearizing the nonlinear potential while also ensuring unconditional energy stability, we adopt this approach to numerically discretize the nonlinear potential.
The IEQ scheme requires the introduction of an auxiliary variable, ( (Ė, ), (Ė, )), to reformulate the nonlinear potential into a quadratic form. The variable is defined as:
where
and 1 , 2 > 0 and > 0 are three constants to be predetermined. Importantly, since the negative terms involving 1 and 2 can be controlled by the positive quartic potential of ( ) and ( ), the function ( , ) is guaranteed to be bounded from below. Thus, a suitable constant can always be chosen to ensure that ( , ) + remains strictly positive. We further incorporate the variable into the equations for chemical potentials (2.9) and (2.11), obtaining their equivalent forms as follows:
where
The above formulation follows the standard IEQ process, which facilitates the linearization of the potential when constructing the numerical scheme by employing implicit-explicit discretization for the nonlinear terms involving and , as well as for and . The extraction of the two quadratic terms associated with 1 and 2 is intended to ensure the 1 -stability of and (see (3.42)).
However, to achieve full decoupling, an additional auxiliary variable, , is further needed to be introduced to handle the nonlinear couplings between -ē, -ē, and --. We relate the nonlocal variable ( ) to an ordinary differential equation (ODE) with a particular structure, which is expressed as:
It is straightforward to observe that the above ODE is equivalent to ( ) 1, given that = 0 and the initial condition | =0 = 1. By incorporating all the aforementioned reformulations, the system (2.8)-(2.13) can be rewritten in the following equivalent format:
)
∇ ç ē = 0, (3.18)
with the following initial conditions:
In the formulations (3.12)-(3.19), beyond incorporating the -reformulated potentials as presented in (3.7)-(3.9) and the -ODE defined in (3.11), we further multiply the factor to all nonlinear terms. This reformulation is valid given the fact that 1, ensuring that the structure of the original system remains unchanged. The design of -ODE is based on the ''zero-energycontribution'' property satisfied by coupling terms. Consequently, the final reformulated system (3.12)-(3.19) is mathematically equivalent to the original system (2.8)-(2.13). Furthermore, this transformed system (3.12)-(3.19) also preserves the energy dissipation law, which follows a similar argument to the energy stability proof of the fully discrete scheme presented in Section 3.3. For simplicity, we omit the proof of the energy law for the PDE system (3.12)- (3.19) and only present the discrete energy law for our proposed algorithm.
We now proceed to develop the fully discrete scheme for the reformulated system (3.12)-(3.19), the weak form of which is expressed as: Find ( , , , , , ē
Computer Methods in Applied Mechanics and Engineering 436 (2025) 117733
Let > 0 be the time step size, and define = for 0 d d = + ∕ ,. We denote by ℎ the approximation of the function at time = . For the discrete formulation, in the subdivision element IGA space as described in Section 3.1, we choose the test functions ℎ , ℎ , ℎ , ℎ * 1 (ā ℎ ), ℎ * 2 (ā ℎ ), Ĕ ℎ * 1 (ā ℎ ) 3 , and ℎ * 2 0 (ā ℎ ) K 1 (ā ℎ ). Given the known approximation solutions at previous time steps = -1, , i.e., ℎ , ℎ , ℎ , ℎ , ℎ , ℎ , ē ℎ , ℎ , we now seek to compute the approximation solutions at the next time step
ℎ , using the second-order backward difference formula (BDF2):
where
(3.38)
Remark 3.1. We now outline our discretization strategy of constructing the above scheme. The key aspect is that the nonlinear term involving is discretized explicitly, while the factor itself is discretized implicitly. For the equation of , the term inside the 2 -product is treated using an implicit-explicit approach: terms to the left of the comma are discretized explicitly, while those
on the right are discretized implicitly. This implicit-explicit combination method ensures both energy stability and efficiency in the numerical scheme.
Remark 3.2. For the incompressible fluid flow system, a projection-based method is employed, enabling the decoupling of the pressure calculation from the velocity field. The convergence analysis for both the velocity and pressure fields, specifically for the Navier-Stokes equations, can be found in [91]. In this approach, ē +1 ℎ represents the intermediate velocity field, while ē +1 ℎ is the final velocity field. It is also important to note that ē +1 ℎ satisfies the discrete divergence-free condition, which can be derived by taking the 2 inner product of (3.37) with ∇ ℎ , where ℎ * 2 0 (ā ℎ ) K 1 (ā ℎ ), yielding
where the result follows from applying (3.36).
In this subsection, we present the energy stability proof for the discrete scheme (3.29)- (3.37). This analysis demonstrates that the discrete scheme preserves the fundamental property of energy dissipation.
The following two identities will be frequently used throughout the analysis, and are given by:
(3.40)
.41) Theorem 3.1. The fully discrete scheme (3.29)-(3.37) is unconditionally energy stable, namely, the following discrete energy dissipation law holds:
where, for an integer e 0, the discrete energy ̃ is denoted as
(3.43)
Proof. By taking ℎ = 2 +1 ℎ in (3.29), we get
(3.45) By taking ℎ = 2 +1 ℎ in (3.31), we get
(3.48)
Computer Methods in Applied Mechanics and Engineering 436 (2025) 117733
By multiplying (3.33) with 2 +1 and applying (3.41), we get
(3.49)
By combining (3.44)-(3.49), we get
(3.51)
We rewrite (3.37) as
By taking the 2 inner product of the above equality with ē +1 ℎ and using (3.39), we derive
Meanwhile, from (3.39) and (3.52), we deduce
Utilizing (3.53) and (3.54), we derive
(3.55)
We rewrite (3.36) as
Then, by taking the square of (3.56) and multiplying the result by 2 2 3 , we deduce
We rewrite (3.36) again as
By taking the 2 inner product of (3.58) with 3 ē +1 ℎ and applying (3.39) and (3.40), we obtain
.59) Computer Methods in Applied Mechanics and Engineering 436 (2025) 117733 Q. Pan et al. By combining (3.51), (3.55), (3.57) and (3.59), we get
(3.60)
By combining (3.50) and (3.60), and using (3.41), we get
where
We finally obtain the result (3.42) after dropping the positive term in (3.61). ¦
At first glance, the scheme (3.29)-(3.35) appears to be coupled; however, this coupling exhibits a unique structure. In the original PDE system, there are direct couplings between variables, such as -and --ē. In our proposed scheme, these direct couplings between variables are eliminated and replaced by indirect couplings between each variable and the auxiliary variable . Consequently, if we can design a method to decouple from the other variables, a fully-decoupled scheme can be constructed. Through leveraging the nonlocal properties of , we establish such a decoupled computational framework, as demonstrated below.
Step 1: We utilize +1 to decompose +1 ℎ , +1 ℎ , +1 ℎ , +1 ℎ , and +1 ℎ into linear combination forms, which are expressed as:
By applying the linear form provided in (3.62), we can decompose the scheme (3.29)-(3.33) into the following six independent systems:
Computer Methods in Applied Mechanics and Engineering 436 (2025) 117733
Each system is independent of the others and is linear. More importantly, each system in (3.63)-(3.66) consists of constantcoefficient elliptic equations, which facilitates the ease of solving the systems and ensures computational efficiency.
Step 2: we use +1 to split ē +1 ℎ into the simple combination form as
Applying this split form, according to +1 , we decompose the scheme (3.34) as the following two independent linear equations:
)
)
Step 3:
ℎ and ē +1 ℎ with = 1, 2 into (3.35), we obtain the following linear equation for +1 , that reads as
where 1 and 2 are given as
We need to demonstrate the solvability of the linear Eq. (3.72) by ensuring that the condition 3 2 -2 0 holds. First, by setting ℎ = +1 2ℎ , ℎ = 3
2 +1 2ℎ in (3.64), and ℎ = +1 2ℎ , ℎ = 3 2 +1
2ℎ in (3.66), and combine the two obtained equalities, we get
). The combination of (3.74)-(3.76) ensures that -2 e 0, which implies that the condition 3 2 -2 0 is always satisfied. Consequently, the linear Eq. (3.72) is guaranteed to be solvable.
Step 4: We update ē +1 To summarize, the four steps outlined above leverage the nonlocal property of and ensure that, at each time step, we achieve fully-decoupled computations. Each resulting equation is either a linear elliptic equation with constant coefficients or a linear algebraic equation. This decoupling structure highlights the practical efficiency of our proposed scheme, making it particularly well-suited for computational implementation.
In this section, we present several numerical examples to illustrate the effectiveness of our proposed numerical scheme (3.29)- (3.37). We begin by demonstrating the convergence rate and energy stability of the scheme, followed by a benchmark simulation of the Saffman-Taylor instability on closed surfaces, driven by rotational forces.
In this example, the computational domain is defined as a spherical surface. Specifically, we use a sphere with a radius of 1 as the test model, which is expressed as:
The initial conditions are specified as follows: ⎧ ⎪ ⎨ ⎪ ⎩ 0 = cos( ) cos( ) cos( ), 0 = 0.2 + 0.01 sin(2 ) sin( 2) sin( ), ē 0 = (cos( ) sin( ) cos( ), sin( ) cos( ) cos( ), -cos( ) sin( ) sin( )) , 0 = sin( ) sin( ) sin( ).
(4.2) The model parameters are set as: { 1 = 2 = 1.6e-3, 1 = 2 = 5e-1, = 1, = 0.4, = 3.2e-2, = 1e5, = 1, 1 = 2 = 8e-3, = 1e-5, 1 = 2 = 1, = 1e-2, = 3e3, 1 = 2 = 0.4. (4.3)
We generate four refined meshes using Loop subdivision, with the total numbers of patches/vertices for these surface models being 512∕258, 2048∕1026, 8192∕4098, and 32768∕16386, respectively. These meshes are denoted as ā ℎ , where ℎ represents the surface mesh size at the th refinement level for = 0, 1, 2, 3, and the corresponding mesh sizes are ℎ 0 = 0.2543, ℎ 1 = 0.1189, ℎ 2 = 0.0622, and ℎ 3 = 0.0326. The spatial approximation accuracy of the subdivision element IGA method was proved to be secondorder in our previous work [47], thus we set the time step sizes as = ℎ , where * (0, 1). The system has no exact solutions, so that we can set the achieved solutions on a significantly finer mesh as the reference solutions for evaluating numerical errors.
We test the convergence rate using progressively finer Loop subdivision meshes. In Fig. 4.1, we plot the 2 -norm error of all variables , , ē, and at = 2.56e-1. It is evident that the variables , and ē exhibit second-order convergence, while the variable demonstrates first-order convergence. These results are consistent with the theoretically established convergence rates.
In this example, we evaluate the energy stability of our proposed scheme by performing a spinodal decomposition simulation. Spinodal decomposition is a process in which a homogeneous mixture becomes unstable and spontaneously separates into distinct phases. The initial state is set as a uniformly perturbed binary mixture, and the system evolves from a near-uniform state to a stable configuration with two distinct phases.
In this simulation, the computational domain remains the spherical surface ā 1 , as defined in (4.1). The surface model has 131072 patches and 65538 control vertices, with vertex valence ranging from 4 to 6. The initial conditions are set as 0 = ̄ 0 + 0.01r and( , , ), 0 = 0.2, ē 0 = ÿ, 0 = 0, (4.5)
We examine two different initial values for ̄ 0 : ̄ 0 = 0.0 and ̄ 0 = 0.3. The dynamic evolution of the phase-field variable and the surfactant concentration at several time points for these two cases are shown in Fig. 4.2 and Fig. 4.3, respectively. In the first case, with ̄ 0 = 0.0, the system evolves to a final steady state that is characterized by a banded structure, where alternating Fig. 4.5. Saffman-Taylor fingering pattern with the absence of surfactant effect ( 1 = 2 = 0). Snapshots of the profiles of are taken at = 1, 3, 7, 10.
regions of the two phases form elongated bands across the domain. By contrast, for ̄ 0 = 0.3, the system develops a steady state with circular phase domains, suggesting a more localized phase separation where distinct droplets of one phase are surrounded by the other. In both simulations, the surfactant concentration profile reveals a common trend: the surfactant accumulates at the interface between the two-phase fluids, leading to a higher concentration at the interface compared to the bulk region. Moreover, using three different time steps, = 0.01, 0.005, and 0.0025, for these two simulation examples, we also present the temporal evolution trend of the total free energy, as shown in Fig. 4.4 . The curves depicting the energy evolution consistently show a monotonic decrease over time, confirming the unconditional stability of our proposed algorithm.
In this subsection, we present the Saffman-Taylor instability simulation, a classical fluid dynamics benchmark problem that illustrates the evolution and form of fine finger-like pattern structures at the interface of two fluids. When a less viscous fluid pushes a more viscous fluid between narrowly spaced plates in a Hele-Shaw cell surface, the Saffman-Taylor instability occurs, then leading to the development of an unstable interface characterized by intricate finger patterns. This phenomenon is significant in understanding interfacial instabilities in multiphase flow systems and has been widely applied in fields such as porous media flow and enhanced oil recovery. In our simulation, we implement this instability on a surface, driven by rotational forces [92,93], and examine how surfactants affect surface tension by observing the number of fingers formed. This approach helps to validate the surfactant's role in weakening interfacial tension, as observed through the dynamics of the evolving interface.
To achieve the rotational state in the Darcy equations, we need introduce an external rotational force. Specifically, we modify the right-hand side of (2.12), replacing it with the following form (as seen in [92,93]):
where the right-hand side represents the rotational force. Here, ð and are positive constants to be pre-specified, Đ is a given vector to be given, and Č denotes the outward normal vector on the surface. This rotational force induces rotational flow necessary to simulate the instability and the resulting fingering patterns. The Hele-Shaw cell surface domain is a spherical surface model with a radius of 2, expressed as: (4.9)
We first simulate the dynamics of the fingers formed in the absence of surfactant effects and then compare them to simulations with surfactant effects applied. In Fig. 4.5, with 1 = 2 = 0 (no surfactant effect), we plot the profiles of the phase-field variable at several time points. The results show that the initially smooth interface develops small protrusions, which gradually grow and extend into the characteristic viscous fingering pattern. This behavior illustrates the natural progression of instability without the stabilizing influence of surfactants.
Next, we introduce the effect of surfactants by setting the parameters 1 = 2 = 0.0004. In Fig. 4.6, the profiles of the phase-field variable are plotted at various time points. The results show that, with the presence of surfactants, the number of fingers formed significantly increases compared to the case without surfactants, highlighting the surfactants' role in reducing surface tension and enhancing finger formation. In Fig. 4.6, the profiles of the surfactant concentration variable are shown at the same time points as . The results reveal that surfactants accumulate at the interface, which promotes interface destabilization and the formation of finer, more intricate fingering patterns by reducing interfacial tension effectively. This underscores the critical role of surfactants in modulating interfacial dynamics and facilitating more complex morphologies.
To further investigate the influence of surfactants, we continue to perform the simulations by varying the coupling parameters 1 and 2 in Fig. 4.7, where the profiles of the phase-field variable are compared at the same time, = 10. We observe that, as 1 and 2 increase, the system exhibits increasingly refined and complex finger patterns, suggesting that stronger surfactant effects intensify the instability and accelerate the growth of the fingering structures. This highlights the critical role of surfactants in modulating interfacial dynamics and determining the morphology of the evolving fingers.
Finally, we perform the rotational Hele-Shaw cell surface simulation on a more complicated computational surface-a human head model, which is represented as: