A ternary mixture model with dynamic boundary conditions

Shuang Liu; Yue Wu; Xueping Zhao; Shuang Liu; Yue Wu; Xueping Zhao

doi:10.3934/mbe.2024091

Mathematical Biosciences and Engineering

2024, Volume 21, Issue 2: 2050-2083. doi: 10.3934/mbe.2024091

Previous Article Next Article

Research article Special Issues

A ternary mixture model with dynamic boundary conditions

1.
Department of Mathematics, University of North Texas, 1155 Union Circle, Denton, Texas 76203-5017, USA
2.
Department of Mathematical Sciences, University of Nottingham Ningbo China, Taikang East Road 199, Ningbo 315100, China

Academic Editor: Jia Zhao

Received: 31 October 2023 Revised: 14 December 2023 Accepted: 21 December 2023 Published: 08 January 2024

The influence of short-range interactions between a multi-phase, multi-component mixture and a solid wall in confined geometries is crucial in life sciences and engineering. In this work, we extend the Cahn-Hilliard model with dynamic boundary conditions from a binary to a ternary mixture, employing the Onsager principle, which accounts for the cross-coupling between forces and fluxes in both the bulk and surface. Moreover, we have developed a linear, second-order and unconditionally energy-stable numerical scheme for solving the governing equations by utilizing the invariant energy quadratization method. This efficient solver allows us to explore the impacts of wall-mixture interactions and dynamic boundary conditions on phenomena like spontaneous phase separation, coarsening processes and the wettability of droplets on surfaces. We observe that wall-mixture interactions influence not only surface phenomena, such as droplet contact angles, but also patterns deep within the bulk. Additionally, the relaxation rates control the droplet spreading on surfaces. Furthermore, the cross-coupling relaxation rates in the bulk significantly affect coarsening patterns. Our work establishes a comprehensive framework for studying multi-component mixtures in confined geometries.

Keywords:

Citation: Shuang Liu, Yue Wu, Xueping Zhao. A ternary mixture model with dynamic boundary conditions[J]. Mathematical Biosciences and Engineering, 2024, 21(2): 2050-2083. doi: 10.3934/mbe.2024091

Related Papers:

[1]	Lorena Bociu, Giovanna Guidoboni, Riccardo Sacco, Maurizio Verri . On the role of compressibility in poroviscoelastic models. Mathematical Biosciences and Engineering, 2019, 16(5): 6167-6208. doi: 10.3934/mbe.2019308
[2]	Gangqi Han, Bing Wang, Mengli Jia, Shuxin Ding, Wenxuan Qiu, Yuxuan Mi, Zhimei Mi, Yuhao Qin, Wenxing Zhu, Xinli Liu, Wei Li . Optimization and evaluation of resveratrol amorphous solid dispersions with a novel polymeric system. Mathematical Biosciences and Engineering, 2022, 19(8): 8019-8034. doi: 10.3934/mbe.2022375
[3]	Rachel Clipp, Brooke Steele . An evaluation of dynamic outlet boundary conditions in a 1D fluid dynamics model. Mathematical Biosciences and Engineering, 2012, 9(1): 61-74. doi: 10.3934/mbe.2012.9.61
[4]	Longting Ding, Yuan Li, Zhanchuang Han, Mengyuan Zhang, Xuancang Wang, Lu He . The effect of RAP content on fatigue damage property of hot reclaimed asphalt mixtures. Mathematical Biosciences and Engineering, 2024, 21(2): 3037-3062. doi: 10.3934/mbe.2024135
[5]	József Z. Farkas, Peter Hinow . Physiologically structured populations with diffusion and dynamic boundary conditions. Mathematical Biosciences and Engineering, 2011, 8(2): 503-513. doi: 10.3934/mbe.2011.8.503
[6]	Sunmi Lee, Chang Yong Han, Minseok Kim, Yun Kang . Optimal control of a discrete-time plant–herbivore/pest model with bistability in fluctuating environments. Mathematical Biosciences and Engineering, 2022, 19(5): 5075-5103. doi: 10.3934/mbe.2022237
[7]	Rinaldo M. Colombo, Mauro Garavello . Stability and optimization in structured population models on graphs. Mathematical Biosciences and Engineering, 2015, 12(2): 311-335. doi: 10.3934/mbe.2015.12.311
[8]	J. Leonel Rocha, Sandra M. Aleixo . An extension of Gompertzian growth dynamics: Weibull and Fréchet models. Mathematical Biosciences and Engineering, 2013, 10(2): 379-398. doi: 10.3934/mbe.2013.10.379
[9]	A. Q. Khan, M. Tasneem, M. B. Almatrafi . Discrete-time COVID-19 epidemic model with bifurcation and control. Mathematical Biosciences and Engineering, 2022, 19(2): 1944-1969. doi: 10.3934/mbe.2022092
[10]	Wenzhang Huang, Maoan Han, Kaiyu Liu . Dynamics of an SIS reaction-diffusion epidemic model for disease transmission. Mathematical Biosciences and Engineering, 2010, 7(1): 51-66. doi: 10.3934/mbe.2010.7.51

Abstract

1. Introduction

The physics of multi-component mixtures, including spontaneous phase separation, nucleation and growth and coarsening have been well-studied both theoretically and experimentally. The question of how the multi-component, multi-phase system interacts with solid walls in confined geometries has gained much attention due to its wide applications in life science and engineering, e.g., membrane-less organelle formation and wetting on the membrane-bound organelles in living cells^[1,2,3]; micro-fluid close to a solid interface in microfluidic devices ^[4]; thin films of polymer blends in slab geometry ^[5,6,7] and so on.

Numerous methods have been developed to study interface problems and their interaction with solid walls. For instance, the lattice Boltzmann model has garnered considerable attention as a tool to facilitate the simulation of multi-phase and multi-component systems with complex boundary conditions, as demonstrated by several studies ^{[8,9,10,11,12]}. Additionally, the exploration of wetting properties in ternary mixture models has been a significant focus ^[13,14,15]. For example, Liang et al. ^[16] formulated a wetting boundary condition for ternary fluids interacting with a solid substrate, applying this to the lattice Boltzmann model. The volume-of-fluid (VOF) technique is particularly effective in scenarios involving solid obstacles, as it integrates the contact angle by applying appropriate boundary conditions to the fluid-function's gradient at solid interfaces. Therefore, VOF has been extensively used in fluid dynamics analyses within complex porous structures ^{[17,18,19,20,21]}. However, among all of these methods, the phase-field method stands out as the most effective tool for modeling and simulating multi-phase systems. Its inherent flexibility in handling complex topological changes in the phases, coupled with its ability to seamlessly integrate with various boundary conditions, makes it a superior choice to address the complexities of interface problems in multi-component systems.

The Cahn-Hilliard equation is a fundamental phase-field model. Since it was first proposed in the field of materials science to describe the pattern formation evolution of microstructures during the phase separation process in binary alloys ^[22,23], the Cahn-Hilliard equation and its variants have been successfully applied to model a wide variety of segregation-like phenomena in science; see, for instance, ^{[24,25,26,27,28,29,30,31]} and the references therein. In this study, we denote the domain as $\Omega$ and the boundary of the domain $\Gamma$ . We introduce $\phi$ , which represents the volume fraction of a specific component within the mixture. $\phi \in (0, 1)$ is dimensionless and quantifies the proportion of the total volume occupied by a given component, providing insight into the concentration and distribution of that component within the system. The total free energy of the system is given as

$\begin{equation} F = \int_{\Omega} \Big[ f_b(\phi) + \frac{\kappa}{2} |\nabla \phi|^2\Big] d\boldsymbol{x}, \end{equation}$

(1.1)

where $f_b(\phi)$ is the bulk free-energy density function and the term $\frac{\kappa}{2} |\nabla \phi|^2$ represents the interfacial energy between phases. The Cahn-Hilliard equation has the following format:

$\begin{align} \partial_t \phi = \nabla M \cdot \nabla \mu \, , \end{align}$

(1.2)

where mobility $M$ can be either a constant or concentration-dependent function and $\mu$ is the chemical potential of the component $\phi$ . The boundary conditions include the following two categories:

1) periodic boundary conditions,

2) physical boundary conditions, e.g., homogeneous Neumann boundary conditions.

However, in some systems, e.g., phase separation in confined geometries, the condensate will interplay with boundaries; thus more generic boundary conditions, which describe the interaction between wall and condensates, are introduced. In this work, we will focus on the study of generic boundary conditions.

The influence of boundaries (solid walls) on the phase separation process of binary mixtures has attracted a lot of attention from scientists. For instance, Fischer et al. ^[27] introduced the following dynamic boundary conditions in 1997:

$\begin{align} \partial_n \mu|_{\Gamma} & = 0, \end{align}$

(1.3a)

