The Cahn-Hilliard equation, as one of the classic representative systems in the diffusive phase-field approach, has been extensively used in the field of fluid dynamics to capture a variety of complex interfacial phenomena in multiphase flows [1][2][3]. Combining the Cahn-Hilliard equation with various flow regimes, e.g., coupling with the incompressible Navier-Stokes equations to form the Cahn-Hilliard-Navier-Stokes system, can be used to simulate the dynamic phenomena of a liquid droplet [4], or the ascent of a gas bubble under the force of gravity [5][6][7], etc. The so-called Cahn-Hilliard-Darcy system [8][9][10], which is a combination of the Cahn-Hilliard equations and the Darcy equations, where the coupling is typically realized by the surface tension and advection, has been extensively utilized to simulate the dynamics of a two-phase incompressible fluid system that is confined within a Hele-Shaw chamber, which refers to a setup consisting of two parallel plates with a narrow gap, filled an immiscible viscous fluid comprising two phases. While it is worth noting that the vast majority of Hele-Shaw flow studies have mostly focused on the investigation of two-dimensional flows between flat plates, in recent years, an increasing number of experimental and theoretical/numerical works have been devoted to the study of pattern-forming phenomena in parallel nonplanar surfaces, such as spherical, cylindrical, conical, toroidal, and saddle-shaped Hele-Shaw chambers [11][12][13][14][15][16][17][18][19][20][21][22], to comprehend how unstable viscous fingering (also known as Saffman-Taylor instability) responds to changes in the geometrical properties of the nonplanar Hele-Shaw chambers, or how the special geometrical features of the nonplanar Hele-Shaw chambers activate the fluid fingering instability, and so on. Moreover, the numerical study of interfacial dynamics using two-phase viscous fluid flow models on complex surfaces to simulate the Saffman-Taylor instability is not only of scientific importance, but also of wide practical engineering applications, such as the pressure molding of molten metal and polymer materials [23], and the formation of coating defects during the drying of coating films [24], etc. Hence, the goal of this article is to propose an efficient numerical algorithm to solve the Cahn-Hilliard-Darcy system, and then simulate the fingering instability in a nonplanar Hele-Shaw cell, on the complex curved surfaces.
Designing efficient, fully discrete numerical algorithms for complex, nonlinearly coupled systems, such as the Cahn-Hilliard-Darcy system, is by no means a simple task, particularly when the computational domain is a surface with complex topology and it is expected to produce a scheme with simple implementation, high order accuracy, and unconditional energy stability. The first primary numerical challenge encountered is spatial discretization, i.e., devising efficient and accurate methods for discretizing complex surfaces. In this article, we discretize the curved space by adopting the subdivision-based isogeometric analysis (IGA) method [25][26][27][28][29][30][31][32]. The concept of IGA is to incorporate geometric computer-aided design (CAD) modeling and finite element (FEM) approximating using spline-type basis, e.g. B-splines and non-uniform rational B-splines (NURBS), which have been extensively applied in CAD but rarely used in FEM analysis. The methodology of the IGA based on the B-spline basis function has the following benefits throughout the geometric design and analysis workflow: first, it eliminates errors in finite element mesh generation due to human intervention, and second, it accelerates the data exchange between design and analysis, which is essential for engineering design optimization. Without interactively communicating with the CAD system, numerical accuracy can be easily improved while maintaining geometric invariance of the computational region through its basic knot inserting and/or order elevating methods. In addition, surface subdivision is a successful technology of generating smooth surfaces with arbitrary topological control meshes using simple refinement schemes [33][34][35][36] and has B-spline refinement capabilities, so it can be a spontaneous fit for finite element applications as a method with higher-order properties in engineering practice [27,28,37,38].
Time discretization is the second numerical difficulty because the Cahn-Hilliard-Darcy system suffers from high stiffness due to its extensive coupling and nonlinear terms (including nonlinear double-well potential, surface tensions, and advection terms), and so the key issue in developing a numerical scheme for time discretization is how to discretize these nonlinearly coupled terms in order to achieve both energy stability and higher-order time accuracy. For the double-well potential, we recall that there are numerous discretization approaches, such as the implicit derivative scheme [39], the quadrature scheme [40][41][42], the convex splitting method [8,[43][44][45][46], the linear stabilization method [47,48], the Invariant Energy Quadratization (IEQ) method [49][50][51][52][53][54][55][56][57] and its various forms of the Scalar Auxiliary Variable (SAV) [58,59] method, etc. In this article, we adopt the explicit-IEQ method, which is also a variant form of the IEQ method, and it has the advantage of easy implementation, energy stability, second-order temporal accuracy, etc. For nonlinear terms such as surface tension and advection terms, it is well-known that a simple explicit or implicit discretization of them rarely results in a second-order fully-decoupled linear scheme that satisfies energy stability, but at best results in, for instance, a nonlinear scheme [60], a second-order linear but coupled scheme [8,61], or a decoupled scheme but only first-order accurate in time [62,63]. To obtain a fully decoupled second-order scheme, we use the so-called Zero-Energy-Contribution (ZEC) method [64][65][66][67], the main idea of which is to design a nonlocal type auxiliary variable and to construct a special form of ODE, followed by a nonlocal splitting scheme to obtain the desired fully-decoupled computation of energy stability and second-order accuracy in time.
By combining the aforementioned methods and techniques for spatial and temporal discretizations with the second-order projection scheme for the incompressible Navier-Stokes equations, we arrive at an effective numerical strategy to solve the Cahn-Hilliard-Darcy system on the complex curved space, which maintains unconditional energy stability and linearity. It is also very effective to solve, as only a few decoupled elliptic equations with constant coefficients need to be solved at each time step. In addition, we rigorously demonstrate the scheme's unconditional energy stability and solvability. Various attractive numerical examples, such as the Saffman-Taylor fingering instability on various complex curved surfaces caused by rotation, are simulated to verify the robustness and exactness of our scheme. In addition, the proposed fully discrete scheme provides a framework for a geometrically precise decoupling scheme to solve various flow-coupled phase-field systems on complex surfaces.
The remainder of this paper is structured as follows. In Section 2, we present the Cahn-Hilliard-Darcy system confined in the Hele-Shaw cell on complex curved surfaces and deduce its dissipation law. In Section 3, we describe the fully discrete time-space scheme that integrates the subdivision-based IGA method for spatial discretization with the temporal marching techniques EIEQ, ZEC, and projection methods. We also provide rigorous proof of the unconditional stability of energy and the solvability of each time step. In Section 4, we simulate some numerical implementations to demonstrate the scheme's accuracy and efficiency. Section 5 provides some concluding remarks.
In this section, we present the Cahn-Hilliard-Darcy phase-field model for a two-phase Hele-Shaw flow system on a closed complex surface.
Q. Pan et al.
Some basic notation for surfaces is given here. Consider the interest domain is a surface, we need to proceed the parameterization on the surface. 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
We have the coefficients of the first fundamental form on a surface, = ïĖ , Ė ð, , = 1, 2, and the important surface normal
We need the denotation [ ] = [ ] -1 and = det[ ] before introducing the 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 = div (∇ ).
(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
For any functions , * 1 (ā), two classical inner products on the surface ā are described as
On a smooth and orientable surface ā, we introduce (Ė, ) as a phase-field variable to mark the volume proportion of the two incompatible fluid parts in their mixture, i.e.,
where a slight smooth layer with the transition span ( ) is adopted to link the two distinct regions of (Ė, ).
The total free energy is postulated as
which involves Ginzburg-Landau type double-well functional ( ) = 1 4 2 ( 2 -1) 2 and the gradient potential |∇ | 2 . The gradient entropy is indicative of the hydrophilic nature of interactions, reflecting a propensity for mixing. Conversely, the nonlinear double-well potential characterizes hydrophobic interactions, signifying a tendency towards separation.
The fluid motion is hypothesized to conform to the mechanical principles observed in porous media, and adhere to the generalized Fick's law, which states that the mass flux is directly proportional to the gradient of chemical potential. The phase-field Cahn-Hilliard equation, linked with Darcy's flow (referred to as Cahn-Hilliard-Darcy model), can be expressed based on the free energy functional (2.8), as follows (see also in [8][9][10]65]),
= (-+ ( )), (2.10)
where is the mobility, ē is the seepage velocity, is the pressure, is the chemical potential, ∇ is the surface tension term, ∇ ç (ē ) is the advection, is a positive parameter, denotes the dimensionless hydraulic conductivity, and is the pressure. In the phase-field equation (2.9)-(2.10), the variational gradient flow approach in the -1 space is used. The Darcy equation (2.11)-(2.12), includes the consideration of the temporal derivative of the seepage velocity ē for incompressible flows occurring within a porous media (see also [8,68,69]).
The Cahn-Hilliard-Darcy system (2.9)-(2.12) adheres to the principle of energy dissipation, which may be obtained by the following procedure. By taking the 2 inner product of (2.9) with -, and using integration by parts, we get
(2.13)
By multiplying the 2 inner product (2.10) with , and using integration by parts, we obtain
By taking the 2 inner products of (2.11) with ē, and using the integration by parts and (2.12), we get
By combining the obtained three equations, we get the subsequent principle of energy dissipation as
where
(2.17)
Note that when deriving the energy law (2.16), the two inner product terms involving the surface tension and advection are canceled, as
This property is the so-called ''zero-energy-contribution'' property and can be used to construct the fully-decoupled scheme as [65][66][67].
In this section, we present the methodology for constructing a spatiotemporally fully-discrete scheme for the Cahn-Hilliard-Darcy system (2.9)-(2.12). To discretize the spatial domain, the subdivision-based IGA framework is employed by us (See [56,57]). For the temporal discretization, we utilize EIEQ method in combination with the ZEC method and the projection method for time marching. The rigorous proof of unconditional energy stability is presented, along with a comprehensive explanation of the methodology for implementing the decoupled type computation.
The framework of IGA merges geometric representation and finite element analysis by using spline-type basis functions, where we use the same set of shape functions for geometric CAD modeling and finite element analysis. In practical applications, traditional finite element shape functions used for finite element analysis are not universal enough for complex geometric CAD modeling. The idea of IGA provides us with a promise that a unified representation can facilitate different geometric modeling and analysis tools currently available. The subdivision surface is a successful geometric modeling method in computer graphics and computeraided design. Subdivision uses B-spline curves to generate smooth surface reprocessing through continuous refinement on control mesh patches with the tensor product method, or, in other words, the relationship between two proportions and the average. The reprocessing ability of B-spline curves enables the shape functions on the coarse mesh to be represented as a linear combination of shape functions on finer meshes. As an optimism, subdivision surfaces provide an ideal isogeometric description idea for geometric CAD modeling and numerical analysis. ) 2 ] and is its valence. (c): Mask of an edge.
The quadric box spline curves are defined on a translation-invariant 3D mesh, i.e., the control mesh is composed only of triangular patches, and each control vertex is surrounded by six adjacent triangular patches. Loop subdivision scheme (c.f. [34]) extends the quartic box splines to meshes with any topological structure. As opposed to the topological structure of quartic box-splines, In the rule of Loop subdivision, the number of control triangular patches surrounding each control vertex could be arbitrary, and the generating refiner surfaces have the property of at least 1 continuity. In each refinement step, the triangular control mesh is refined by quartering its control patch, and the positions of each control vertex are recalculated as a weighted average of the positions of its neighboring control points. The weight coefficients are represented in the elements of subdivision masks. In Fig. 3.1(a) and (b), masks for control vertices are used to recalculate the coordinates of existing control vertices on the coarse control mesh (see Fig. 3.1(a) and (b)). A mask for an edge is applied to compute the coordinates of control points created on the edge in the process of each quartering (see Fig. 3.1(c)). It should be noted that the weight coefficients are only related to the connected way of the mesh and not to the real coordinates of control points.
The position and curvature at arbitrary position and other related quantities are needed to be computed. In the procedure of Loop subdivision, box-splines are used to represent patches with twelve regular control points, which is
The parameters ( 1 , 2 , 1-1 -2 ) are the barycentric coordinates of the unit parametric triangle, Ė are denoted as the corresponding two-ring neighboring control points for the interest triangle (See Fig. 3.2(a)), and is the quartic box-splines basis functions (cf. [35]). Any patch with irregular control points whose valence is unequal to six can also be parameterized with a box spline. In the procedure of Loop subdivision (see Fig. 3.2(b) and (c)), only regular points are generated through the original coarse triangular mesh divided by quartered separating, therefore we can see that all newly generated control points has valence of six. Obviously, it will create increasingly regular patches in the procedure of repeated Loop subdivision. The number of irregular control points of the mesh is not changed, therefore their influence region converges to the limit of zero. On the basis of this observation, a fast parameterization method for the Loop subdivision surface is proposed (cf. [35]), which agrees with the rule of box-spline regular patches.
In this subsection, we introduce the smoothness and convergence properties for Loop subdivision surfaces. For each control point of the mesh, we only consider its one-ring adjacent control points. Denote an assembly matrix
where the matrix Ĕ 0 is ( + 1) × 3, is the number of adjacent edges of the control vertex Ė 0 , and the superscript 0 is the subdivision level. We can compute the matrix Ĕ at the th subdivision level by use of
which means that the initial vector of coordinates Ĕ 0 is updated by times multiplication of the subdivision matrix ď. The fast algorithm about the subdivision matrix ď can be referred to as [35], where it depends only on the valence of the control point and is unrelated to the control point coordinates.
We can achieve a limit surface when the subdivision level ³ @, and the explicit position for every control point, which is described in the following Lemma 3.1. The limit surface of Loop subdivision arrive at 2 everywhere except at the extraordinary points where it is 1 .
as the subdivision time ³ @ (see [34] for the proof).
We consider the Loop subdivision surface, to be more precise, the subdivision shape functions equivalent to quartic box splines, for the discretization scheme for the weak form of the problem. Hence, it is imperative to utilize smooth, or more technically 2 -conforming, shape functions to achieve a conforming finite element discretization. On the limit surface evaluated using (3.4), the parametric coordinates like at Gaussian quadrature points refer to the unit parametric triangle on the control patch, therefore we need establish a mapping from the unit parametric triangle to the control surface patch where we perform the actual evaluation, see also [56,57] .
First, based on the linearization approach of the IEQ approach [49-51,55-57], we define a local-type auxiliary variable as
where the two positive constants and are specified as follows. The negative term -2 2 2 can be always bounded from below for any , because of the existence of the fourth-order polynomial ( ), when we apply a positive constant here so as to make the sum of all terms in the square root is always positive.
Using the new variable , one can reformulate the Cahn-Hilliard-Darcy equations as follows:
)
where the new term ( ) is given as
Second, based on the ZEC property of the surface tension and advection term, we introduce another nonlocal auxiliary variable ( ), and design its associated ODE system, that reads as
(3.12)
Obviously, from (2.18), we get = 0. Using the initial condition | =0 = 1, we get the exact solution of (3.12) is ( ) 1. Therefore, with the help of the new variables and , we rewrite the system (2.9)-(2.12) into the following equivalent form
∇ ç ē = 0, (3.17)
The initial conditions read as Remark 3.2. The above scheme differs significantly from our previous work in [70], where we used a combination of SAV and ZEC methods to discretize the system and consider only temporal semi-discretization without considering the spatial discretization. In the current numerical scheme constructed above, we use the EIEQ with the ZEC method for the temporal discretization and consider the spatial discretization on the complex surface. The above-developed fully-discrete scheme not only ensures all desired time marching properties-fully-decoupled, second-order temporal accuracy, linear, and unconditional energy stability-but also operates on complex manifolds of arbitrary shape and topology, a significant difference from the temporal semi-discrete scheme detailed in [70]. For complex manifolds of arbitrary shape and topology, the difficulty of parameterizing differential geometric operators on manifolds is encountered since there is no explicit expression for the surfaces. Although only the necessary operators are listed in Section 2.1, the fact that a large number of differential geometric operators on manifolds is very complex may even form a systematic research area. In addition, if the classical finite element method is used for the spatial discretization of manifolds, it will result in large geometric errors, which can be avoided by the IGA method used in this paper. Therefore, integrating the spatial discretization functionality provided by the IGA framework with multiple temporal discretization methods, and ultimately forming a unified and comprehensive full discretization framework, is the focus of this study.
The new transformed system (3.13)-(3.18) also retains the law of energy dissipation, which is described as the following theorem.
Theorem 3.1. The transformed equivalent system (3.13)-(3.18) holds the law of the energy dissipation as
where
Proof. By taking the 2 inner product of (3.13) with , and using the integration by parts, we get
By taking the 2 inner product of (3.14) within 2 space, and using the integration by parts, we get
By taking the 2 inner product of (3.15) with 2 , we obtain
By taking the 2 inner product of (3.16) with ē, and using integration by parts and applying (3.17
Now, we are ready to design the fully discrete scheme for solving the transformed system (3.13)- (3.18). We denote the time step size as > 0 and = for 0 d d = [ ∕ ]. Take the test functions ℎ , ℎ , ℎ , ℎ , ℎ * 2 (ā ℎ ), which is our subdivision-based IGA finite element space described in Section 3.1.
Suppose ℎ , ℎ , ℎ , , ē ℎ , ē ℎ for = , -1 and e 1 are known, we compute +1 ℎ ,
ℎ by the following second-order difference formula (BDF2):
)
where
(3.34) Remark 3.3. In the above-designed algorithm, a straightforward discretization approach has been utilized. Specifically, the linear component is treated implicitly, while the nonlinear component, which involves the nonlocal variable , is handled through a combination of explicit and implicit methods. This means that, in addition to the implicit treatment of , all other components are processed explicitly. The decoupling method for the above scheme, which appears apparently to be of the coupled type, has to take full advantage of the nonlocal character of the Q variable. Next, we will prove the energy stability of the proposed scheme, followed by the implementation details for decoupling.
Remark 3.4. Concerning the incompressible fluid dynamical system, we employ a projection-based method. Such an algorithm can decouple the pressure calculation from that of the velocity field. It is demonstrated that the velocity field and the pressure can achieve second-order and first-order accuracy in time, respectively, where rigorous error estimates were given in [71]. The reason for the loss of first-order accuracy of pressure is that an artificial homogeneous Neumann boundary condition is imposed on the pressure [72].
The intermediate velocity field is ē ℎ , and the final velocity field is ē ℎ , which can be shown to satisfy the divergence-free condition in a discrete sense, i.e., by taking the 2 inner product of (3.33) with ∇ ℎ , ℎ * ℎ , we obtain
Hence, from (3.32), we derive
We show the unconditional energy stability of the discrete scheme (3.27)- (3.33) in the following theorem. We will use the following two identities repeatedly:
(3.37)
Q. Pan et al. where, for an integer e 0, the discrete energy ̃ is defined as
(3.40)
Proof. By taking ℎ = 2 +1 ℎ in (3.27), we get
where (3.38) is used. By multiplying (3.31) with 2 +1 and using (3.38), we obtain
We combine (3.41)-(3.44) to get the following
(3.45)
By taking the inner product of (3.30) with 2 ē +1 ℎ in the 2 space, we get
We rewrite (3.33) as
Taking the 2 inner product of the above equality with ē ℎ with = + 1, , -1, we derive
and
(3.50)
We rewrite (3.33) as
Taking the 2 inner product of the above equation with itself and multiply the result with 3 2 , we derive
We rewrite (3.33) as
By taking the 2 inner product of the above equation with 3ē +1 ℎ , using (3.36) and (3.37), we obtain
Then by integrating (3.46), (3.50), (3.52) and (3.54), we get
By combining (3.45) and (3.55), and using (3.38), we obtain
(3.56)
Finally, we obtain the desired result (3.39) after dropping some positive terms of (3.56). ¦
In this section, we introduce a nonlocal splitting method for obtaining the decoupling implementation process of the proposed scheme (3.27)- (3.31).
Step 1: Since +1 is a nonlocal variable, we use it to split the local variables +1 ℎ , +1 ℎ and +1 ℎ into the linear combination forms, that read as
(3.57)
By applying the linear form given in (3.57), we can split the scheme (3.27)-(3.29) as follows:
Hence, in this step, one compute (3.58)-(3.61) to obtain +1 ℎ , +1 ℎ , +1 ℎ for = 1, 2.
Step 2: we use +1 to split ē +1 ℎ into the linear combination form that reads as
By applying the linear form given in (3.62), we can split the scheme (3.30) as follows, (
)
)
where 1 and 2 are given as
).
(3.66)
Note that we must establish the solvability of (3.65) by demonstrating that 3 2 -2 0, as follows. First, by taking ℎ = +1 2ℎ and ℎ = -3 2 +1 2ℎ in (3.59), we get
).
(3.67)
Finally we choose ℎ = ē +1 2ℎ in (3.64) to achieve Remark 3.5. In this study, we exclusively focus on conducting a stability analysis for the proposed scheme. It is worth noting the extensive research available on error estimates for the Cahn-Hilliard-Darcy system within regular regions (refer to [73][74][75][76]). However, to the best of the authors' knowledge, there is currently no existing relevant work on error estimates specifically related to the numerical scheme of the Cahn-Hilliard-Darcy system on manifolds. Addressing this gap will be a focus of our future research directions.
In this section, we implement some numerical simulations to test the energy stability, and verify the convergence rates and effectiveness of our developed numerical scheme (3.27)- (3.33), where various benchmark phenomena of the spinodal decompositions and the Saffman-Taylor fingering stability induced by the rotation force are carried out. A robust GMRES iterative solver is applied to solve the final systems generated from the variational format because these systems are highly sparse, where a small enough tolerance is set to obtain the proper convergence.
In the first numerical test, we verify the convergence rates for our proposed numerical scheme, where the computed domain is set as a spherical surface with radius 1,
The initial conditions are set as { 0 = cos( ) cos( ) cos( ), 0 = sin( ) sin( ) sin( ), ē 0 = (cos( ) sin( ) cos( ), sin( ) cos( ) cos( ), -cos( ) sin( ) sin( )) . We use four sets of refined meshes (ā ℎ , = 1, & , 4) generated by Loop subdivision, where ℎ denote the th level surface mesh size, with ℎ 1 = 0.2543, ℎ 2 = 0.1189, ℎ 3 = 0.0622, ℎ 4 = 0.0326. Since the subdivision-based IGA method has been proved to be secondorder accurate in [25], we set the time step sizes to = ℎ where is the given constant with * (0, 1). Note that the exact solutions are unknown, therefore the approximate solutions obtained by employing a sufficiently fine mesh size are assumed to be the exact solutions in order to observe numerical errors. In Fig. 4.1, we plot the 2 errors at the time = 2.56e-1 for all unknown variables , ē, and . The second order convergence rate for and ē, and the first order convergence rate for can be observed, which are in good agreement with the theoretical expected rate.
In this example, we perform the numerical simulations of phase separation (spinodal decomposition) on a sphere using the developed scheme, which is a standard numerical test of the Cahn-Hilliard related system in which the initial conditions are generally set to be uniformly perturbed within a some scope of binary mixtures. The system can change from a homogeneous situation to a stable dissociation of the two phases as a result of an intrinsic increase in concentration fluctuations.
The surface domain ā 1 is set as a closed sphere with radius 2, which reads as
In this example, we simulate the viscous finger-pattern instability, i.e., the Saffman-Taylor instability, which is the signature problem of the interfacial instability exhibited by a two-phase Darcy fluid system in a Hele-Shaw cell. There are typically two methods to achieve viscous fingering, one is to use a less viscous fluid to displace another more viscous fluid between narrowly spaced plates, and the other is to introduce the rotational forces into the system, which is adopted here. At the initial state, a droplet is placed in a rotating Hele-Shaw cell, and under the centrifugal force of the rotating motion, the interface of the droplet develops a vivid rotating finger-like structure, transitioning from a droplet-like structure to filamentous arms with extended arms, or a membrane-like skeleton structure with nearly constant fingertip width [77][78][79][80][81].
To impose a rotational force, we replace the momentum equation (2.11) with the following form:
where a rotational force is added on the right hand side, and Č is the outward normal direction at the surface, and are two pre-assigned positive constants, and Đ is a given vector.
We carry out two simulations on two surfaces, the first of which is a spherical domain with the radius 1, that reads as
It is discretized with 2 097 152 Loop subdivision limit subdivision patches by use of 1 048 578 control vertices, and the span of the vertex valence is 4 to 6. The initial conditions are set as 0 = tanh( 0 - In Fig. 4.4, we show the snapshots of the variable at various times, where red-blue color map is used. We can see that over time, the originally smooth interface produces tiny tentacles that slowly grow and elongate, forming the typical viscous fingering pattern. To detect the effect of different model parameters on the obtained fingering pattern, in Fig. 4.5, we vary the surface tension parameter with = 5e-3, 1.5e-3, 5e-4, 5e-5 for comparison of the differences in the fingering patterns obtained. When the surface tension parameter is large, the number of finger generated is small and the shape is thick. When the parameter is small, the number of finger becomes more and more, and the shape becomes thin and long, consistent to other experimental/numerical results performed on the flat surface [77,78,81].
We conduct a second simulation with a similar rotation in another, but with the computational region switched to the more complex surface of a human head, that reads as In Fig. 4.6, snapshots of the variable at various times are shown. For clarity, we show two different views for each 3D head at any given time. Similar to the previous simulation, tiny fingers initially appear on the interface and then slowly grow longer and denser across the entire surface of the head, forming viscous fingering patterns densely spread across the computed curved head surface.
In this paper, we study a highly nonlinear and coupled system, namely, the Cahn-Hilliard-Darcy phase-field model for two-phase flow in the Hele-Shaw cell, on a complex curved surface. We establish a novel numerical approach within the IGA framework that integrates subdivision surface technology with a number of efficient time marching methods, including explicit-IEQ linearization and ZEC decoupling. Subdivision technology can bring high smoothness and hierarchic refineability, precisely representing the convoluted surface with arbitrary topology. Effective time marching and decoupling methods enable us to discretize nonlinear terms without introducing additional nonlinear iterations and any form of coupling computation in a simple and straightforward manner. The optimal combination of time marching and space discretization permits the development of numerical methods with all of the desired properties: geometric exactness of the domains of interest, linearity, second-time accuracy in time, decoupling, and unconditional energy stability. We rigorously demonstrate the unconditional energy stability and provide numerous numerical examples of complex surface models, in which the accuracy, stability, and some benchmarking fingering pattern examples produce a captivating visual effect that numerically validates the robustness of our established numerical method.