In this study, we investigate numerical approximations of the hydrodynamically coupled Cahn-Hilliard phase-field model [1], pertaining to the incompressible two-phase fluid flow system, on complex curved surfaces. The Cahn-Hilliard equation, recognized as the fundamental model within the phase-field framework, is generally derived through an energy variational approach applied to the total free energy, encompassing both linear and nonlinear potentials. We note that extensive research has been dedicated to the theoretical analysis and computational techniques associated with phase-field equations, along with their adaptable applications in various systems, such as coupling with hydrodynamics to form multiphase flow systems, coupling with electromagnetic equations to form multiphase magnetic fluids or ferromagnetic fluids, and coupling with heat equations to form dendritic crystal growth systems, etc. However, it is also evident that the advancement of numerical methods and simulations for phase-field-related systems has predominantly focused on traditional rectangular or box-shaped regions, with comparatively limited exploration in the context of regions featuring complex closed manifolds or curved surfaces.
It is worth emphasizing, however, that numerical simulations of hydrodynamic systems within surface geometries are important for understanding the mechanical properties of key physical systems, including, but not limited to, soap bubble membranes [2][3][4][5], lipid bilayers [6][7][8][9][10], and interface-embedded colloidal systems [11][12][13][14][15], etc. Numerical simulations of hydrodynamical systems on curved surfaces have been relatively scarce thus far. The majority of efforts have focused on the development of spatial operator discretization methods, rather than the pursuit of achieving a numerical algorithm characterized by high spatio-temporal accuracy and unconditional energy stability. As a result, in the case of multiphase fluid flow systems coupled by the Cahn-Hilliard and hydrodynamic equations on curved surfaces -systems that are notably more complex and nonlinear -there is a significant deficiency in research efforts dedicated to developing numerical schemes and conducting computational simulations. It is important to note that the model studied in this paper serves as the foundation and source for numerous phase-field type models. Therefore, investigating the performance of this model on curved surfaces holds substantial guiding significance. In fact, we have identified only one recently introduced numerical scheme [16], specifically designed for a hydrodynamically coupled Cahn-Hilliard system on a curved surface. However, it is imperative to underscore that the numerical scheme detailed in that study suffers from a lack of energy stability and its inherent nonlinear nature, potentially incurring substantial computational expenses.
Therefore, the goal of this article is to propose an efficient numerical algorithm to solve the hydrodynamically coupled Cahn-Hilliard phase-field model, and then carry out the numerical simulation of drop dynamics under the gravity force, on the complex curved surfaces. For this purpose, we consider the discretization in the spatial and temporal directions separately. Recognizing the non-trivial nature of tackling computational challenges associated with surfaces featuring complex topologies, our objective is to devise an efficient and accurate approximation method for discretizing these surfaces, thereby effectively reducing geometric errors inherent in parametric approximations of such complexity.
The methodology of isogeometric analysis (IGA) has been witnessed to integrate the geometric simulation and the finite element solution [17,18]. It employs the same spline-type basis functions, such as non-uniform rational B-splines (NURBS) or B-splines, etc., for both the geometric representation of computer-aided design (CAD) models and the solution functions of the finite element method (FEM). This uniform representation brings us some superiority over the classical finite element basis functions. The geometric exactness can be maintained because it removes the geometric errors caused by mesh generation. It possesses high-order smoothness for basis functions, providing higher accuracy per degree of freedom through its basic knot inserting and/or order elevating method. This allows us to enhance the accuracy, efficiency, and quality of the solution method in computer-aided engineering (CAE). Subdivision is a mature technique for generating smooth surfaces/volumes in computer graphics. It can support more flexible geometry with any topology. It is well-known that subdivision surfaces can generate smooth surfaces with arbitrary topological control meshes using the quartering refinement technique [19][20][21]. Additionally, subdivision surfaces possess B-spline refinement capabilities, making them a superior choice for finite element analysis with higher-order properties in computer-aided engineering (CAE) [22][23][24][25][26].
When considering the time discretization for the system, one of the widely recognized numerical difficulties lies in the inherent stiffness of the Cahn-Hilliard equation, which is introduced by the nonlinear double-well potential incorporating the interfacial width parameter. The choice between fully implicit and explicit discretizations constitutes a well-acknowledged dilemma. A simple fully implicit scheme entails significant computational costs and potential instability, often necessitating certain time-step constraints [27]. On the other hand, a fully explicit scheme offers convenience in implementation but tends to result in computational instability [28]. Hence, a common strategy for discretizing nonlinear terms is to combine explicit and implicit methods, effectively mitigating or removing the time step's stability constraint. Prominent approaches in this context encompass nonlinear convex splitting methods [29][30][31], linear Invariant Energy Quadratization (IEQ) methods [32][33][34][35] or its variant version of Scalar Auxiliary Variable (SAV) method [36], etc. Frequently, when developing linear numerical discrete solution methods for the Cahn-Hilliard equation with a focus on achieving unconditional energy stability, the choice of the latter two methods often provides a convenient linear numerical approach. Importantly, they also ensure at least second-order time accuracy. Furthermore, as the Cahn-Hilliard equations are extended to integrate with the hydrodynamic system, giving rise to a coupled multiphase fluid flow system, more intricate numerical challenges come into play. These challenges include addressing the discretization of the nonlinear coupling terms arising from the integration of velocity and phase field variables to ensure the preservation of linearity, stability, and secondorder accuracy. To tackle this numerical challenge, we can employ the so-called ZEC decoupling method which is quite effective in handling the nonlinear couplings, see [37][38][39][40]. This method entails introducing a nonlocal auxiliary variable and devising a new ordinary differential equation (ODE) for it. The ODE is used to transform the original system into an equivalent one that is more amenable to discretization, making it ''algorithm-friendly''. In this way, the coupling terms can be simply discretized using the explicit type methods to achieve the desired numerical scheme.
Therefore, through the integration of the explicit-IEQ (EIEQ) method [32,33,[41][42][43] for linearizing the nonlinear potential, the ZEC method [37][38][39][40] for managing coupled nonlinear terms to achieve decoupled implementations, the projection method [44] for separating the computation of the velocity field from the pressure, and the IGA method for the spatial discretization, we establish a fully-discrete numerical scheme tailored for solving the hydrodynamically coupled Cahn-Hilliard system on complex curved surfaces. The numerical scheme is highly efficient, as it only necessitates solving a small number of decoupled elliptic equations with constant coefficients at each time step. We also provide a rigorous demonstration of the scheme's unconditional energy stability and solvability, and carry out several compelling numerical examples of drop dynamics on complex curved surfaces to verify the robustness and accuracy of the scheme. The developed algorithm additionally serves as a framework for a geometrically accurate decoupling scheme to solve a wide range of flow-coupled phase-field systems on complicated surfaces, while maintaining second-order time accuracy, unconditional energy stability, decoupling structure, and linearity.
The rest of this paper is organized as follows: In Section 2, we introduce the phase-field Cahn-Hilliard equation coupled with the Navier-Stokes (CH-NS) system on complex curved surfaces, and achieve the law of energy dissipation. In Section 3, we describe the fully discrete time-space scheme, which combines the subdivision-based IGA method with temporal marching techniques including EIEQ, ZEC, and projection methods. We also provide rigorous proofs of the unconditional stability of energy and the solvability of the final linear system achieved through the utilization of the nonlocal variable splitting technique. In Section 4, we present various engaging numerical examples. We give the conclusion in Section 5.
In this section, we provide a brief introduction to the CH-NS model on a complex curved surface of a closed manifold. This highly nonlinear system comprises the Cahn-Hilliard equation and the Navier-Stokes equations, commonly employed to represent the dynamics of a two-phase incompressible fluid flow system. To lay the groundwork, we initially introduce some fundamental concepts related to surfaces in the next subsection.
Some basic notations for surfaces are provided here. Assuming the domain of interest is a surface, we proceed with surface parameterization. For a sufficiently smooth and orientable surface ∶= {𝐱(𝑢 1 , 𝑢 2 ) ∈ R 3 ∶ (𝑢 1 , 𝑢 2 ) ∈ 𝒟 ⊂ R 2 }, with the help of its parameterization (𝑢 1 , 𝑢 2 ), we denote the following derivative formula
(2.1)
We have the coefficients of the first fundamental form on a surface, 𝑔 𝛼𝛽 = (𝐱 𝑢 𝛼 , 𝐱 𝑢 𝛽 ), 𝛼, 𝛽 = 1, 2, and the important surface normal
where (⋅, ⋅), ⋅ × ⋅ and ‖ ⋅ ‖ stand for the usual inner product, cross product and Euclidean norm in R 3 respectively. We need denote [𝑔 𝛼𝛽 ] = [𝑔 𝛼𝛽 ] -1 and 𝑔 = det[𝑔 𝛼𝛽 ] before introducing the following classical differential geometric operators.
On the function space 𝐶 1 (), for any 𝑓 ∈ 𝐶 1 (), the tangential gradient operator ∇ 𝑠 is defined as
For a vector-valued function 𝐟 = [𝑓 1 , … , 𝑓 𝑘 ] 𝑇 ∈ 𝐶 1 () 𝑘 , the gradient ∇ 𝑠 is defined as
On the smooth vector field 𝐯 ∈ [𝐶 1 ()] 3 , the divergence operator div 𝑠 is defined as
(2.4)
For any 𝑓 ∈ 𝐶 2 (), the Laplace-Beltrami operator (LBO) 𝛥 𝑠 is defined as
(2.5)
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
Suppose be a sufficiently smooth and orientable surface. We use 𝜙(𝐱, 𝑡) as a phase-field variable to mark the volume fraction of the two different fluid components in the fluid mixture, i.e.,
These two distinct regions of 𝜙(𝐱, 𝑡) are smoothly attached by a thin spanned region with a width of 𝑂(𝜀), 𝜀 ≤ 1. To simplify the system, we make the assumption that the two fluids share the same density. Introducing the velocity field 𝐮 in the fluid mixture, we assume the total free energy to be
which is a summation of the gradient potential |∇ 𝑠 𝜙| 2 and Ginzburg-Landau type double-well functional 𝐹 (𝜙) = 1 4𝜀 2 (𝜙 2 -1) 2 . The model parameter 𝜆 depends on the surface tension, which is defined as
𝜆 𝜀 [45]. The gradient entropy represents the hydrophilic type of interaction or the tendency to mix, and the other double-well potential describes the hydrophobic type of interaction or the tendency to separation.
The CH-NS system is characterized by adopting the 𝐻 -1 gradient flow approach for the phase-field variable. Simultaneously, the system assumes that the fluid is incompressible, adhering to the generalized Fick's law, which dictates that the mass flux is proportional to the chemical potential gradient. This combination of conditions results in the following formulation for the CH-NS system:
)
)
)
where 𝑀 is the mobility parameters, 𝜔 is the chemical potential, 𝜙∇ 𝑠 𝜔 is the surface tension term, ∇ 𝑠 (𝐮𝜙) is the advection, 𝜈 is the fluid viscosity, 𝑝 is the pressure, and 𝑓 (𝜙) = 𝐹 ′ (𝜙) = 1 𝜀 (𝜙 3 -𝜙) and the variational derivative 𝜔 = 𝛿𝐸 𝛿𝜙 . In this article we only treat closed surface domains, and remove all boundary conditions.
Next, we will show the character of the energy dissipation for the CH-NS dynamical system (2.9)-(2.12), as follows.
Lemma 2.1. The CH-NS dynamical system (2.9)-(2.12) satisfies
where 𝐸(𝐮, 𝜙) is defined as (2.8).
Proof. By taking the 𝐿 2 inner products of (2.9) with 𝐮, considering (2.10) and using the integration by parts, we get
where we apply (𝐮 𝑡 , 𝐮) = 1 2 𝑑 𝑑𝑡 ‖𝐮‖ 𝟐 . By taking the 𝐿 2 inner products of (2.11) with -𝜔 and integration by parts, we get
(2.15)
By taking the 𝐿 2 inner product of (2.12) with 𝜙 𝑡 and integration by parts, we obtain
By summing up the aforementioned three equalities and consolidating all three terms associated with 𝑑 𝑑𝑡 on the left-hand side, we derive (2.13). □
In the preceding proof, the cancellation of inner products involving nonlinear terms is attributed to the properties stated as
(2.17)
These crucial zero equalities, referred to as the 'zero-energy-contribution' (ZEC) property (cf. [46,47]), arise from integration by parts, causing these nonlinear terms to vanish and not contribute to the total energy. This property serves as a key element guiding our approach to decoupling strategies.
and 𝑛 is its valence. (c): Edge mask.
In this section, we present the spatiotemporally fully discrete scheme for the model system (2.9)-(2.12). The IGA paradigm unifies the shape parameterization methods adopted in CAD and the discretization techniques used in FEM. We consider the domains of geometrically and topologically complex surface. Subdivision surfaces originate from CAD (see [19][20][21]) and are widely applied in computer animation applications. The subdivision technique can generate smooth surfaces on unstructured meshes, which is similar to B-splines or NURBS on structured meshes. In the FEM method, subdivision surfaces have superior conditions for square integrability with second-order derivatives. We adopt the recently developed subdivision-based IGA framework for spatial discretization and introduce the EIEQ method, merging the ZEC concept and the projection approach for time marching [41][42][43]. We prove that the discrete scheme is unconditionally energy stable and present the corresponding decoupling implementation approach to obtain high computational efficiency.
The well-known IGA methodology adopts the NURBS/T-splines [48][49][50][51], which uses the same set of basis functions to describe the solution domain and implement the numerical computation of PDEs. The framework of IGA not only develops the seamless integration between FEM and CAD, but also has higher numerical exactness than FEM through 𝑝-refinement, ℎ-refinement, and even 𝑘-refinement with the knot insertion and/or order elevation. Therefore we can improve the accuracy of the numerical simulation without destroying the original geometric properties, which removes the interactive communication with the CAD system.
Subdivision technology was proposed in computer graphics area, which can construct smooth surfaces from arbitrary topological meshes by designing a set of simple and efficient refinement schemes [19,20], and handle complicated geometric models while maintaining original characteristics near boundaries through straightforward extensions, such as concave/convex angles and sharp/smooth creases. Subdivision is compatible with NURBS, which has the capability of the refineability of B-spline techniques. Both Loop subdivision [24,26] and Catmull-Clark subdivision [25] have been utilized in IGA. The subdivision-based IGA methodology can be viewed as the natural choice for higher-order FEM in engineering practice, see [22][23][24]26] as well.
In this work, we apply loop subdivision surfaces based on triangle meshes, or more exactly, the essential subdivision shape functions, for the discretization method of the weak form. The surface subdivision technique possesses superior conditions of highorder smoothness and suitability for any topological structure, which is widely used in high-order model such as shell model [22,23]. We can generate spline-class surfaces using repeated refinement on a given initial mesh, which is essentially a discretization method. In addition to the flexible expression of geometric shapes, the high-order smoothness of subdivision surfaces also makes them ideal for solving high-order PDEs.
We adopt Loop subdivision approach [20], which creates quartic box-splines defined on three-directional triangular meshes with arbitrary topological structure. We know that the quartic box-splines are defined on a translationally invariant three-dimensional mesh, which means that the mesh is only composed of triangular patches, and each control vertex is related to adjacent triangular patches. As different from quartic box-splines, for Loop scheme, each control vertex can have any number of triangular patches, and the resulting refiner surfaces can be at least 𝐶 1 continuous. In each refinement step, the refined meshes are achieved by repeating quadrilateral segmentation of the previous mesh, with each refinement step subdividing the triangle into four smaller triangles. The position of each control vertex is recalculated as a weighted average of itself and its neighboring control vertices. The weight coefficients are represented in the form of subdivisions stencils. Vertex stencils are introduced to recalculate the coordinates of existing control vertices of the coarse mesh (see Fig. 3.1(𝑎) and (𝑏)). Edge stencils are introduced to compute the coordinates of control vertices generated on the edge during each quartering (see Fig. 3.1(𝑐)). It should be noted that the weight coefficients only depend on the connect relationship of the mesh and not on the actual positions of the control vertices. We need to compute quantities, such as the tangential vector and curvature, at any location on the parameterized subdivision surface. Loop subdivision can generate box-splines on patches supported only by regular control vertices, i.e.,
We use (𝑢 1 , 𝑢 2 , 1 -𝑢 1 -𝑢 2 ) denote the barycentric coordinates of the unit parametric triangle and 𝐱 𝑘 𝑖 denote its 2-ring surrounding neighbor control vertices of the control triangle (see Fig. 3.2(𝑎)), and use 𝐵 𝑖 as the quartic box-splines (see the details in [21]). A patch containing irregular control vertices whose valence is different from 6 can also be parameterized with box spline shape functions. It is crucial to recognize that in the process of refinement, only regular vertices are generated through the original coarse triangular mesh divided by quartered separating, which means that the valence of all new vertices is 6 (see Fig. 3.2(𝑏) and (𝑐)). Obviously, the total number of irregular control vertices in the mesh remains unchanged, and their influence domain finally closes to the limit of zero. Therefore, in the process of repeated subdivision and refinement, we can obtain more and more regular patches. We should note that Stam proposed a fast parameterization method for subdivision surfaces, which coincides with the spline surface on the regular patch [21].
We next analyze the smoothness and convergence properties of loop subdivision schemes. For the convergence properties of subdivision schemes at each control vertex, we should consider their one-ring adjacent control vertices. We denote an assembly matrix
where 𝐗 0 is a (𝑛 + 1) × 3 matrix, and valence 𝑛 denotes the number of its connected edges. The superscript 0 indicates the subdivision level. We can compute the positions of 𝐗 𝑘 , where the initial positions 𝐗 0 is multiplied 𝑘 times with the subdivision matrix 𝐒
We can refer the fast algorithm about the subdivision matrix 𝐒 𝑘 as [21] because it depends only on the valence of the control vertex and unrelated to the control vertex positions. A limit surface can be achieved as the subdivision level 𝑘 → ∞, whose every control vertex has the explicit position, as described in the following Lemma 3.1. The limit surface of Loop subdivision is 𝐶 2 everywhere except at the extraordinary vertices where it is 𝐶 1 .
Lemma 3.1. Let 𝐱 𝑘 0 be a control vertex of valence 𝑛 on the mesh 𝛺 𝑘 ℎ , and 𝐱 𝑘 𝑗 , 𝑗 = 1, … , 𝑛, be its 1-ring neighbor-controlled vertices. All these vertices converge to a single position
as the subdivision time 𝑘 → ∞ (see [20] for the proof).
A quartic box spline shape function contains two rings of neighboring elements as its support, which consists of piecewise polynomial surface patches that run smoothly across the surrounding patch boundaries. The finite element solution procedure is implemented on the limit surface. Loop basis functions naturally satisfy the 𝐻 2 -smoothness requirement of the fourth-order PDEs and exactly describe arbitrary topology, which provides a uniform mathematical description for representing both the geometry Q. Pan et al. and the physics character of the complex geometry. We perform actual calculations at the Gaussian integral points on the limit surface (3.4), so we need to establish a mapping from the unit parameter triangle to the control surface patch, where we can refer it to [24][25][26].
Firstly, according to the initial version of IEQ strategy, we only need to define a local-type auxiliary variable 𝑈 as
Here, the negative term -𝑆 2𝜀 2 𝜙 2 can be bounded for any 𝑆, and the constant 𝐵 needs to be positive to make sure the summation of all terms in the square root is always positive. The goal of designing the variable 𝑈 is to guarantee the original energy (2.8) to be quadratic. Using the new variable 𝑈 , the new CH-NS system can be rewritten as
)
)
)
)
where the new term
Secondly, we modify the IEQ method to the explicit-IEQ (EIEQ) method combined with the ZEC method by introducing an auxiliary nonlocal variable and design a special type but trivial ODE for it. We introduce another nonlocal variable 𝑄(𝑡) and an ODE system related to it that reads as Obviously, from (2.17), we get 𝑄 𝑡 = 0. Using the initial condition 𝑄| 𝑡=0 = 1, we get the exact solution of (3.12) is 𝑄(𝑡) ≡ 1. The last two terms in (3.12) correspond to the EIEQ technique [32,33,[41][42][43], while the first three terms in (3.12) belong to the ZEC decoupling method [37][38][39][40]. Note that the objective of designing 𝑄 in this manner is to achieve the decoupled type scheme while ensuring the energy stability, as detailed in Theorem 3.2 and Section 3.2.4. Therefore, by combining the two new variables 𝑈 and 𝑄, we rewrite the system (3.6)-(3.10) into the following equivalent form:
)
)
)
)
)
The initial conditions read as
Note that we multiply the surface tension term 𝜙∇ 𝑠 𝜔 and the advection term ∇ 𝑠 ⋅ (𝐮𝜙) by the nonlocal variable 𝑄. However, it is important to emphasize that (3.13)-(3.18) remains unaffected as 𝑄(𝑡) ≡ 1. In addition, integrating (3.17) with respect to time 𝑡 and incorporating the initial condition (3.19) enables us to readily obtain the original system (2.12).
The inherent characteristic of the hydrodynamically coupled phase-field system is its adherence to the energy dissipation law. Consequently, the newly transformed system (3.13)-(3.18) also preserves the structure of the energy dissipation law, as stated in the following Theorem 3.1.
In this subsection, we describe the fully discrete scheme for solving the transformed dynamical system (3.13)- (3.18). In what follows, the time step size is denoted as 𝛿𝑡 > 0 and 𝑡 𝑛 = 𝑛𝛿𝑡. The test functions are chosen as 𝜉 ℎ , 𝜃 ℎ , 𝜗 ℎ , 𝜍 ℎ , 𝜁 ℎ ∈ 𝐻 2 ( ℎ ), which is the finite element space supported by Loop subdivision elements described in Section 3.1.
Assuming that ũ𝑘 ℎ , 𝜙 𝑘 ℎ , 𝜔 𝑘 ℎ , 𝑈 𝑘 ℎ , 𝑄 𝑘 , 𝐮 𝑘 ℎ for 𝑘 = 𝑛, 𝑛 -1, 𝑛 -2, and 𝑛 ≥ 2 are known, we proceed with the update:
(3.27) Furthermore, we calculate ũ𝑛+1 ℎ , 𝜙 𝑛+1 ℎ , 𝜔 𝑛+1 ℎ , 𝑈 𝑛+1 ℎ , 𝑄 𝑛+1 , 𝐮 𝑛+1 ℎ by the following second-order difference formula (BDF2):
and
Q. Pan et al.
Remark 3.1. The time derivative is discretized using the BDF2 method. However, achieving a high-order formula like BDF3/4 poses significant challenges, primarily due to stability issues. Additionally, a direct discretization method is used in the above scheme. Specifically, linear components are implicitly treated, while nonlinear components involving a nonlocal variable 𝑄 are treated through the explicit-implicit schemes. It means that we implicitly process the term 𝑄, and explicitly process all other terms. The decoupling method of the above scheme is obviously coupled and must fully utilize the nonlocal characteristics of the 𝑄 variable.
Next, we will demonstrate the energy stability of the developed scheme, followed by the implementation details of decoupling.
Remark 3.2. We adopt a projection-based approach for the incompressible fluid dynamic system. This algorithm can decouple the calculation of pressure and velocity fields. The results show that the velocity field and pressure field can obtain second-order and first-order accuracy in time, respectively, and strict error estimates are provided [52]. The reason for the first-order accuracy loss of pressure is the application of artificial homogeneous Neumann boundary conditions to the pressure [53].
The final velocity field is denoted as 𝐮 𝑛 ℎ , and the intermediate velocity field is denoted as ũ𝑛 ℎ , which satisfy the divergence-free condition in a discrete sense. For ∇𝑣 ℎ , 𝑣 ℎ ∈ 𝑂 ℎ , we take the 𝐿 2 inner product of (3.35) to obtain
Therefore, from (3.34), we derive
(3.37)
We aim to demonstrate that the proposed discrete scheme (3.28)- (3.35) unconditionally satisfies the energy dissipation law, akin to the continuous system (3.13)- (3.18). Prior to presenting the proof in Theorem 3.2, we introduce two identities that will be repeatedly used:
where, for an integer 𝑘 ≥ 0, the discrete energy 𝐸 𝑘 is defined as
(3.41)
Proof. By taking the test function 𝝃 ℎ = 2𝛿𝑡ũ 𝑛+1 ℎ in (3.28), we get
We rewrite (3.35) as
Taking the 𝐿 2 inner product of the above equality with 𝐮 𝑘 ℎ with 𝑘 = 𝑛 + 1, 𝑛, 𝑛 -1, we derive
and
where (3.37) is used. We reformulate the first term in (3.42) as
where we add the underlined term -3𝐮 𝑛+1 ℎ + 3𝐮 𝑛+1 ℎ , which is zero.
For the first split term (3ũ 𝑛+1 ℎ -3𝐮 𝑛+1 ℎ , ũ𝑛+1 ℎ ) in (3.46), we deduce
where the added underlined term is actually zero, as can be seen from (3.44).
For the second split term (3𝐮 𝑛+1 ℎ -4𝐮 𝑛 ℎ + 𝐮 𝑛-1 ℎ , ũ𝑛+1 ℎ ) in (3.46), we deduce
where the underlined term is zero due to (3.45), and the last equality is derived using (3.39). Thus, the first term in (3.42) can be written as
(3.49)
We rewrite (3.35) as
Taking the 𝐿 2 inner product of the above equation with itself and multiply the result with 3 2 , we derive
From (3.37), we get 2𝛿𝑡(𝐮 𝑛+1 ℎ , ∇ 𝑠 𝑝 𝑛+1 ℎ ) = 0, thus (3.51) can be rewritten as
We rewrite (3.35) as
By taking the 𝐿 2 inner product of the above equation with 3𝐮 𝑛+1 ℎ , and using (3.37), we obtain
Then we get by (3.54)
(3.55)
By combining (3.42), (3.46), (3.52), and (3.55), we get
By taking 𝜃 ℎ = 2𝛿𝑡𝜔 𝑛+1 ℎ in (3.29), we get
By taking 𝜍 ℎ = 4𝜆𝛿𝑡𝑈 𝑛+1 ℎ in (3.31), we get
where (3.39) is used. By multiplying (3.32) with 2𝛿𝑡𝑄 𝑛+1 and using (3.39), we obtain
We combine (3.57)-(3.60) to get the following
For the first two terms of (3.61), by using (3.39), we have
(3.62)
By combining (3.56), (3.61) and (3.62), we obtain
In this section, we introduce our decoupling computational steps by using a splitting technique for treating the systems (3.28)- (3.32).
Step 1: we use 𝑄 𝑛+1 to split ũ𝑛+1 ℎ into the linear combination form that reads as
By applying the linear form given in (3.64), we can split the scheme (3.28) as follows:
).
In this section, we demonstrate the effectiveness of the numerical scheme (3.28)-(3.35) developed through a series of numerical simulations. We validate its energy stability and convergence properties. The simulations cover scenarios including spinodal decomposition and droplet dripping dynamics on a complex curved surface, driven by gravitational forces. We perform the spatial discretization by use of the surface subdivision method in Section 3.1. A robust GMRES iterative solver is applied to solve the final systems derived from the variational format because these systems are highly sparse, where a small enough tolerance is set to obtain the proper convergence.
The considered surface is defined as a sphere with a radius of 1, as follows: We generate four refined meshes by Loop subdivision, where the total numbers of patches/vertices for these surface models are 512∕258, 2048∕1026, 8192∕4098, 32768∕16386 respectively. We denote them as ℎ 𝑖 , where ℎ 𝑖 represents the surface mesh size at the 𝑖th refinement level for 𝑖 = 0, 1, 2, 3, and ℎ 0 = 0.2543, ℎ 1 = 0.1189, ℎ 2 = 0.0622, ℎ 3 = 0.0326. The subdivision-based IGA method has been shown to be second-order accurate in our previous work [24], hence we set the time step sizes as 𝛿𝑡 𝑖 = 𝐶ℎ 𝑖 with 𝐶 ∈ (0, 1).
Since the exact solutions of this equation system are unavailable, we may set the solutions using a significantly small mesh size and consider them as reference solutions for computing numerical errors. For all unknown variables 𝜙, 𝐮 and 𝑝, we plot their 𝐿 2 numerical errors at the time 𝑡 = 2.56e-1 in Fig. 4.1, where the two variables 𝜙 and 𝐮 have the second-order convergence rate and the variable 𝑝 has the first-order convergence rate. These results align well with the theoretically anticipated convergence rate.
This simulation focuses on the investigation of the dynamical modeling of the spinodal decomposition, which can show the effectiveness of our developed scheme. We set the initial state as a uniformly perturbed binary mixture within a certain range. The system goes through a transition from a homogeneous state to a stable situation with two distinct phases, driven by the essential growth of concentration fluctuations.
The interest domain is a closed sphere that reads as
The total numbers of patches/vertices for the surface model is 131 072/65 538, and the span of the vertex valence is 4 to 6. The initial conditions are given by 𝜙 0 = φ0 + 0.01rand(𝑥, 𝑦, 𝑧), 𝐮 0 = 𝟎, 𝑝 0 = 0,
where the term rand(𝑥, 𝑦, 𝑧) is the random number in [-1, 1] with zero mean. The model parameters are set as
In Fig. 4.2, two simulation results corresponding to initial values of φ0 = 0.0 and φ0 = 0.3 are presented, respectively. We plot snapshots of 𝜙 on the spherical surface over a period of time. Two final stable states correspond to a banded and a circular equilibrium state, respectively. Additionally, in Fig. 4.3, we plot the total free energy functionals over time for the two performed simulations, where we choose three different time steps 𝛿𝑡 = 0.0025, 0.005, 0.01 to test the energy decay situation. All of the evolving curves of the total energy monotonically decrease over time which confirm the unconditional stability of our method.
In this simulation, we study the drop-dripping dynamics under the gravitational force on the complex surfaces. To mimic the gravitational force, we use the Boussinesq approximation by augmenting the fluid momentum equation (2.9) with an extra force that reads as (cf. [44,54])
where we assign 𝑔 as a gravity constant, 𝜌 1 and 𝜌 2 are the density of the heavier drop and lighter medium, respectively. Note that the Boussinesq approximation is effective only for small density differences. When dealing with a large ratio in density variance between two phases, a different model must be employed, and the corresponding algorithms need to be redesigned due to the presence of distinct nonlinear terms, see [41,55].
We simulate the droplet's dripping behavior under the influence of gravity on the complex surface domain of the human head. We denote it as 2 , which is embedded in 𝛺 = (-0.346, 0.346) × (-0.497, 0.497) × (-0.495, 0.495). The total numbers of patches/vertices for the surface model is 268 686/134 345, and the span of the vertex valence is 4 to 12. Two droplets, one initially circular and the other elliptical, are simulated to investigate how the dripping process evolves.
The initial conditions for these two drops are set as follows: In Fig. 4.4(a) and (b), we plot snapshots of the variable 𝜙 for two initially differently shaped droplets at different times. It can be seen that the initially circular droplet essentially descends as a cohesive whole, whereas the initially elliptical droplet forms a more elongated filament. This observation suggests that the initial shape of the droplet can indeed have a notable impact on the dynamics of its dripping process on the surface. We denote another surface as 3 , which is a splayed surface, and embedded in 𝛺 = (-0.61, 0.61) × (-1.27, 1.27) × (-0.26, 0.26). The total numbers of patches/vertices for the surface model is 98 304/49 150, and the span of the vertex valence is 4 to 8.
We set two initial conditions for this model, one of which reads as: In Figs. 4.5 and 4.6, snapshots of 𝜙 at different times are plotted. For these two scenarios, we observe distinct droplet dripping dynamics influenced by both the geometric shape of the surface and the initial position of the droplet. In the first scenario, at the initial moment, the droplet is positioned to the upper right of the surface of the first ring. In contrast, in the second case, the droplet starts right in the middle of the first ring. For the first case, the droplet consistently descends along the right side, while in the second case, it falls along both sides of the ring. During the descent in both cases, the center portion of the droplet undergoes significant elongation, forming a long filament. However, in the second case, the droplet falling along both sides is pinched off into two parts by gravity. Upon reaching the second ring, the two separated parts rejoin and accumulate at the bottom of the surface, resulting in a final fluid interface appearing perfectly horizontal in both cases.
In this work, we solve the hydrodynamically coupled Cahn-Hilliard phase-field model for a two-phase fluid flow system on complex curved surfaces. We conduct the spatial discretization using the IGA methodology, specifically employing the Loop subdivision technique. This technique possesses high smoothness and hierarchical refinement properties, enabling accurate representation of complex surfaces with arbitrary topologies. We integrate the EIEQ method, the ZEC method, and the projection method to discretize