$\begin{align} \partial_t \phi|_{\Gamma} & = - \Gamma_s \Big[ \kappa \, \partial_{n} \phi + {f_s}' - \kappa_{\Gamma} \Delta_{\Gamma} \phi_{\Gamma} \Big], \end{align}$

(1.3b)

where $f_s$ is the surface free-energy density, which represents the short-range interaction between mixture components and the solid wall, $\Delta_{\Gamma}$ stands for the Laplace–Beltrami operator on the boundary surface $\Gamma$ , $\Gamma_s$ defines a surface kinetic coefficient and the term $\kappa \, \partial_n \phi$ is due to the surface contribution that comes from the variation of the bulk free energy $f_b$ . The boundary (1.3a) is simply the condition that no current can flow through the surface, while the boundary (1.3b) can be derived by requiring the system to tend to minimize its surface free energy as follows:

$\begin{align} F_s = \int_{\Omega} \Big[ \frac{\kappa_{\Gamma}}{2} |\nabla_{\Gamma} \phi_{\Gamma}|^2 + f_s\Big]. \end{align}$

(1.4)

Boundary conditions of similar form as (1.3) have also been derived from a semi-infinite Ising model with Kawasaki spin exchange dynamics ^[32,33]. In 2011, Goldstein et al. derived a Cahn-Hilliard model in a domain with non-permeable walls ^[29]. In this model, they assumed that the total mass in the bulk and on the boundary, i.e.,

$\begin{align} \int_\Omega \phi(\boldsymbol{x}) d\boldsymbol{x} + \int_{\Gamma} \phi(\boldsymbol{x}) dS, \end{align}$

(1.5)

is conserved. Its boundary conditions are as follows:

$\begin{align} \partial_t \phi|_{\Gamma} & = \nabla_{\Gamma} M \cdot \nabla_{\Gamma} \mu_{\Gamma} - \beta M \partial_n \mu, \end{align}$

(1.6a)

$\begin{align} \mu|_{\Gamma} & = \beta \Big[ \kappa \, \partial_{n} \phi + {f_s}' - \kappa_{\Gamma} \Delta_{\Gamma} \phi_{\Gamma} \Big] . \end{align}$

(1.6b)

In addition, they proved the existence and uniqueness of weak solutions and studied their asymptotic behavior as time goes to infinity.

Liu and Wu ^[30] introduced a Liu-Wu model in 2019 with non-flux boundary conditions for chemical potential and Cahn–Hilliard-type gradient flow on the boundary, i.e.,

$\begin{align} \partial_n \mu & = 0, \end{align}$

(1.7a)

$\begin{align} \partial_t \phi|_{\Gamma} & = \nabla_{\Gamma} M \cdot \nabla_{\Gamma} \Big[ \kappa \, \partial_{n} \phi + {f_s}' - \kappa_{\Gamma} \Delta_{\Gamma} \phi_{\Gamma} \Big] . \end{align}$

(1.7b)

In 2021, Patrik Knopf and collaborators proposed a general model, i.e., the KLLM model ^[31]:

$\begin{align} L\partial_{n} \mu & = -\mu|_{\Gamma} + \beta \Big[ \kappa \, \partial_{n} \phi + {f_s}' - \kappa_{\Gamma} \Delta_{\Gamma} \phi_{\Gamma} \Big] , \end{align}$

(1.8a)

$\begin{align} \partial_t \phi|_{\Gamma} & = \nabla_{\Gamma} M \cdot \nabla_{\Gamma} \Big[ \kappa \, \partial_{n} \phi + {f_s}' - \kappa_{\Gamma} \Delta_{\Gamma} \phi_{\Gamma} \Big]-\beta M \partial_{n} \mu. \end{align}$

(1.8b)

Garcke et al. ^[34] performed nonlinear analysis to evaluate the long-time dynamics of the Cahn–Hilliard equation with kinetic rate dependent dynamic boundary conditions for the KLLM model. We refer the reader to ^[35] for a review of the model and analysis.

Correspondingly, there are numerical studies for binary mixture models with dynamic boundary conditions. For example, Kenzler et al. ^[28] proposed an implicit numerical scheme to solve the Cahn–Hilliard equation with (1.3). The authors discretized the solution in space by applying the finite-difference method to edge points. A method of stability analysis was derived and a variable time-stepping strategy was employed to simulate the system over a long period of time. Furthermore, the numerical scheme is conditionally gradient-stable, which means that the free energy decreases monotonously in time. In 2017, Fukao et al. ^[36] constructed a structure-preserving finite-difference scheme for the Cahn–Hilliard equation with dynamic boundary conditions in the one-dimensional case for Goldstein model (1.6). Unlike the space discretization in ^[28], the authors applied the finite-difference method to center points. The existence of the solution and the error estimate were also obtained. In addition, the laws of mass conservation and energy dissipation were satisfied at the discrete level. Meng et al. ^[37] developed a second-order stabilized semi-implicit scheme for Liu-Wu model (1.7). The corresponding energy stability and convergence analysis of the scheme were derived theoretically.

Besides the binary mixture models, phase-field models for ternary mixture systems have been developed due to their frequent real-world applications, as indicated in the literature ^[38,39,40]. Wetting behavior^[41,42] has also been proposed based on the phase-field models. For example, a phase-field model for multi-component Cahn-Hilliard systems in complex domains was developed in ^[41] by considering contact angle or no mass flow boundary conditions. Building on this, Yang et al. ^[42] extended the geometric formulation of wetting conditions by using the weighted contact angles defined within the implementation of wetting conditions, following the concept of the diffusive interface ^[43]. This innovative phase-field model efficiently describes the dynamic behavior of compound droplets in contact with solid objects.

However, we note that dynamic boundary conditions are missing in the aforementioned multi-component phase-field models. Moreover, to the best of our knowledge, a ternary mixture model with the force-flux cross-coupling relations has not been explored yet. In summary, based on the above discussion, we focus on the following aims in this work:

1) Derive a model for a ternary mixture with the force-flux cross-coupling relations in a confined geometry with dynamic boundary conditions based on the Onsager principle and irreversible thermodynamics ^{[44,45,46,47,48,49]}.

2) Develop a linear second-order unconditionally energy-stable numerical scheme for the derived model based on the invariant energy quadratization (IEQ) method ^{[50,51,52,53]}.

3) Investigate equilibrium and non-equilibrium properties of the system based on the numerical simulations.

Section 2 details the derivation of a kinetic model. This model, which addresses a ternary mixture within confined geometries and its interaction with adjacent solid walls, is formulated by using principles of irreversible thermodynamics. It features a conservation of total mass and a decrease in total energy over time. The mathematical model undergoes an equivalent reformulation based on the IEQ method in Section 3, where we also introduce a second-order, linear, unconditional energy-stable numerical scheme for the revised model. Section 4 examines how spontaneous phase separation and surface wettability are influenced by various model parameters through the use of numerical studies. The paper concludes in Section 5, where we summarize our findings. For clarity, we list abbreviations and notations in Tables 1 and 2.

Table 1. Abbreviations.

Abbreviations	Explanations
DBC	Dynamic boundary condition
IEQ	Invariant energy quadratization

| Show Table

DownLoad: CSV

Table 2. Notations.

Notations	Explanations
$\chi_{ij}$	interaction strength between the components $i$ and $j$
$\gamma$	coupling interaction strength between components on the surface
$\Gamma$	boundary of domain
$\Gamma_i$	relaxation parameter of dynamic boundary condition, $i$ =1, 2, 12
$\Gamma_s$	surface kinetic coefficient
$\kappa_B$	Boltzmann constant
$\kappa_i$ , $\kappa_{i\Gamma}$	gradient coefficients associated with interface free energies, $i$ =1, 2, 12
$\mu_i$	chemical potential in the bulk for components $i$ =1, 2
$\mu_{i\Gamma}$	chemical potential on the surface for components $i$ =1, 2
$\Delta_\Gamma$	Laplace–Beltrami operator on the boundary surface
$\nu$	molecular volume of component 3
$\nu_i$	molecular volume of component $i$ , $i$ =1, 2
$\Omega$	domain
$\omega_i$	internal free energy coefficient for component $i$ in the bulk
$\phi_i$	volume fraction of three components, $i$ =1, 2, 3
$\phi_{i\Gamma}$	volume fraction of a component, $i$ =1, 2
$D_i$	diffusion coefficients, $i$ =1, 2, 12
$f_b$	bulk free-energy density function
$f_s$	surface free-energy density function
$g_i$	molecule-molecule interaction strength of each component near the wall
$h_i$	interaction strength proportional to the component volume fraction, $i$ =1, 2
$J_i$	mass flux, $i$ =1, 2
$M_i$	mobility function, $i$ =1, 2, 12
T	temperature of the isotropic system

| Show Table

DownLoad: CSV

2. Mathematical model derivation

In this section, we will extend the binary mixture model with dynamic boundary conditions to a ternary mixture model, based on the Onsager principles^{[44,45,46,47,48]}. We denote the volume fractions of three components as $\phi_1$ , $\phi_2$ and $\phi_3$ , where $\phi_3 = 1-\phi_1 - \phi_2$ based on the incompressibility assumption. The total free energy of the system involves two contributions: one is the free energy in the bulk $\Omega$ , and another on the surface $\Gamma$ , i.e.,

$\begin{align} E = E_{\text{bulk}} + E_{\text{surf}}, \end{align}$

(2.1)

where the bulk free energy $E_{\text{bulk}}$ and surface free energy $E_{\text{surface}}$ are as follows:

$\begin{align} E_{\text{bulk}} = \int_{\Omega} \Big[ f_b(\phi_1, \phi_2) + \frac{\kappa_1}{2}|\nabla \phi_1|^2 + \frac{\kappa_2}{2}|\nabla \phi_2|^2 + \kappa_{12} \nabla \phi_1 \cdot \nabla \phi_2 \Big], \end{align}$

(2.2)

$\begin{align} E_{\text{surf}} = \int_{\Gamma} \Big[ f_s(\phi_{1 \Gamma}, \phi_{2 \Gamma}) + \frac{\kappa_{1\Gamma}}{2}|\nabla_{\Gamma} \phi_{1 \Gamma}|^2 + \frac{\kappa_{2\Gamma}}{2}|\nabla_{\Gamma} \phi_{2 \Gamma}|^2 + \kappa_{12\Gamma}\nabla_{\Gamma}\phi_{1\Gamma} \cdot \nabla_{\Gamma} \phi_{2\Gamma} \Big]. \end{align}$

(2.3)

$\kappa_{i}$ , $\kappa_{i\Gamma}$ , $i = 1, 2, 12$ denote the gradient coefficients associated with interface free energies. For simplicity, we set $\kappa_{12} = 0$ and $\kappa_{12\Gamma} = 0$ in the following work. $\nabla_{\Gamma}$ denotes the gradient operator on the surface. The bulk free-energy density function $f_b(\phi_1, \phi_2)$ and surface free-energy density function $f_s(\phi_{1 \Gamma}, \phi_{2 \Gamma})$ are defined as the Flory-Huggins free energy and a second-order polynomial based on application in the polymer field ^[54]:

$\begin{align} f_b(\phi_1, \phi_2) & = \frac{k_B T}{\nu} \Big[ \frac{\phi_1}{n_1} \ln \phi_1 + \frac{\phi_2}{n_2} \ln \phi_2 + (1-\phi_1 - \phi_2) \ln(1-\phi_1 - \phi_2) + \chi_{12} \phi_1 \phi_2 \\ & \quad + \chi_{13} \phi_1 (1-\phi_1 -\phi_2) + \chi_{23} \phi_2 (1-\phi_1 - \phi_2) + \omega_1 \phi_1 + \omega_2 \phi_2 \Big], \end{align}$

(2.4)

$\begin{align} f_s(\phi_{1\Gamma}, \phi_{2\Gamma)} & = \frac{k_B T}{\nu_s} \Big[ h_1 \phi_{1\Gamma} + g_1 \phi_{1\Gamma}^2 + h_2 \phi_{2\Gamma} + g_2 \phi_{2\Gamma}^2 + \gamma \phi_{1\Gamma} \phi_{2\Gamma} \Big], \end{align}$

(2.5)

where $\nu_i = n_i \nu$ represents the molecular volume of component $i$ , $i = 1, 2$ , and $\nu$ represents the molecular volume of component 3. For simplicity, we assume the molecular volumes of each component are equivalent, i.e., $\nu_1 = \nu_2 = \nu$ and $n_1 = n_2 = 1$ . $k_B$ is the Boltzmann constant and $T$ is the temperature of the isotropic system. The term $\chi_{ij} \phi_i \phi_j$ refers to the interaction strength between the component $i$ and the component $j$ . $\omega_i \phi_i$ denotes the internal free energy for component $i$ in the bulk $\Omega$ . On the surface, $h_i\phi_i$ , $i = 1, 2$ is the linear order term, which denotes the interaction strength as is proportional to the component volume fraction. The quadratic $g_i \phi_i^2$ depicts the molecule-molecule interactions of each components near the wall, while $\gamma \phi_1 \phi_2$ describes the coupling interaction between components on the surface.

We assume that there is no chemical reaction in the system, that is, each component of the mixture conserves the total mass. We write the conservation laws as follows:

$\begin{align} \partial_t \phi_1 +\nabla \cdot \boldsymbol{J_1} & = 0, \end{align}$

(2.6)

$\begin{align} \partial_t \phi_2 + \nabla \cdot \boldsymbol{J_2} & = 0. \end{align}$

(2.7)

2.1. Irreversible thermodynamics

We consider a closed system in this work. According to irreversible thermodynamics, the entropy product rate is inversely proportional to the energy change rate.

$\begin{align} - T \partial_t S & = \frac{dE}{dt} = \frac{\partial E}{\partial \phi_1} \frac{\partial \phi_1}{\partial t}+ \frac{\partial E}{\partial \phi_2} \frac{\partial \phi_2}{\partial t} + \frac{\partial E}{\partial \phi_{1 \Gamma} }\frac{\partial \phi_{1 \Gamma}}{\partial t} + \frac{\partial E}{\partial \phi_{2 \Gamma}} \frac{\partial \phi_{2 \Gamma}}{\partial t} \\ & \quad + \frac{\partial E}{\partial \nabla \phi_1} \frac{\partial \nabla \phi_1}{\partial t}+ \frac{\partial E}{\partial \nabla \phi_2} \frac{\partial \nabla \phi_2}{\partial t} + \frac{\partial E}{\partial \nabla_{\Gamma} \phi_{1 \Gamma} }\frac{\partial \nabla_{\Gamma} \phi_{1 \Gamma}}{\partial t} + \frac{\partial E}{\partial \nabla_{\Gamma} \phi_{2 \Gamma}} \frac{\partial \nabla_{\Gamma} \phi_{2 \Gamma}}{\partial t} \end{align}$

(2.8)

$\begin{align} & = \int_{\Gamma} dS \Big[ \boldsymbol{n}\cdot \kappa _1 \nabla \phi_1 + (\frac{\partial f_s}{\partial \phi_{1 \Gamma}} -\kappa_{1 \Gamma} \Delta_{\Gamma} \phi_1) \Big] \frac{\partial \phi_{1\Gamma}}{\partial t} \\ & \quad + \int_{\Gamma} dS \Big[ \boldsymbol{n}\cdot \kappa_2 \nabla \phi_2 + (\frac{\partial f_s}{\partial \phi_{2 \Gamma}} -\kappa_{2 \Gamma} \Delta_{\Gamma} \phi_2) \Big] \frac{\partial \phi_{2\Gamma} }{\partial t} \\ & \quad - \int_{\Gamma} dS \Big[\frac{\mu_1}{\nu_1} \boldsymbol{n} \cdot \boldsymbol{J_1} + \frac{\mu_2}{\nu_2} \boldsymbol{n} \cdot \boldsymbol{J_2}\Big] + \int_{\Omega} d\boldsymbol{x} \Big[\nabla \frac{\mu_1}{\nu_1} \cdot \boldsymbol{J_1} + \nabla \frac{\mu_2}{\nu_2} \cdot \boldsymbol{J_2}\Big], \end{align}$

(2.9)

where the chemical potentials in the bulk for components 1 and 2 are respectively defined as

$\begin{align} \mu_1 & = \nu_1 \frac{\delta E}{\delta \phi_1} = \nu_1 \Big[\frac{\partial f_b}{\partial \phi_1} - \kappa_1 \Delta \phi_1 \Big], \end{align}$

(2.10)

$\begin{align} \mu_2 & = \nu_2 \frac{\delta E}{\delta \phi_2} = \nu_2 \Big[\frac{\partial f_b}{\partial \phi_2} - \kappa_2 \Delta \phi_2 \Big]. \end{align}$

(2.11)

We obtain the conjugate fluxes and forces as follows:

$\begin{align} \boldsymbol{J_1} & \longleftrightarrow - \nabla\mu_1\, , & x \in \Omega\, , \end{align}$

(2.12a)

$\begin{align} \boldsymbol{J_2} & \longleftrightarrow - \nabla\mu_2\, , & x \in \Omega\, , \end{align}$

(2.12b)

$\begin{align} \partial_t \phi_{1 \Gamma} & \longleftrightarrow - \left( \boldsymbol{n}\cdot \kappa _1 \nabla \phi_1 + \frac{\partial f_s}{\partial \phi_{1 \Gamma}} -\kappa_{1 \Gamma} \Delta_{\Gamma} \phi_1 \right) \, , & x \in \Gamma \, , \end{align}$

(2.12c)

$\begin{align} \partial_t \phi_{2 \Gamma} & \longleftrightarrow - \left( \boldsymbol{n}\cdot \kappa _2 \nabla \phi_2 + \frac{\partial f_s}{\partial \phi_{2 \Gamma}} -\kappa_{2 \Gamma} \Delta_{\Gamma} \phi_2 \right) \, , & x \in \Gamma \, . \end{align}$

(2.12d)

Based on the Onsager principle, we define a symmetric positive definite mobility function between conjugate fluxes and forces in both the bulk $\Omega$ and the surface $\Gamma$ , i.e.,

$\begin{align} \boldsymbol{J_1} & = - M_1(\phi_1, \phi_2) \nabla \mu_1 - M_{12}(\phi_1, \phi_2) \nabla \mu_2 , \end{align}$

(2.13a)

$\begin{align} \boldsymbol{J_2} & = - M_{12}(\phi_1, \phi_2) \nabla \mu_1 - M_2(\phi_1, \phi_2) \nabla \mu_2 , \end{align}$

(2.13b)

$\begin{align} \partial_t \phi_{1 \Gamma} & = -\Gamma_1 \Big[\boldsymbol{n} \cdot \kappa_1 \nabla \phi_1 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}} \Big] - \Gamma_{12} \Big[\boldsymbol{n} \cdot \kappa_2 \nabla \phi_2 - \kappa_{2 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}} \Big] , \end{align}$

(2.13c)

$\begin{align} \partial_t \phi_{2 \Gamma} & = - \Gamma_{12} \Big[\boldsymbol{n} \cdot \kappa_1 \nabla \phi_1 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}} \Big] -\Gamma_2 \Big[\boldsymbol{n} \cdot \kappa_2 \nabla \phi_2 - \kappa_{2 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}}\Big] , \end{align}$

(2.13d)

where $\Gamma_1$ , $\Gamma_2$ and $\Gamma_{12}$ are the relaxation parameters of the dynamic boundary conditions. $M_{i}(\phi_1, \phi_2)$ and $M_{12}(\phi_1, \phi_2)$ are defined as

$\begin{align} M_1(\phi_1, \phi_2) & = m_1 \phi_1 (1 - \phi_1 - \phi_2), \end{align}$

(2.14)

$\begin{align} M_2(\phi_1, \phi_2) & = m_2 \phi_2 (1 - \phi_1 - \phi_2), \end{align}$

(2.15)

$\begin{align} M_{12}(\phi_1, \phi_2) & = m_{12} \phi_1 \phi_2 (1 - \phi_1 - \phi_2) \end{align}$

(2.16)

to remove the singularity from the Flory-Huggins free-energy density functions. $M_{i}(\phi_1, \phi_2)$ and $\Gamma_{i}$ , $i = 1, 2$ are the diagonal elements in the mobility matrix, depicting the direct effects of forces on the corresponding fluxes, while $M_{12}(\phi_1, \phi_2)$ and $\Gamma_{12}$ are the cross-coupling terms denoting the cross interplay of forces and fluxes. The cross-coupling relaxation rate refers to the rate at which perturbations in one variable are coupled to changes in another variable in the system. We note that our model is an extension of the binary model developed in ^[27]. Our work provides a general framework that can extend all existing binary mixture models ^{[29,30,31,35]} with dynamic boundary conditions to a ternary model.

Combining the definition of fluxes derived from the Onsager principle with the conservation law (2.6), we obtain the governing equations of the system as follows:

$\begin{align} \ \partial_t \phi_1 & = \nabla M_1(\phi_1, \phi_2) \cdot \nabla \mu_1 + \nabla M_{12}(\phi_1, \phi_2) \cdot \nabla \mu_2 , \end{align}$

(2.17a)

$\begin{align} \partial_t \phi_2 & = \nabla M_2(\phi_1, \phi_2) \cdot \nabla \mu_2 + \nabla M_{12}(\phi_1, \phi_2) \cdot \nabla \mu_1, \end{align}$

(2.17b)

with the following boundary conditions:

(2.18a)

$\begin{align} \partial_t \phi_{2 \Gamma} & = -\Gamma_2 \Big[\boldsymbol{n} \cdot \kappa_2 \nabla \phi_2 - \kappa_{2 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}}\Big] - \Gamma_{12} \Big[\boldsymbol{n} \cdot \kappa_1 \nabla \phi_1 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}} \Big], \end{align}$

(2.18b)

$\begin{align} \partial_{\boldsymbol{n}} \mu_1|_{\Gamma} & = 0, \end{align}$

(2.18c)

$\begin{align} \partial_{\boldsymbol{n}} \mu_2|_{\Gamma} & = 0. \end{align}$

(2.18d)

We denote the diffusion coefficients as $D_i = m_i k_B T$ , $i = 1, 2$ , and $D_{12} = m_{12} k_B T$ . In this work, we use the dynamic boundary conditions given by (2.18) in the $x$ direction and the periodic boundary condition for all variables in the $y$ direction.

Our model has the following properties:

Theorem 2.1. The total mass of each component in the domain $\Omega$ is conserved, i.e.,

$\begin{align} \int_{\Omega} \phi_i(x, y, t) d\boldsymbol{x} = \int_{\Omega} \phi_i(x, y, 0) d\boldsymbol{x}, \quad i = 1, 2. \end{align}$

(2.19)

proof: For $i = 1$ ,

$\begin{align} \frac{d}{dt}\int_{\Omega} \phi_1 dx & = \int_{\Omega} \partial_{t} \phi_1 dx = -\int_{\Omega} \nabla \cdot \boldsymbol{J_1} dx = \int_{\Gamma} \boldsymbol{J_1}\cdot\boldsymbol{n} dS \\ & = \int_{\Gamma} \Big[- M_1(\phi_1, \phi_2) \nabla \mu_1 - M_{12}(\phi_1, \phi_2) \nabla \mu_2\Big]\cdot\boldsymbol{n} dS = 0, \end{align}$

(2.20)

with the boundary conditions $\partial_{\boldsymbol{n}} \mu_1|_{\Gamma} = 0$ and $\partial_{\boldsymbol{n}} \mu_2|_{\Gamma} = 0$ . It is similar for $i = 2$ . Thus, (2.19) is obtained.

Theorem 2.2. The total free energy is dissipative, i.e.,

$\begin{align} \partial_t E = \partial_t (E_{\mathit{\text{surf}}} + E_{\mathit{\text{bulk}}}) \leq 0. \end{align}$

(2.21)

proof: For simplicity, we assume that $\nu_1 = \nu_2 = \nu$ , and we define

$\begin{align} \mu_{1\Gamma} = \boldsymbol{n} \cdot \kappa_1 \nabla \phi_1 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}}, \end{align}$

(2.22)

$\begin{align} \mu_{2\Gamma} = \boldsymbol{n} \cdot \kappa_2 \nabla \phi_2 - \kappa_{2 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}}. \end{align}$

(2.23)

The total energy dissipation rate is given by

$\begin{align} \partial_t E & = \int_{\Omega} dx \Big[ (\frac{\partial f_b}{\partial \phi_1} - \kappa_1 \Delta \phi_1)\partial_t \phi_1 + (\frac{\partial f_b}{\partial \phi_2} - \kappa_2 \Delta \phi_2 )\partial_t \phi_2 \Big] \\ & \quad + \int_{\Gamma} dS \Big[ \boldsymbol{n}\cdot \kappa _1 \nabla \phi_1 + (\frac{\partial f_s}{\partial \phi_{1 \Gamma}} -\kappa_{1 \Gamma} \Delta_{\Gamma} \phi_1) \Big] \frac{\partial \phi_{1\Gamma}}{\partial t} + \int_{\Gamma} dS \Big[ \boldsymbol{n}\cdot \kappa_2 \nabla \phi_2 + (\frac{\partial f_s}{\partial \phi_{2 \Gamma}} -\kappa_{2 \Gamma} \Delta_{\Gamma} \phi_2) \Big] \frac{\partial \phi_{2\Gamma} }{\partial t} \\ & = -\int_{\Omega} dx \frac{1}{\nu} \Big[ M_{1}(\phi_1, \phi_2) |\nabla\mu_1|^2 + 2M_{12}(\phi_1, \phi_2)\nabla\mu_1\cdot\nabla\mu_2 + M_{2}(\phi_1, \phi_2) |\nabla\mu_2|^2 \Big] \\ & \quad - \int_{\Gamma} dS \Big[ \Gamma_1 \mu_{1\Gamma}^2 + 2\Gamma_{12}\mu_{1\Gamma}\mu_{2\Gamma} + \Gamma_2 \mu_{2\Gamma}^2 \Big]. \end{align}$

(2.24)

Since the mobility matrices ${\left[\begin{array}{cc} M_1 & M_{12}\\ M_{12} & M_2\\ \end{array}\right]}$ , ${\left[\begin{array}{cc} \Gamma_1 & \Gamma_{12}\\ \Gamma_{12} & \Gamma_2\\ \end{array}\right]}$ are positive definite, $\partial_t E \leq 0$ holds.

2.2. Non-dimensionalization

We set the characteristic length scale as the molecular length $l_0 = \nu^{1/3}$ and characteristic timescale $t_0 = \nu^{2/3}/D_1$ such that $\tilde{D}_1 = D_1 t_0/l_0^2 = 1$ . Rescaling $\tilde{\boldsymbol{x}} = \boldsymbol{x}/ l_0$ and $\tilde{t} = t /t_0$ , our model has the following non-dimensional parameters:

$\begin{align} \tilde{D}_{2} = \frac{D_2}{D_1} \, , \quad \tilde{D}_{12} = \frac{D_{12}}{D_1} \, , \quad \tilde{\Gamma_i} = \Gamma_i \frac{k_{\text{B}}T}{\nu_{s}} t_0 \, , \quad {\bar{ \kappa}_i} = {\kappa}_i \frac{\nu_{s}}{k_{\text{B}}T}\frac{1}{l_0} = \frac{\nu^{1/3}}{l_0} \, , \end{align}$

(2.25)

$\begin{align} \tilde{f_s} = f_s \frac{\nu_{s}}{k_{\text{B}}T} \, , \quad { \tilde{ \kappa }_{i\Gamma} = \frac{1}{l_0^{2/3}} \kappa_{i\Gamma} \frac{\nu_{s}}{k_{\text{B}}T} \, , \quad \tilde{ \kappa }_i = \frac{1}{l_0^{2/3}} \kappa_{i} \frac{\nu}{k_{\text{B}}T} }\, , i = 1, 2 . \end{align}$

(2.26)

For brevity, we skip the tildes in the following equations. The dimensionless equations governing the kinetics of the system are given as

$\begin{align} \partial_t \phi_1 & = \nabla \phi_1 (1-\phi_1-\phi_2) \cdot \nabla {\mu}_1 + \nabla D_{12} \phi_1 \phi_2 (1-\phi_1-\phi_2) \cdot \nabla {\mu}_2 , \end{align}$

(2.27a)

$\begin{align} \partial_t \phi_2 & = \nabla D_{12} \phi_1 \phi_2 (1-\phi_1-\phi_2) \cdot \nabla {\mu}_1 + \nabla D_2 \phi_2 (1-\phi_1-\phi_2) \cdot \nabla {\mu}_2, \end{align}$

(2.27b)

with the dimensionless boundary conditions written as

$\begin{align} \partial_t \phi_{1 \Gamma} & = -\Gamma_1 \Big[\boldsymbol{n} \cdot \bar{\kappa}_1 \nabla \phi_1 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}} \Big] - \Gamma_{12} \Big[\boldsymbol{n} \cdot \bar{\kappa}_2 \nabla \phi_2 - \kappa_{2 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}} \Big] , \end{align}$

(2.28a)

$\begin{align} \partial_t \phi_{2 \Gamma} & = -\Gamma_2 \Big[\boldsymbol{n} \cdot \bar{\kappa}_2 \nabla \phi_2 - \kappa_{2 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}}\Big] - \Gamma_{12} \Big[\boldsymbol{n} \cdot \bar{\kappa}_1 \nabla \phi_1 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}} \Big], \end{align}$

(2.28b)

$\begin{align} \partial_{\boldsymbol{n}} \mu_1|_{\Gamma} & = 0, \end{align}$

(2.28c)

$\begin{align} \partial_{\boldsymbol{n}} \mu_2|_{\Gamma} & = 0. \end{align}$

(2.28d)

In the next section, we will develop the numerical scheme based on the IEQ method ^[51,52,53].

3. Numerical scheme

We first reformulate the mathematical model into an equivalent form. Then, we discretize the reformulated system in time and space, respectively. In this work, we consider the 2D domain for simplicity.

3.1. Mathematical reformulation

We reformulate the total free energy in a quadratic form:

$\begin{align} q_1 = \sqrt{ \frac{\phi_1}{n_1} \ln \phi_1 + \frac{\phi_2}{n_2} \ln \phi_2 + (1-\phi_1 - \phi_2) \ln(1-\phi_1 - \phi_2) + \chi_{12}\phi_1\phi_2 + \chi_{13} \phi_1\phi_3 + \chi_{23} \phi_2\phi_3 + C}, \end{align}$

(3.1)

where $C$ is a constant such that the value under the square root is positive; $\phi_3 = 1-\phi_1-\phi_2$ . Then, the total energy becomes

$\begin{align} E = E_{\text{surf}} + E_{\text{bulk}} = E_{\text{surf}} + \int_{\Omega} \Big[ q_1^2 + \omega_1 \phi_1 + \omega_2 \phi_2 + \frac{\kappa_1}{2}|\nabla \phi_1|^2 + \frac{\kappa_2}{2}|\nabla \phi_2|^2 \Big]d\boldsymbol{x}; \end{align}$

(3.2)

the chemical potentials become

$\begin{align} \mu_1 & = 2 q_1 \frac{\partial q_1}{\partial \phi_1} + \omega_1 - \kappa_1 \Delta \phi_1 , \end{align}$

(3.3a)

$\begin{align} \mu_2 & = 2 q_1 \frac{\partial q_1}{\partial \phi_2} + \omega_2- \kappa_2 \Delta \phi_2. \end{align}$

(3.3b)

The reformulated governing equations of the system are as follows:

$\begin{align} \partial_t \phi_1 & = \nabla M_1(\phi_1, \phi_2) {\cdot} \nabla \mu_1 + \nabla M_{12}(\phi_1, \phi_2) {\cdot} \nabla \mu_2 , \end{align}$

(3.4a)

$\begin{align} \partial_t \phi_2 & = \nabla M_2(\phi_1, \phi_2) {\cdot} \nabla \mu_2 + \nabla M_{12}(\phi_1, \phi_2) {\cdot} \nabla \mu_1, \end{align}$

(3.4b)

$\begin{align} \partial_t q_1 & = \frac{\partial q_1}{\partial \phi_1} \partial_t \phi_1 + \frac{\partial q_1}{\partial \phi_2} \partial_t \phi_2, \end{align}$

(3.4c)

with the following boundary conditions:

$\begin{align} \partial_t \phi_{1 \Gamma} & = -\Gamma_1 \Big[\boldsymbol{n} \cdot \kappa_1 \nabla \phi_1 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}} \Big] - \Gamma_{12} \Big[\boldsymbol{n} \cdot \kappa_2 \nabla \phi_2 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}} \Big] , \end{align}$

(3.5a)

$\begin{align} \partial_t \phi_{2 \Gamma} & = -\Gamma_2 \Big[\boldsymbol{n} \cdot \kappa_2 \nabla \phi_2 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}}\Big] - \Gamma_{12} \Big[\boldsymbol{n} \cdot \kappa_1 \nabla \phi_1 - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}} \Big], \end{align}$

(3.5b)

$\begin{align} \partial_n \mu_1|_{\Gamma} & = 0, \end{align}$

(3.5c)

$\begin{align} \partial_n \mu_2|_{\Gamma} & = 0. \end{align}$

(3.5d)

The reformulated system has the following properties.

Theorem 3.1. The total mass of each component in the domain $\Omega$ is conserved, i.e.,

$\begin{align} \int_{\Omega} \phi_i(x, y, t) d\boldsymbol{x} = \int_{\Omega} \phi_i(x, y, 0) d\boldsymbol{x}, \quad i = 1, 2. \end{align}$

(3.6)

Theorem 3.2. The total free energy is dissipative, i.e.,

$\begin{align} \partial_t E = \partial_t (E_{\mathit{\text{surf}}} + E_{\mathit{\text{bulk}}}) \leq 0. \end{align}$

(3.7)

We skip the proof here. Next, we will discretize the reformulated governing equation in both the time and space directions.

3.2. Semi-discrete scheme in time

We use the Crank-Nicolson method in the time direction. $\Delta {t}$ is the time step. $(\cdot)^n$ represents the solution at the $n^{th}$ time step, i.e., $t_n = n\Delta t$ . We denote

$\begin{align} \partial_t^{n+1} (\cdot) = \frac{(\cdot)^{n+1}-(\cdot)^n}{\Delta t}, \quad (\cdot)^{n+1/2} = \frac{(\cdot)^{n} + (\cdot)^{n+1}}{2}, \quad \bar{(\cdot)}^{n+1/2} = \frac{3(\cdot)^{n} - (\cdot)^{n-1}}{2}. \end{align}$

(3.8)

The governing equations of the system are as follows:

$\begin{align} \partial_t^{n+1} \phi_1 & = \nabla M_1(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2}) {\cdot} \nabla \mu_1^{n+1/2} + \nabla M_{12}(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2}) {\cdot} \nabla \mu_2^{n+1/2} , \end{align}$

(3.9a)

$\begin{align} \partial_t^{n+1} \phi_2 & = \nabla M_2(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2}) {\cdot} \nabla \mu_2^{n+1/2} + \nabla M_{12}(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2}) {\cdot} \nabla \mu_1^{n+1/2}, \end{align}$

(3.9b)

$\begin{align} \partial_t^{n+1} q_1 & = \bar{\frac{\partial q_1}{\partial \phi_1}}^{n+1/2} \partial_t^{n+1} \phi_1 +\bar{\frac{\partial q_1}{\partial \phi_2}}^{n+1/2} \partial_t^{n+1} \phi_2, \end{align}$

(3.9c)

with the following boundary conditions:

$\begin{align} \partial_t^{n+1} \phi_{1 \Gamma} & = -\Gamma_1 \Big[\boldsymbol{n} \cdot \bar{\kappa}_1 \nabla \phi_1^{n+1/2} - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma}^{n+1/2} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}}^{n+1/2} \Big] \\ & \quad - \Gamma_{12} \Big[\boldsymbol{n} \cdot \bar{\kappa}_2 \nabla \phi_2^{n+1/2} - {\kappa}_{1 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma}^{n+1/2} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}}^{n+1/2} \Big] , \end{align}$

(3.10a)

$\begin{align} \partial_t^{n+1} \phi_{2 \Gamma} & = -\Gamma_2 \Big[\boldsymbol{n} \cdot \bar{\kappa}_2 \nabla \phi_2^{n+1/2} - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{2 \Gamma}^{n+1/2} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}}^{n+1/2}\Big] \\ & \quad - \Gamma_{12} \Big[\boldsymbol{n} \cdot \bar{\kappa}_1 \nabla \phi_1^{n+1/2} - \kappa_{1 \Gamma} \Delta_{\Gamma} \phi_{1 \Gamma}^{n+1/2} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}}^{n+1/2} \Big] , \end{align}$

(3.10b)

$\begin{align} \partial_{\boldsymbol{n}} \mu_1^{n+1}|_{\Gamma} & = 0, \end{align}$

(3.10c)

$\begin{align} \partial_{\boldsymbol{n}} \mu_2^{n+1}|_{\Gamma} & = 0, \end{align}$

(3.10d)

where

$\begin{align} \mu_1^{n+1/2} & = 2 q_1^{n+1/2} \bar{\frac{\partial q_1}{\partial \phi_1}}^{n+1/2} + \omega_1 - {\kappa}_1 \Delta \phi_1^{n+1/2} , \end{align}$

(3.11a)

$\begin{align} \mu_2^{n+1/2} & = 2 q_1^{n+1/2} \bar{\frac{\partial q_1}{\partial \phi_2}}^{n+1/2} + \omega_2 - {\kappa}_{2} \Delta \phi_2^{n+1/2}. \end{align}$

(3.11b)

3.3. Fully discrete scheme

We use the central finite-difference method such that there is second-order accuracy in the space. Specifically, we apply a uniform mesh in 2D space $[0, L] \times [0, L]$ (see ). In the $x$ direction, the domain is divided into $N_x$ equal-sized subintervals. Similarly, in the $y$ direction, the domain is divided into $N_y$ subintervals, each of equal size. We discretize the scalar functions given as $\phi_i$ and $\Delta \phi_i$ values at the center point $(x_i, y_j)$ , where

$\begin{align} x_i & = \left (i - \frac{1}{2} \right )\frac{L}{N_x}, \qquad i = 0, 1, \cdots, N_x+1, \end{align}$

(3.12)

$\begin{align} y_j & = \left (j - \frac{1}{2} \right )\frac{L}{N_y}, \qquad j = 0, 1, \cdots, N_y+1, \end{align}$

(3.13)

and we discretize the vector functions, e.g., $\nabla \phi_i$ , at the edge points $(x_{i+\frac{1}{2}}, y_{j})$ or $(x_{i}, y_{j+\frac{1}{2}})$ ^[55], where

$\begin{align} x_{i+\frac{1}{2}} & = i\frac{L}{N_x}, \qquad i = 0, 1, \cdots, N_x, \end{align}$

(3.14)

$\begin{align} y_{j+\frac{1}{2}} & = j\frac{L}{N_y}, \qquad j = 0, 1, \cdots, N_y. \end{align}$

(3.15)

At the boundary $\Gamma$ , we use the average value at adjacent discrete center points, e.g.,

$\begin{align} \phi_{1\Gamma}^n|_{\frac{1}{2}, j} = \frac{\phi_1^n|_{0, j}+\phi_1^n|_{1, j}}{2}, \qquad \phi_{1\Gamma}^n|_{N_x+\frac{1}{2}, j} = \frac{\phi_1^n|_{N_x, j}+\phi_1^n|_{N_x+1, j}}{2}, \end{align}$

(3.16)

$j = 1, \cdots N_y$ , is the index in the $y$ direction.

Figure 1. Sketch of the staggered mesh in 2D space. The blue line is the boundary

$\Gamma$ of the 2D domain

$\Omega = [0, L]\times[0, L]$ in the

$x-y$ plane. The black solid points are the center points inside the bulk, the black empty points are the ghost points adjacent to the boundary

$\Gamma$ and the crosses denote the edge points.

DownLoad: Full-Size Img PowerPoint

We respectively denote the edge-to-center average and difference operator as $a_x$ and $d_x$ in the $x$ direction and $a_y$ and $d_y$ in the $y$ direction, which are defined as follows:

$\begin{align} a_xu_{i, j}: = \frac{1}{2}(u_{i+\frac{1}{2}, j}+u_{i-\frac{1}{2}, j}), \quad &d_xu_{i, j}: = \frac{1}{h_x}(u_{i+\frac{1}{2}, j}-u_{i-\frac{1}{2}, j}), \end{align}$

(3.17)

$\begin{align} a_yv_{i, j}: = \frac{1}{2}(v_{i, j+\frac{1}{2}}+v_{i, j-\frac{1}{2}}), \quad &d_yu_{i, j}: = \frac{1}{h_y}(v_{i, j+\frac{1}{2}}-v_{i, j-\frac{1}{2}}). \end{align}$

(3.18)

We respectively denote the center-to-edge average and difference operators by $A_x$ and $D_x$ in the $x$ direction. Analogously, the center-to-edge average in the $y$ direction and difference operators are respectively denoted as $A_y$ and $D_y$ . Their definitions are as follows:

$\begin{align} A_x\phi_{i+\frac{1}{2}, j}: = \frac{1}{2}(\phi_{i+1, j}+\phi_{i, j}), \quad &D_x\phi_{i+\frac{1}{2}, j}: = \frac{1}{h_x}(\phi_{i+1, j}-\phi_{i, j}), \end{align}$

(3.19)

$\begin{align} A_y\phi_{i, j+\frac{1}{2}}: = \frac{1}{2}(\phi_{i, j+1}+\phi_{i, j}), \quad &D_y\phi_{i, j+\frac{1}{2}}: = \frac{1}{h_y}(\phi_{i, j+1}-\phi_{i, j}). \end{align}$

(3.20)

The standard 2D discrete Laplace operator is defined as ${\Delta}_h$ in the bulk and ${\Delta}_{\Gamma h}$ on the surface.

$\begin{equation} {\Delta}_h\phi : = d_x(D_x\phi)+d_y(D_y\phi), \quad {\Delta}_{\Gamma h }\phi_{\Gamma} : = d_y(D_y\phi_{\Gamma}), \quad \text{or}\quad {\Delta}_{\Gamma h }\phi_{\Gamma} : = d_x(D_x\phi_{\Gamma}). \end{equation}$

(3.21)

The discrete norm in the bulk and surface are defined as follows:

$\begin{align} ||\phi ||_c = h_x h_y \sum\limits_{i = 1}^{N_x} \sum\limits_{j = 1}^{N_y} \phi_{i, j}, \qquad ||\nabla \phi ||_v = h_x h_y \sum\limits_{i = 1}^{N_x} \sum\limits_{j = 1}^{N_y} \Big[a_xD_x \phi_{i+1/2, j} + a_yD_y \phi_{i, j+1/2}\Big], \end{align}$

(3.22)

$\begin{align} ||\phi ||_{c, \Gamma} = h_y \sum\limits_{j = 1}^{N_y} \phi_{j\Gamma} , \qquad ||\nabla \phi ||_{v, \Gamma} = h_y \sum\limits_{j = 1}^{N_y} \Big[ a_yD_y \phi_{i, j+1/2}\Big]. \end{align}$

(3.23)

The fully discrete numerical scheme:

$\begin{align} \Big \{\partial_t^{n+1} \phi_1 & = d_x A_x(M_1(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2}) D_x \mu_1^{n+1/2} + d_y A_y(M_1(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2})) D_y \mu_1^{n+1/2} \\ & + d_x A_x(M_{12}(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2}) D_x \mu_2^{n+1/2} + d_y A_y(M_{12}(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2})) D_y \mu_2^{n+1/2} \Big \}\Big|_{i, j}, \end{align}$

(3.24)

$\begin{align} \Big \{\partial_t^{n+1} \phi_2 & = d_x A_x(M_2(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2}) D_x \mu_2^{n+1/2} + d_y A_y(M_2(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2})) D_y \mu_2^{n+1/2} \\ & + d_x A_x(M_{12}(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2}) D_x \mu_1^{n+1/2} + d_y A_y(M_{12}(\bar{\phi_1}^{n+1/2}, \bar{\phi_2}^{n+1/2})) D_y \mu_1^{n+1/2} \Big \}\Big|_{i, j}, \end{align}$

(3.25)

$\begin{align} \Big \{\partial_t^{n+1} q_1 & = \bar{\frac{\partial q_1}{\partial \phi_1}}^{n+1/2} \partial_t^{n+1} \phi_1 +\bar{\frac{\partial q_1}{\partial \phi_2}}^{n+1/2} \partial_t^{n+1} \phi_2\Big \}\Big|_{i, j} \, , \end{align}$

(3.26)

where $i = 1, \cdots N_x$ , $j = 1, \cdots N_y$ , with the following boundary conditions:

$\begin{align} \Big \{ \partial_t^{n+1} \phi_{1 \Gamma} & = -\Gamma_1 \Big[\boldsymbol{n} \cdot \bar{\kappa}_1 d_x \phi_1^{n+1/2} - \kappa_{1 \Gamma} \Delta_{\Gamma h} \phi_{1 \Gamma}^{n+1/2} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}}^{n+1/2} \Big] \\ & \quad - \Gamma_{12} \Big[\boldsymbol{n} \cdot \bar{\kappa}_2 d_x \phi_2^{n+1/2} - {\kappa}_{2 \Gamma} \Delta_{\Gamma h} \phi_{2 \Gamma}^{n+1/2} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}}^{n+1/2} \Big] \Big \}\Big|_{\frac{1}{2}\text{or} (N_x+\frac{1}{2}) , j}, \end{align}$

(3.27a)

$\begin{align} \Big \{ \partial_t^{n+1} \phi_{2 \Gamma} & = -\Gamma_2 \Big[\boldsymbol{n} \cdot \bar{\kappa}_2 d_x \phi_2^{n+1/2} - \kappa_{2 \Gamma} \Delta_{\Gamma h} \phi_{2 \Gamma}^{n+1/2} + \frac{\partial f_s}{\partial \phi_{2 \Gamma}}^{n+1/2}\Big] \\ & \quad - \Gamma_{12} \Big[\boldsymbol{n} \cdot \bar{\kappa}_1 d_x \phi_1^{n+1/2} - \kappa_{1 \Gamma} \Delta_{\Gamma h} \phi_{1 \Gamma}^{n+1/2} + \frac{\partial f_s}{\partial \phi_{1 \Gamma}}^{n+1/2} \Big] \Big \}\Big|_{\frac{1}{2} \text{or} (N_x+\frac{1}{2}) , j}, \end{align}$

(3.27b)

$\begin{align} \mu_{k}|_{0, j} & = \mu_{k}|_{1, j}, \qquad \mu_{k}|_{N_x, j} = \mu_{k}|_{N_x+1, j}, \qquad k = 1, 2. \end{align}$

(3.27c)

$\begin{align} \mu_{k}|_{i, 0} & = \mu_{k}|_{i, 1}, \qquad \mu_{k}|_{i, N_y} = \mu_{k}|_{i, N_y+1}, \qquad k = 1, 2, \end{align}$

(3.27d)

where

$\begin{align} \mu_1^{n+1/2}|_{i, j} & = 2 q_1|_{i, j}^{n+1/2} \bar{\frac{\partial q_1}{\partial \phi_1}}^{n+1/2}|_{i, j} + \omega_1 - {\kappa}_1 \Delta_h \phi_1^{n+1/2}|_{i, j} , \end{align}$

(3.28a)

$\begin{align} \mu_2^{n+1/2}|_{i, j} & = 2 q_1|_{i, j}^{n+1/2} \bar{\frac{\partial q_1}{\partial \phi_2}}^{n+1/2}|_{i, j} + \omega_2 - {\kappa}_{2} \Delta_h \phi_2^{n+1/2}|_{i, j}. \end{align}$

(3.28b)

The fully discrete numerical scheme satisfies the conditions of the total mass conservation law and the total energy dissipation rate at the discrete level.

Theorem 3.3. Mass conservation law at discrete level: Based on the discrete boundary conditions given by (3.27),

$\begin{align} \sum\limits_{i = 1}^{N_x}\sum\limits_{j = 1}^{N_y} \frac{\phi_1^{n+1}|_{i, j} - \phi_1^n|_{i, j}}{\Delta t} = \sum\limits_{i = 1}^{N_x}\sum\limits_{j = 1}^{N_y} \Big [ \nabla \cdot M_1 \nabla \mu_1^{n+1/2}|_{i, j} + \nabla \cdot M_{12} \nabla \mu_2^{n+1/2}|_{i, j} \Big] = 0; \end{align}$

(3.29)

it is the same as for the total mass of $\phi_2$ .

Theorem 3.4. Unconditional energy stability at discrete level: The discrete energy of the system is defined as

$\begin{align} E^{n+1} = E_{\mathit{\text{bulk}}}^{n+1} + E_{\mathit{\text{surf}}}^{n+1}, \end{align}$

(3.30)

where

$\begin{align} E_{\mathit{\text{bulk}}}^{n+1} & = || f_b(\phi_1, \phi_2)^{n+1}||_c + \frac{\kappa_1}{2}||\nabla \phi_1^{n+1}||_v^2 + \frac{\kappa_2}{2}||\nabla \phi_2^{n+1}||_v^2 , \end{align}$

(3.31)

$\begin{align} E_{\mathit{\text{surf}}}^{n+1} & = ||f_s(\phi_{1 \Gamma}, \phi_{2 \Gamma})^{n+1}||_{c, \Gamma} + \frac{\kappa_{1\Gamma}}{2}||\nabla \phi_{1 \Gamma}^{n+1}||_{v, \Gamma}^2 + \frac{\kappa_{2\Gamma}}{2}||\nabla \phi_{2 \Gamma}^{n+1}||_{v, \Gamma}^2. \end{align}$

(3.32)

The total energy at the discrete level decreases with respect to time, i.e.,

$\begin{align} E^{n+1} \leq E^n. \end{align}$

(3.33)

proof: Using the definition of reformulated total energy at the discrete level given by (3.30) and the formula $(a^2 - b^2) = (a + b) (a -b)$ , we obtain

$\begin{align} \frac{E^{n+1}-E^n}{\Delta t} = &|| 2\cdot \frac{q_1^{n+1}-q_1^n}{\Delta t} \, \frac{q_1^{n+1} + q_1^n}{2} + \omega_1 \frac{\phi_1^{n+1}-\phi_1^n}{\Delta t} + \omega_2 \frac{\phi_2^{n+1}-\phi_2^n}{\Delta t}||_c \\ & + \frac{\kappa_1}{2} ||2\cdot \frac{\nabla \phi_1^{n+1}-\nabla \phi_1^n}{\Delta t} \, \frac{\nabla \phi_1^{n+1} + \nabla \phi_1^n}{2}||_v + \frac{\kappa_2}{2} ||2\cdot \frac{\nabla \phi_2^{n+1}-\nabla \phi_2^n}{\Delta t} \, \frac{\nabla \phi_2^{n+1} + \nabla \phi_2^n}{2}||_v \\ & + \frac{\kappa_{1\Gamma}}{2} ||2\cdot \frac{\nabla_{||} \phi_{1\Gamma}^{n+1}-\nabla_{||} \phi_{1\Gamma}^n}{\Delta t} \, \frac{\nabla_{||} \phi_{1\Gamma}^{n+1} + \nabla_{||} \phi_{1\Gamma}^n}{2}||_{v, \Gamma} \\ & + \frac{\kappa_{2\Gamma}}{2} ||2\cdot \frac{\nabla_{||} \phi_{2\Gamma}^{n+1}-\nabla_{||} \phi_{2\Gamma}^n}{\Delta t} \, \frac{\nabla_{||} \phi_{2\Gamma}^{n+1} + \nabla_{||} \phi_{2\Gamma}^n}{2}||_{v, \Gamma} \\ & + h_1 ||\frac{\phi_1^{n+1}-\phi_1^n}{\Delta t} ||_c + g_1 ||2\cdot \frac{\phi_1^{n+1}-\phi_1^n}{\Delta t} \frac{\phi_1^{n+1} + \phi_1^n}{2}||_c \\ & + h_2||\frac{\phi_2^{n+1}-\phi_2^n}{\Delta t}||_c + g_2 ||\frac{\phi_2^{n+1}-\phi_2^n}{\Delta t} \frac{\phi_2^{n+1} + \phi_2^n}{2} ||_c\\ & + \gamma ||\frac{\phi_1^{n+1} + \phi_1^n}{2} \, \frac{\phi_2^{n+1} - \phi_2^n}{\Delta t}||_c + \gamma ||\frac{\phi_2^{n+1} + \phi_2^n}{2} \, \frac{\phi_1^{n+1} - \phi_1^n}{\Delta t}||_c. \end{align}$

(3.34)

According to the discrete chemical potentials defined in (3.28) and discrete boundary conditions given by (3.27), we have

$\begin{align} \frac{E^{n+1}-E^n}{\Delta t} & = ||\mu_1^{n+1/2}\cdot(\nabla M_1 {\cdot} \nabla \mu_1^{n+1/2} + \nabla M_{12} {\cdot} \nabla \mu_2^{n+1/2}) ||_{v} \\ & \quad + ||\mu_2^{n+1/2}\cdot(\nabla M_2 {\cdot} \nabla \mu_2^{n+1/2} + \nabla M_{12} {\cdot} \nabla \mu_1^{n+1/2}) ||_{v} \\ & \quad + ||\mu_{1\Gamma}^{n+1/2}\cdot(-\Gamma_{1} \mu_{1\Gamma}^{n+1/2} - \Gamma_{12} \mu_{2\Gamma}^{n+1/2})||_{v, \Gamma} \\ & \quad + ||\mu_{2\Gamma}^{n+1/2}\cdot(-\Gamma_{2} \mu_{2\Gamma}^{n+1/2} - \Gamma_{12} \mu_{1\Gamma}^{n+1/2})||_{v, \Gamma} \\ & = - || M_1 |\nabla \mu_1^{n+1/2}|^2 + 2M_{12}\nabla \mu_1^{n+1/2}\cdot \nabla \mu_2^{n+1/2} + M_2|\nabla \mu_2^{n+1/2}|^2 ||_{v} \\ & \quad - || \Gamma_{1} |\mu_{1\Gamma}^{n+1/2}|^2 + 2\Gamma_{12}\mu_{1\Gamma}^{n+1/2}\mu_{2\Gamma}^{n+1/2} + \Gamma_{2} |\mu_{2\Gamma}^{n+1/2}|^2 ||_{v, \Gamma}^2 \leq 0, \end{align}$

(3.35)

i.e., (3.33) is obtained.

4. Physical phenomenon investigation by numerical simulations

Taking advantage of the ability to vary parameters in a controlled manner and over a wide range of scales in numerical simulations, we investigate several physical phenomena in this section, including the wettability of multi-component droplets on solid surfaces, the effects of wall-mixture interactions on patterns in the bulk and the role of cross-coupling relaxation rates in controlling kinetic processes in both the bulk and surface.

According to the Gibbs rule, a ternary mixture system can have at most three coexisting phases in a system (see Appendix for details). In this section, we will evidence our points in both two-phase and three-phase coexisting scenarios.

4.1. Contact angle is determined by the additive effects of wall-mixture interaction

The contact angle is a measure of the wetting behavior of a liquid on a solid surface. It is determined by the balance of forces between the solid phase and two other liquid phases on the top of the solid surface. The wall-mixture interaction plays a significant role in determining the contact angle.

In the case of a multi-component liquid mixture, the contact angle is determined by the combined effects of the interactions between the individual component of the mixture and the solid surface. To grasp these additive effects, we define the surface free energy as the summation of interactions between individual components and the solid wall.

Specifically, the surface free energy is given by

$\begin{equation} f_s(\phi_{1\Gamma}, \phi_{2\Gamma}) = \frac{k_B T}{\nu_s}\Big[ h_1 \phi_{1\Gamma} + g_1 \phi_{1\Gamma}^2 + h_2 \phi_{2\Gamma} + g_2 \phi_{2\Gamma}^2 + \gamma \phi_{1\Gamma} \phi_{2\Gamma} \Big]. \end{equation}$

(4.1)

Here, $h_i \phi_i$ depicts the interaction between the solid wall and $i$ -th component and $g_i \phi_i^2$ represents the interaction among $i$ -th component molecules near the solid surface. $\gamma \phi_1 \phi_2$ denotes the interaction between different kinds of molecules.

If all of the interactions mentioned above are mutual, i.e., $h_i = g_i = \gamma = 0$ , the contact angle of droplets on the surface is $90^{\text{o}}$ (see ). If the interactions between each component and the solid wall are all attractive, i.e., $h_i < 0$ , $g_i < 0$ or $\gamma < 0$ , the contact angle will decrease until complete wetting is achieved (see Figures 2–4, 5g and , ), while repulsive interactions will lead to a larger contact angle, i.e., $\theta > 90^{\text{o}}$ (see Figures 2– and ). If interactions between individual components and the solid wall are different, e.g., the interaction between $\phi_1$ and the solid wall is attractive, while the interaction between $\phi_2$ and the solid wall is repulsive, the effective interaction between the droplet and solid wall is determined by the dominant component with a higher concentration in the droplet (see Figures 2–4 and 5a, i). These observations can be verified in the coexistence of the three phases (see Appendix) wetting phenomena (see Figures 7 and 8).

Figure 2. Profiles of

$\phi_1$ with a linear interaction between solid walls and mixture components along the contact angle of a two-phase coexisting mixture in equilibrium. The snapshots denote the profiles of component 1, i.e.,

$\phi_1$ , with different values of

$h_1$ and

$h_2$ . At each snapshot, the yellow area depicts the condensates and a white line is used to denote the tangent line of the spherical interface at the contact point on the surface

$\Gamma$ . The contact angle of a droplet between the tangent line and surface is denoted by

$\theta$ . In each column (row), we use the same

$h_1$ (

$h_2$ ) values. The basic parameter values used in this case are as follows:

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 0$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$g_1 = g_2 = \gamma = 0$ ,

$\chi_{12} = -1, \chi_{23} = 0, \chi_{13} = 2.5$ . Since the droplet is a composite of

$\phi_1$ , mainly, and the contact angle is determined by the value of

$h_1$ in this study.

DownLoad: Full-Size Img PowerPoint

Figure 3. Profiles of

$\phi_2$ with a linear interaction between solid walls and mixture components along the contact angle of a two-phase coexisting mixture in equilibrium. The snapshots denote the profiles of component 2, i.e.,

$\phi_2$ , with different values of

$h_1$ and

$h_2$ . The profile of

$\phi_2$ changes rapidly near the surface with respect to different

$h_2$ values. However, the contact angle is determined by the value of

$h_1$ since the the droplet mainly contains

$\phi_1$ . The basic parameter values used in this case are as follows:

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 0$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$g_1 = g_2 = \gamma = 0$ ,

$\chi_{12} = -1, \chi_{23} = 0, \chi_{13} = 2.5$ .

DownLoad: Full-Size Img PowerPoint

Figure 4. Profiles of

$\phi_1$ with a quadratic interaction between solid walls and mixture components along the contact angle of a two-phase coexisting mixture in equilibrium. The snapshots denote the profiles of component 1, i.e.,

$\phi_1$ , with different values of

$g_1$ and

$g_2$ . We observe that the interaction between

$\phi_2$ and a solid wall mainly influences the contact angle. The basic parameter values used in this case are as follows:

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 0$ ,

$\kappa_1 = \kappa_2 = 0$ ,

$h_1 = h_2 = \gamma = 0$ ,

$\chi_{12} = -1, \chi_{23} = 0, \chi_{13} = 2.5$ .

DownLoad: Full-Size Img PowerPoint

Figure 5. Profiles of

$\phi_2$ with a quadratic interaction between solid walls and mixture components along the contact angle of a two-phase coexisting mixture in equilibrium. The snapshots denote the profiles of component 2, i.e.,

$\phi_2$ , with different values of

$g_1$ and

$g_2$ . We find that the interaction between

$\phi_1$ and a solid surface changes the profiles of

$\phi_2$ near the surface, while the contact angle is determined by the interaction between

$\phi_1$ and a solid wall. The basic parameter values used in this case are as follows:

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 0$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$h_1 = h_2 = \gamma = 0$ ,

$\chi_{12} = -1, \chi_{23} = 0, \chi_{13} = 2.5$ .

DownLoad: Full-Size Img PowerPoint

Figure 6. Profiles of

$\phi_1$ and

$\phi_2$ with a nonlinear coupling interaction between

$\phi_1$ and

$\phi_2$ in equilibrium. Plots (

$\textbf{a}, \textbf{b}, \textbf{c}$ ) represent the concentration profiles of

$\phi_1$ , while plots (

$\textbf{d}, \textbf{e}, \textbf{f}$ ) represent the concentration profiles of

$\phi_2$ . We observe that, as the negative coupling interaction between

$\phi_1$ and

$\phi_2$ becomes stronger, the contact angle decreases until the condensate (yellow) completely wets on the solid surface (plots (

$\textbf{a}, \textbf{d}$ ). Here, we only investigate the negative coupling interaction since the positive coupling interaction will lead to nonphysical solutions. The basic parameter values used in this case are as follows:

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 0$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$h_1 = h_2 = g_1 = g_2 = 0$ ,

$\chi_{12} = -1, \chi_{23} = 0, \chi_{13} = 2.5$ .

DownLoad: Full-Size Img PowerPoint

Figure 7. Profiles of

$\phi_1$ with a quadratic interaction between the solid wall and mixture components along the contact angle of a three-phase coexisting mixture at equilibrium. With the interaction parameter values

$\chi_{12} = \chi_{23} = \chi_{13} = 3$ , the ternary mixture has three coexisting phases (respectively denoted by red numbers 1, 2, 3 in plot e). We changed the strength of the quadratic interaction, i.e., the values of

$g_i, i = 1, 2$ . We observe that the contact angle of droplets 1 and 2 are mainly determined by the interaction between the dominant component and the solid wall, i.e.,

$\phi_1$ in droplet 1 and

$\phi_2$ in droplet 2. Other basic parameter values used in this case are as follows:

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 0$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$h_1 = h_2 = \gamma = 0$ .

DownLoad: Full-Size Img PowerPoint

Figure 8. Profiles of

$\phi_2$ with a quadratic interaction between the solid wall and mixture components along the contact angle of a three-phase coexisting mixture at equilibrium. The basic parameter values used in this case are as follows:

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 0$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$h_1 = h_2 = \gamma = 0$ ,

$\chi_{12} = \chi_{23} = \chi_{13} = 3$ .

DownLoad: Full-Size Img PowerPoint

4.2. Wall-mixture interaction changes the condensates patterns far away from the surface

The influence of wall-mixture interactions in the formation of condensate patterns in systems exhibiting three-phase coexistence has yet to be extensively investigated. This study aims to bridge this gap by delving into the phase separation and condensate coarsening processes within such systems, paying special attention to the varying types of wall-mixture interactions.

Our numerical studies revealed that interactions between the wall and the mixture exert a profound effect on the condensate patterns in the bulk, even in regions distant from the surface. We observed distinct variations in the condensate configurations that are contingent on the nature of the wall-mixture interactions. For instance, in cases in which strong attractive or repulsive forces were at play, the resultant patterns manifested as strips running parallel to the solid wall, as evidenced by Figures 10 and 11a–d. Conversely, when the interactions were weak, the patterns that emerged were less structured and more sporadic, as can be seen in Figures 10 and 11e. These findings underscore the significance of wall-mixture interactions in determining the arrangement of condensates within the bulk.

Figure 9. Relaxation of dynamic boundary conditions controls the droplet spreading on the surface. Plots (a–i) are the snapshots of component 1, i.e.,

$\phi_1$ . We do not consider the cross-coupling interaction in this simulation, i.e.,

$\Gamma_{12} = 0$ . The relaxation rate

$\Gamma_i$ ,

$i = 1, 2$ decreases from 1 to

$10^{-4}$ , and the spreading process slows down correspondingly. Other basic parameter values used in this case are as follows:

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 0$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$h_1 = 0, h_2 = 0, g_1 = -1, g_2 = 0, \gamma = 0$ ,

$\chi_{12} = -1$ ,

$\chi_{23} = 0$ ,

$\chi_{13} = 2.5$ .

DownLoad: Full-Size Img PowerPoint

Figure 10. Profiles of

$\phi_1$ during the spontaneous phase separation of the ternary mixture with three-phase coexistence with various types of wall-mixture interactions. We show the snapshots of

$\phi_1$ at

$T = 10,000$ . Under the condition of neutral wall-mixture interaction (plot e, i.e.,

$g_i = 0, i = 1, 2$ ), the condensates (yellow) distribute in the bulk homogeneously. However, under the condition of attractive/repulsive wall-mixture interaction, the condensates distribute horizontally, parallel with the solid wall (see plots a, b, c, d, f, g, i). The basic parameter values used in this case are as follows:

$h_1 = h_2 = 0$ ,

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 1$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$\gamma = 0$ ,

$\chi_{12} = \chi_{23} = \chi_{13} = 3$ .

DownLoad: Full-Size Img PowerPoint

Figure 11. Profiles of

$\phi_2$ during the spontaneous phase separation of the ternary mixture with three-phase coexistence with various types of wall-mixture interactions. The basic parameter values used in this case are as follows:

$h_1 = h_2 = 0$ ,

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 1$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$\gamma = 0$ ,

$\chi_{12} = \chi_{23} = \chi_{13} = 3$ .

DownLoad: Full-Size Img PowerPoint

4.3. Relaxation rate controls the kinetic process in the bulk and surface

The relaxation rate serves as a pivotal factor that influences the kinetic processes occurring within the bulk and on the surface of a system. To elucidate the impact of kinetic rates in the spreading of droplets on a solid surface, we placed a droplet atop the surface and varied the kinetic rates, denoted as $\Gamma_i$ , from $10^{-4}$ to 1. Our observations revealed that the surface relaxation rates play a decisive role in dictating the kinetics of droplet spreading. Specifically, a faster surface relaxation rate accelerates the spreading of the droplet, as depicted in Figure 9.

The exploration of the effects of cross-coupling coefficients has been limited in previous research. In this study, we altered the cross-coupling coefficient, $m_{12}$ , in the bulk to investigate its impact on the phase separation process therein. Our simulations indicate that changes in cross-coupling can indeed alter the kinetic processes (see and for the snapshots). The total masses of $\phi_1$ and $\phi_2$ (see ) with different cross-coupling coefficient $m_{12}$ values are constants, as proved by Theorem 3.3. From the total energy profiles (see ), it is evident that variations in $m_{12}$ influence the total energy evolution during the initial stages of spontaneous phase separation. Further, we adjusted the cross-coupling coefficient on the surface, $\Gamma_{12}$ , to scrutinize its effects on the kinetics of droplet spreading on the solid surface. Our findings suggest that cross-coupling exerts a minimal effect on droplet spreading and coarsening subsequent to the phase separation process.

Figure 12. Profiles of

$\phi_1$ during spontaneous phase separation evolution at different time points with different values of the cross-coupling coefficient

$m_{12}$ . We find that changing the cross-coupling coefficient

$m_{12}$ leads to changes in the phase separation and following coarsening processes. The basic parameter values used in this study are as follows:

$h_1 = h_1 = g_1 = g_2 = \gamma = 0$ ,

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 1$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$\gamma = 0$ ,

$\chi_{12} = \chi_{23} = \chi_{13} = 3$ .

DownLoad: Full-Size Img PowerPoint

Figure 13. Profiles of

$\phi_2$ during spontaneous phase separation evolution at different time points with different values of cross-coupling coefficient

$m_{12}$ . The basic parameter values used in this case are as follows:

$h_1 = h_1 = g_1 = g_2 = \gamma = 0$ ,

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 1$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$\gamma = 0$ ,

$\chi_{12} = \chi_{23} = \chi_{13} = 3$ .

DownLoad: Full-Size Img PowerPoint

Figure 14. Total masses of

$\phi_1$ and

$\phi_2$ in the system during the spontaneous phase separation process. The basic parameter values used in this case are as follows:

$h_1 = h_1 = g_1 = g_2 = \gamma = 0$ ,

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 1$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$\gamma = 0$ ,

$\chi_{12} = \chi_{23} = \chi_{13} = 3$ .

DownLoad: Full-Size Img PowerPoint

Figure 15. Total energy profiles during spontaneous phase separation process with different values of the cross-coupling coefficient

$m_{12}$ . The basic parameter values used in this case are as follows:

$h_1 = h_1 = g_1 = g_2 = \gamma = 0$ ,

$\kappa_{1\Gamma} = \kappa_{2\Gamma} = 1$ ,

$\kappa_1 = \kappa_2 = 1$ ,

$\gamma = 0$ ,

$\chi_{12} = \chi_{23} = \chi_{13} = 3$ .

DownLoad: Full-Size Img PowerPoint

5. Conclusions

In this paper, we have presented a systematic derivation of multi-component, multi-phase mixtures with dynamic boundary conditions. The governing equations in the models are composed of the mass conservation and constitutive equations, which are derived by using the Onsager principle to preserve the energy dissipation rate in time. The gradient flow structure is $H^{-1}(\Omega)$ in the bulk and $L_2(\Gamma)$ on the surface.

Moreover, we have presented a second-order, fully discrete, linear and unconditionally energy-stable numerical scheme for the ternary mixture model with dynamic boundary conditions. First, we reformulate the model by introducing intermediate variables by implementing the IEQ strategy. Using the reformulated model equations, we develop a second-order, energy-stable, semi-discrete numerical scheme in time. Then, we obtain a fully discrete numerical scheme by applying the finite-difference method to the staggered grid in space, which preserves a fully discrete energy dissipation rate and the total mass conservation laws.

Beyond theoretical advancements, our work has substantial engineering implications. Utilizing our efficient and accurate numerical solver, we have investigated the effects of short-ranged interaction between the mixture components and solid surface on the wettability of the surface, from the perspectives of both stationary and kinetic solutions. We found that the contact angle of condensates are determined by the additive effects introduced by the interaction between each component and the solid wall. The droplet spreading can be affected hugely by the surface relaxation rates, while the cross-coupling relaxation rate for the surface does not influence the droplet spreading a lot. Furthermore, the spontaneous phase separation phenomena of ternary mixtures with two phases and three phases, and subject to dynamic boundary conditions, have been explored. Our results suggest that strong attractive/repulsive coupling between the surface and bulk leads to ordered condensate patterns in the bulk. Moreover, the cross-coupling in the bulk changes the phase separation process significantly. This insight is invaluable in fields like semiconductor manufacturing, where surface properties are critical, or in medical device fabrication, where understanding biofluid interactions with surfaces can lead to better product designs.

Furthermore, our model and the scheme can be readily extended to models of $N$ -component mixtures with $N > 3$ . Our work provides a general framework for the study of multi-component mixtures with dynamic boundary conditions. The multi-component mixtures with more complex boundary conditions, e.g., those developed as Goldstein (1.6), Liu-Wu model (1.7) and KLLM model (1.8), can be studied based on this framework. This work not only provides a foundational framework for future study of multi-component mixtures with dynamic boundary conditions for a wide range of applications, it also exerts a significant influence in the study of complex fluid dynamics.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

We thank C. A. Weber and Z. Zhang for insightful discussions. X. Zhao acknowledges the "FoSE New Researchers Grant" of the University of Nottingham Ningbo China for financial support. We also thank the anonymous referees for their critical comments.

Conflict of interest

The authors declare that there is no conflict of interest.

References

[1]	C. Brangwynne, C. R. Eckmann, D. S. Courson, A. Rybarska, C. Hoege, J. Gharakhani, et al., Germline p granules are liquid droplets that localize by controlled dissolution/condensation, Science, 324 (2009), 1729–1732. https://doi.org/10.1126/science.1172046 doi: 10.1126/science.1172046
[2]	J. G. Gall, M. Bellini, Z. Wu, C. Murphy, Assembly of the nuclear transcription and processing machinery: cajal bodies (coiled bodies) and transcriptosomes, Mol. Biol. Cell, 10 (1999), 4385–4402. https://doi.org/10.1091/mbc.10.12.4385 doi: 10.1091/mbc.10.12.4385
[3]	J. Agudo-Canalejo, S. W. Schultz, H. Chino, S. M. Migliano, C. Saito, I. Koyama-Honda, et al., Wetting regulates autophagy of phase-separated compartments and the cytosol, Nature, 591 (2021), 142–146. https://doi.org/10.1038/s41586-020-2992-3 doi: 10.1038/s41586-020-2992-3
[4]	P. Tabeling, Introduction to Microfluidics, Oxford University Press, 2005.
[5]	F. Bruder, R. Brenn, Spinodal decomposition in thin films of a polymer blend, Phys. Rev. Lett., 69 (1992), 624–627. https://doi.org/10.1103/PhysRevLett.69.624 doi: 10.1103/PhysRevLett.69.624
[6]	G. Krausch, C. Dai, E. J. Kramer, F. S. Bates, Real space observation of dynamic scaling in a critical polymer mixture, Phys. Rev. Lett., 71 (1993), 3669–3672. https://doi.org/10.1103/PhysRevLett.71.3669 doi: 10.1103/PhysRevLett.71.3669
[7]	L. Sung, A. Karim, J. F. Douglas, C. C. Han, Dimensional crossover in the phase separation kinetics of thin polymer blend films, Phys. Rev. Lett., 76 (1996), 4368–4371. https://doi.org/10.1103/PhysRevLett.76.4368 doi: 10.1103/PhysRevLett.76.4368
[8]	A. K. Gunstensen, D. H. Rothman, S. Zaleski, G. Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A, 43 (1991), 4320–4327. https://doi.org/10.1103/PhysRevA.43.4320 doi: 10.1103/PhysRevA.43.4320
[9]	X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E, 47 (1993), 1815–1819. https://doi.org/10.1103/PhysRevE.47.1815 doi: 10.1103/PhysRevE.47.1815
[10]	X. Shan, H. Chen, Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E, 49 (1994), 2941–2948. https://doi.org/10.1103/PhysRevE.49.2941 doi: 10.1103/PhysRevE.49.2941
[11]	M. R. Swift, W. R. Osborn, J. M. Yeomans, Lattice Boltzmann simulation of nonideal fluids, Phys. Rev. Lett., 75 (1995), 830–833. https://doi.org/10.1103/PhysRevLett.75.830 doi: 10.1103/PhysRevLett.75.830
[12]	M. R. Swift, E. Orlandini, W. R. Osborn, J. M. Yeomans, Lattice Boltzmann simulations of liquid-gas and binary fluid systems, Phys. Rev. E, 54 (1996), 5041–5052. https://doi.org/10.1103/PhysRevE.54.5041 doi: 10.1103/PhysRevE.54.5041
[13]	S. Schmieschek, J. Harting, Contact angle determination in multicomponent lattice boltzmann simulations, Commun. Comput. Phys., 9 (2011), 1165–1178. https://doi.org/10.4208/cicp.201009.271010s doi: 10.4208/cicp.201009.271010s
[14]	G. Yan, Z. Li, T. Bore, S. A. G. Torres, A. Scheuermann, L. Li, Discovery of dynamic two-phase flow in porous media using two-dimensional multiphase lattice Boltzmann simulation, Energies, 14 (2021), 4044. https://doi.org/10.3390/en14134044 doi: 10.3390/en14134044
[15]	G. Yan, Z. Li, T. Bore, S. A. G. Torres, A. Scheuermann, L. Li, A lattice Boltzmann exploration of two-phase displacement in 2D porous media under various pressure boundary conditions, J. Rock Mech. Geotech. Eng., 14 (2022), 1782–1798. https://doi.org/10.1016/j.jrmge.2022.05.003 doi: 10.1016/j.jrmge.2022.05.003
[16]	H. Liang, J. Xu, J. Chen, Z. Chai, B. Shi, Lattice Boltzmann modeling of wall-bounded ternary fluid flows, Appl. Math. Modell., 73 (2019), 487–513. https://doi.org/10.1016/j.apm.2019.03.009 doi: 10.1016/j.apm.2019.03.009
[17]	A. Ferrari, I, Lunati, Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy, Adv. Water Resour., 57 (2013), 19–31. https://doi.org/10.1016/j.advwatres.2013.03.005 doi: 10.1016/j.advwatres.2013.03.005
[18]	A. Ferrari, I, Lunati, Inertial effects during irreversible meniscus reconfiguration in angular pores, Adv. Water Resour., 74 (2014), 1–13. https://doi.org/10.1016/j.advwatres.2014.07.009 doi: 10.1016/j.advwatres.2014.07.009
[19]	A. Ferrari, J. Jimenez-Martinez, T. L. Borgne, Y. Méheust, I. Lunati, Challenges in modeling unstable two-phase flow experiments in porous micromodels, Water Resour. Res., 51 (2015), 1381–1400. https://doi.org/10.1002/2014WR016384 doi: 10.1002/2014WR016384
[20]	Z. Li, S. Galindo-Torres, A. Scheuermann, L. Li, Mesoscopic approach to fluid-solid interaction: Apparent liquid slippage and its effect on permeability estimation, Phys. Rev. E, 98 (2018), 052803. https://doi.org/10.1103/PhysRevE.98.052803 doi: 10.1103/PhysRevE.98.052803
[21]	Z. Li, J. Li, G. Yan, S. Galindo-Torres, A. Scheuermann, L. Li, Mesoscopic model framework for liquid slip in a confined parallel-plate flow channel, Phys. Rev. Fluids, 6 (2021), 034203. https://doi.org/10.1103/PhysRevFluids.6.034203 doi: 10.1103/PhysRevFluids.6.034203
[22]	J. W. Cahn, J. E. Hilliard, Free energy of a nonuniform system. Ⅰ. interfacial free energy, J. Chem. Phys., 28 (1958), 258–267. https://doi.org/10.1063/1.1744102 doi: 10.1063/1.1744102
[23]	J. W. Cahn, On spinodal decomposition in cubic crystals, Acta Metall., 10 (1962), 179–183. https://doi.org/10.1016/0001-6160(62)90114-1 doi: 10.1016/0001-6160(62)90114-1
[24]	F. Eliot, E. G. Morton, Continuum theory of thermally induced phase transitions based on an order parameter, Physica D, 68 (1993), 326–343. https://doi.org/10.1016/0167-2789(93)90128-N doi: 10.1016/0167-2789(93)90128-N
[25]	F. Bai, C. M. Elliott, A. Gardiner, A. Spence, A. M. Stuart, The viscous Cahn–Hilliard equation. Ⅰ. Computations, Nonlinearity, 8 (1995), 131–160. https://doi.org/10.1088/0951-7715/8/2/002 doi: 10.1088/0951-7715/8/2/002
[26]	M. E. Gurtin, Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance, Physica D, 92 (1996), 178–192. https://doi.org/10.1016/0167-2789(95)00173-5 doi: 10.1016/0167-2789(95)00173-5
[27]	H. P. Fischer, P. Maass, W. Dieterich, Novel surface modes in spinodal decomposition, Phys. Rev. Lett., 79 (1997), 893–896. https://doi.org/10.1103/PhysRevLett.79.893 doi: 10.1103/PhysRevLett.79.893
[28]	R. Kenzler, F. Eurich, P. Maass, B. Rinn, J. Schropp, E. Bohl, et al., Phase separation in confined geometries: Solving the Cahn–Hilliard equation with generic boundary conditions, Comput. Phys. Commun., 133 (2001), 139–157. https://doi.org/10.1016/S0010-4655(00)00159-4 doi: 10.1016/S0010-4655(00)00159-4
[29]	R. G. Gisèle, M. Alain, S. Giulio, A Cahn–Hilliard model in a domain with non-permeable walls, Physica D, 240 (2011), 754–766. https://doi.org/10.1016/j.physd.2010.12.007 doi: 10.1016/j.physd.2010.12.007
[30]	C. Liu, H. Wu, A Cahn–Hilliard model in a domain with non-permeable walls, Arch. Ration. Mech. Anal., 233 (2019), 167–247. https://doi.org/10.1007/s00205-019-01356-x doi: 10.1007/s00205-019-01356-x
[31]	P. Knopf, K. F. Lam, C. Liu, S. Metzger, Phase-field dynamics with transfer of materials: The Cahn–Hilliard equation with reaction rate dependent dynamic boundary conditions, ESAIM Math. Model. Numer. Anal., 55 (2021), 229–282. https://doi.org/10.1051/m2an/2020090 doi: 10.1051/m2an/2020090
[32]	K. Binder, H. L. Frisch, Dynamics of surface enrichment: A theory based on the Kawasaki spin-exchange model in the presence of a wall, Z. Phys. B: Condens. Matter, 84 (1991), 403–418. https://doi.org/10.1007/BF01314015 doi: 10.1007/BF01314015
[33]	S. Puri, K. Binder, Surface effects on spinodal decomposition in binary mixtures and the interplay with wetting phenomena, Phys. Rev. E, 49 (1994), 5359–5377. https://doi.org/10.1103/PhysRevE.49.5359 doi: 10.1103/PhysRevE.49.5359
[34]	H. Garcke, P. Knopf, S. Yayla, Long-time dynamics of the Cahn–Hilliard equation with kinetic rate dependent dynamic boundary conditions, Nonlinear Anal., 215 (2022), 112619. https://doi.org/10.1016/j.na.2021.112619 doi: 10.1016/j.na.2021.112619
[35]	H. Wu, A review on the Cahn–Hilliard equation: classical results and recent advances in dynamic boundary conditions, Electron. Res. Arch., 30 (2022), 2788–2832. https://doi.org/10.3934/era.2022143 doi: 10.3934/era.2022143
[36]	T. Fukao, S. Yoshikawa, S. Wada, Structure-preserving finite difference schemes for the Cahn-Hilliard equation with dynamic boundary conditions in the one-dimensional case, Commun. Pure Appl. Anal., 16 (2017), 1915–1938. https://doi.org/10.3934/cpaa.2017093 doi: 10.3934/cpaa.2017093
[37]	X. Meng, X. Bao, Z. Zhang, Second order stabilized semi-implicit scheme for the Cahn–Hilliard model with dynamic boundary conditions, J. Comput. Appl. Math., 428 (2023), 115145. https://doi.org/10.1016/j.cam.2023.115145 doi: 10.1016/j.cam.2023.115145
[38]	J. M. Park, R. Mauri, P. D. Anderson, Phase separation of viscous ternary liquid mixtures, Chem. Eng. Sci., 80 (2023), 270–278. https://doi.org/10.1016/j.ces.2012.06.017 doi: 10.1016/j.ces.2012.06.017
[39]	A. Lamorgese, R. Mauri, Diffusion-driven dissolution or growth of a liquid drop embedded in a continuous phase of another liquid via phase-field ternary mixture model, Langmuir, 33 (2017), 13125–13132. https://doi.org/10.1021/acs.langmuir.7b02105 doi: 10.1021/acs.langmuir.7b02105
[40]	A. Lamorgese, R. Mauri, Dissolution or growth of a liquid drop via phase-field ternary mixture model based on the non-random, two-liquid equation, Langmuir, 20 (2018), 125. https://doi.org/10.3390/e20020125 doi: 10.3390/e20020125
[41]	Y. Li, J. Choi, J. Kim, Multi-component Cahn–Hilliard system with different boundary conditions in complex domains, J. Comput. Phys., 323 (2016), 1–16. https://doi.org/10.1016/j.jcp.2016.07.017 doi: 10.1016/j.jcp.2016.07.017
[42]	J. Yang, Y. Li, J. Kim, Modified multi-phase diffuse-interface model for compound droplets in contact with solid, J. Comput. Phys., 491 (2023), 112345. https://doi.org/10.1016/j.jcp.2023.112345 doi: 10.1016/j.jcp.2023.112345
[43]	C. Zhang, H. Ding, P. Gao, Y. Wu, Diffuse interface simulation of ternary fluids in contact with solid, J. Comput. Phys., 309 (2016), 37–51. https://doi.org/10.1016/j.jcp.2015.12.054 doi: 10.1016/j.jcp.2015.12.054
[44]	L. Onsager, Reciprocal relations in irreversible processes. Ⅰ., Phys. Rev., 37 (1931), 405–426. https://doi.org/10.1103/PhysRev.37.405 doi: 10.1103/PhysRev.37.405
[45]	L. Onsager, Reciprocal relations in irreversible processes. Ⅱ., Phys. Rev., 38 (1931), 2265–2279. https://doi.org/10.1103/PhysRev.38.2265 doi: 10.1103/PhysRev.38.2265
[46]	X. Yang, J. Li, G. Forest, Q. Wang, Hydrodynamic theories for flows of active liquid crystals and the generalized Onsager principle, Entropy, 18 (2016), 202. https://doi.org/10.1103/PhysRev.38.2265 doi: 10.1103/PhysRev.38.2265
[47]	P. Sheng, J. Zhang, C. Liu, Onsager principle and electrorheological fluid dynamics, Prog. Theor. Phys. Suppl., 175 (2008), 131–143. https://doi.org/10.1143/PTPS.175.131 doi: 10.1143/PTPS.175.131
[48]	M. Doi, Onsager principle in polymer dynamics, Prog. Polym. Sci., 112 (2021), 101339. https://doi.org/10.1016/j.progpolymsci.2020.101339 doi: 10.1016/j.progpolymsci.2020.101339
[49]	F. Jülicher, J. Prost, Generic theory of colloidal transport, Eur. Phys. J. E, 29 (2009), 27–36. https://doi.org/10.1140/epje/i2008-10446-8 doi: 10.1140/epje/i2008-10446-8
[50]	X. Yang, J. Zhao, Q. Wang, J. Shen, Numerical approximations for a three-component Cahn–Hilliard phase-field model based on the invariant energy quadratization method, Math. Models Methods Appl. Sci., 27 (2017), 1993–2030. https://doi.org/10.1142/S0218202517500373 doi: 10.1142/S0218202517500373
[51]	X. Yang, D. Han, Linearly first- and second-order, unconditionally energy stable schemes for the phase field crystal equation, J. Comput. Phys., 333 (2017), 1116–1134. https://doi.org/10.1016/j.jcp.2016.10.020 doi: 10.1016/j.jcp.2016.10.020
[52]	J. Zhao, X. Yang, J. Li, Q. Wang, Energy stable numerical schemes for a hydrodynamic model of Nematic liquid crystals, SIAM J. Sci. Comput., 38 (2016), A3264–A3290. https://doi.org/10.1137/15M1024093 doi: 10.1137/15M1024093
[53]	J. Zhao, X. Yang, Y. Gong, X. Zhao, X. Yang, J. Li, et al., A general strategy for numerical approximations of non-equilibrium models–Part Ⅰ: Thermodynamical systems, Int. J. Numer. Anal. Mod., 15 (2018), 884–918.
[54]	J. W. Cahn, Critical point wetting, J. Chem. Phys., 66 (1977), 3667–3672. https://doi.org/10.1063/1.434402 doi: 10.1063/1.434402
[55]	X. Zhao, Q. Wang, A second order fully-discrete linear energy stable scheme for a binary compressible viscous fluid model, J. Comput. Phys., 395 (2019), 382–409. https://doi.org/10.1016/j.jcp.2019.06.030 doi: 10.1016/j.jcp.2019.06.030

Reader Comments

Your name:*

Email:*
© 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Mathematical Biosciences and Engineering

3.9

Metrics

Article views(1359) PDF downloads(79) Cited by(0)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(15) / Tables(2)

Mathematical Biosciences and Engineering

A ternary mixture model with dynamic boundary conditions

Related Papers:

Abstract

1. Introduction

2. Mathematical model derivation

2.1. Irreversible thermodynamics

2.2. Non-dimensionalization

3. Numerical scheme

3.1. Mathematical reformulation

3.2. Semi-discrete scheme in time

3.3. Fully discrete scheme

4. Physical phenomenon investigation by numerical simulations

4.1. Contact angle is determined by the additive effects of wall-mixture interaction

4.2. Wall-mixture interaction changes the condensates patterns far away from the surface

4.3. Relaxation rate controls the kinetic process in the bulk and surface

5. Conclusions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog

Abstract

1. Introduction

2. Mathematical model derivation

2.1. Irreversible thermodynamics

2.2. Non-dimensionalization

3. Numerical scheme

3.1. Mathematical reformulation

3.2. Semi-discrete scheme in time

3.3. Fully discrete scheme

4. Physical phenomenon investigation by numerical simulations

4.1. Contact angle is determined by the additive effects of wall-mixture interaction

4.2. Wall-mixture interaction changes the condensates patterns far away from the surface

4.3. Relaxation rate controls the kinetic process in the bulk and surface

5. Conclusions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References