3.1. Derivation of the model

We first consider a mixture of two fluids in a domain ${\it\Omega}$ , and we take an arbitrary material volume $V\in {\it\Omega}$ that moves with the mixture. Within the material volume, we define the properties of the binary compressible fluid as

where $M$ , $\boldsymbol{P}$ , $E$ and $S$ are the total mass, momentum, energy and entropy of the mixture, ${\it\rho}(c)$ is the variable density of the mixture, $\boldsymbol{v}$ is the mass-averaged velocity of the mixture, $|\boldsymbol{v}|^{2}/2$ is the kinetic energy per unit mass, $gz$ is the gravitational potential energy per unit mass, $z$ is the $z$ -coordinate, $\hat{u}$ ( ${\hat{s}}$ ) is the internal energy (entropy) per unit mass and $c$ is the phase variable. Substitution of the mass concentration (2.4) into (3.5) gives

where $C$ stands for the constituent mass of fluid $1$ within the material volume $V(t)$ . In phase-field modelling, apart from the classical free energy density for bulk phases, an extra gradient term is typically added into the model to account for the free energy of the diffuse interface (Cahn & Hilliard Reference Cahn and Hilliard1958). Several ways have been suggested to introduce the gradient term into the phase-field model, e.g. by introducing it into the entropy functional (Wang et al. Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993; Antanovskii Reference Antanovskii1995), free energy functional (Lowengrub & Truskinovsky Reference Lowengrub and Truskinovsky1998) or internal energy functional (Anderson et al. Reference Anderson, McFadden and Wheeler1998; Verschueren et al. Reference Verschueren, van de Vosse and Meijer2001). In the present work, as the thermocapillary effects along the interface are investigated, we expect the surface free energy (serving as the surface tension (see § 5.4)) of our phase-field model to be temperature-dependent. Therefore, according to the thermodynamic relations, we introduce the gradient term into both the internal energy and the entropy of our model, such that

where $u$ , $s$ and $f$ stand for the classical parts of the specific internal energy, entropy and free energy separately. Here, $f$ is the Helmholtz free energy. The parts $u^{grad}$ , $s^{grad}$ and $f^{grad}$ are the gradient terms analogous to the Landau–Ginzburg (Ginzburg & Landau Reference Ginzburg and Landau1950) or Cahn–Hilliard (Cahn & Hilliard Reference Cahn and Hilliard1958) gradient energy. It should be noted that these parts are termed as the ‘non-classical’ terms by Anderson & McFadden (Reference Anderson and McFadden1996), who used a phase-field model to study a single compressible fluid with different phases near its critical point. In addition, ${\it\lambda}_{u}$ and ${\it\lambda}_{s}$ are constant parameters, and ${\it\lambda}_{f}(T)$ is a parameter depending on the temperature and will lead to thermocapillary effects along the interface. It should be noted that ${\it\lambda}_{u}$ , ${\it\lambda}_{s}$ and ${\it\lambda}_{f}(T)$ can be further used to relate the surface tension of the phase-field model to that of the sharp-interface model when the phase-field model reduces to its sharp-interface limit (see § 5.4 for details). As $u({\it\rho},s,c)$ is the classical contribution to the specific internal energy $\hat{u}$ , we have the thermodynamic relation

where the subscripts indicate which variables are held constant when the various partial derivatives are taken. This relation states that the heat transfer ( $T\text{d}s$ ), pressure–volume work $(p/{\it\rho}^{2}\text{d}{\it\rho})$ and chemical work ( $(\partial u/\partial c)\text{d}c$ ) all contribute to the changes in the internal energy. Further, we have the thermodynamic relation for the Helmholtz free energy,

Having in mind the relation (3.10), we obtain

such that

Similarly, we assume that the same thermodynamic relations that hold for the classical terms also hold for the general terms, such that

With the relations (3.11) and (3.13), we must also have the relations for the gradient terms,

and for the corresponding coefficients,

For simplicity, we omit all the subscripts in the following derivations. Under the assumptions above, the general forms of the physical balance associated with $M$ , $\boldsymbol{P}$ , $E$ , $S$ and $C$ can be given as follows:

where $\hat{\boldsymbol{n}}$ is the unit outward normal vector of the boundary and $\hat{\boldsymbol{z}}$ is the vertical component of the unit normal vector. Equation (3.17) represents the mass balance of the mixture within the volume. Equation (3.18) represents the momentum balance, stating that the rate of change in total momentum equals the force (surface forces $\unicode[STIX]{x1D662}$ and body forces ${\it\rho}g\hat{\boldsymbol{z}}$ ) acting on the volume. Here, only the gravitational forces are considered. The energy balance equation (3.19) states that the change in total energy equals the rate of work done by the forces ( $\unicode[STIX]{x1D662}$ ) on the boundary plus the energy flux (classical $\boldsymbol{q}_{E}$ and non-classical $\boldsymbol{q}_{E}^{nc}$ internal energy flux) through the boundary. The entropy balance (3.20) states that the rate of change of entropy in the control volume during the process equals the net entropy transfer through the boundary (classical $\boldsymbol{q}_{E}/T$ and non-classical $\boldsymbol{q}_{S}^{nc}$ entropy flux) plus the local entropy generation ( $S_{gen}\geqslant 0$ ) within the control volume (e.g. Moran et al. Reference Moran, Shapiro, Boettner and Bailey2010). Based on the second law of thermodynamics, the local entropy generation is non-negative for a dissipative system (or, say, for an irreversible process), which is key to the thermodynamic frame that we used for the derivations. For the constituent mass balance (3.21), we use (3.6) and Theorem 1 to obtain

Substituting (2.6) into (3.22), we obtain $\boldsymbol{q}_{C}=\tilde{{\it\rho}}_{1}(\boldsymbol{v}_{1}-\boldsymbol{v})$ , where $\boldsymbol{q}_{C}$ stands for the mass flux of fluid 1 with velocity $(\boldsymbol{v}_{1}-\boldsymbol{v})$ through the boundary of the control volume. It should be noted that in the following derivations $q_{C}$ will be related to the chemical potential of the phase field, which is analogous to the standard derivations of the Cahn–Hilliard equations (see, for example, Anderson et al. Reference Anderson, McFadden and Wheeler1998; Lowengrub & Truskinovsky Reference Lowengrub and Truskinovsky1998).

In what follows, we use the definitions (3.1)–(3.5) and the balance laws (3.17)–(3.21) to obtain equations that are expressed in terms of the above unknowns, including $\unicode[STIX]{x1D662}$ , $\boldsymbol{q}_{E}$ , $\boldsymbol{q}_{E}^{nc}$ , $\boldsymbol{q}_{S}^{nc}$ , $\boldsymbol{q}_{C}$ and $S_{gen}$ . We then specify these unknowns with respect to the second law of thermodynamics (ensuring $S_{gen}\geqslant 0$ ) and the concept of thermodynamic consistency of the phase-field model.

For the mass balance (3.17), we use Theorem 1 to obtain

based on which we have the following.

It should be noted that as Theorems 1 and 2 are frequently used, we will not refer to them in the following derivations.

For momentum balance (3.18), we simply have

For energy balance (3.19), we obtain

where (3.10), (3.23) and (3.25) and the following identities are used:

and

Here, ‘ $:$ ’ stands for the double dot product of the stress tensor (e.g. Mase & Mase Reference Mase and Mase1999).

For entropy balance (3.20), we obtain

where, similarly to (3.28), the following identity is used:

For constituent mass balance (3.21), we simply have

We then use (3.29) and (3.31) to substitute the terms ${\it\rho}\text{D}s/\text{D}t$ and ${\it\rho}\text{D}c/\text{D}t$ in (3.26), and use the relation (3.16) to obtain the expression for the entropy generation,

To ensure the non-negativity of the entropy generation $S_{gen}\geqslant 0$ (second law of thermodynamics), we specify the unknown terms in the form

It should be noted that ${\bf\tau}$ is the deviatoric stress tensor from the classical Navier–Stokes equations (e.g. Batchelor Reference Batchelor2000). Here, we use the thermodynamic relation (3.13) to obtain the chemical potential ${\it\mu}_{C}$ . The pressure $p$ can be obtained immediately through the thermodynamic relation (3.13). By substituting the above terms into (3.17)–(3.21), we obtain the system of equations for the phase-field model governing binary compressible flows with thermocapillary effects, It should be noted that the second term in the stress tensor $\unicode[STIX]{x1D662}$ is the extra reactive stress (Ericksen’s stress) to mimic the surface tension. This stress term is associated with the presence of concentration gradient energy (Cahn–Hilliard energy). We note that the temperature-dependent coefficient in $\unicode[STIX]{x1D662}$ is a linear function of temperature, which leads to thermocapillary effects along the interface (see § 4.2 for details). Here, $m_{C}$ is a positive constant standing for the mobility of the diffuse interface. It should be noted that in the non-classical heat (or entropy) flux $\boldsymbol{q}_{E}^{nc}$ (or $\boldsymbol{q}_{S}^{nc}$ ), the term ${\it\rho}{\it\lambda}_{u}\boldsymbol{{\rm\nabla}}c~\text{D}c/\text{D}t$ (or ${\it\rho}{\it\lambda}_{s}\boldsymbol{{\rm\nabla}}c~\text{D}c/\text{D}t$ ) is associated with the gradient energy (or entropy respectively) and is in the direction of the gradient of the phase variable. Similar terms were obtained by Wang et al. (Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993), who used a phase-field model to study the solidification of single material, and by Anderson & McFadden (Reference Anderson and McFadden1996), who used a phase-field model to study a single compressible fluid with different phases near its critical point. In addition, a non-classical energy flux term $m_{C}{\it\mu}_{C}\boldsymbol{{\rm\nabla}}{\it\mu}_{C}$ appears in our energy balance equation (3.40). The same energy flux term was obtained by Gurtin et al. (Reference Gurtin, Poligone and Vinale1996) (see (28)), who re-derived the model H in the framework of classical continuum mechanics. A ‘counterpart’ entropy flux term was identified by Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998) when deriving a phase-field model for binary compressible fluid, where this term is required to keep the model compatible with the second law of thermodynamics. In the latter work, the isothermal fluid flow was studied, so that the temperature $T$ in the entropy flux was treated as constant, whereas in our work, the temperature is not a constant as the thermocapillary effects are considered here. They identified this non-classical term as the entropy flux transported through the boundary by chemical diffusion. Our model agrees with these works well and therefore we identify this non-classical energy flux term as the energy that is carried into the control volume by the chemical diffusion. It should be noted that several phase-field models (e.g. Blesgen Reference Blesgen1999; Allaire, Clerc & Kokh Reference Allaire, Clerc and Kokh2002; Abels & Feireisl Reference Abels and Feireisl2008) have been presented to study binary compressible fluids, where the specifications of the free energy (3.37) that contributes to the compressibility of binary compressible fluids are discussed and provided.

Similarly to the approach that defines the variable density ${\it\rho}(c)$  (2.5), we define the variable viscosity ${\it\mu}(c)$ and the variable thermal diffusivity $k(c)$ for the mixture in the form of the harmonic average,

where ${\it\mu}_{i}$ and $k_{i}$ are the viscosity and thermal conductivity of the $i$ th fluid.

3.2. Thermodynamic consistency and Galilean invariance

As our phase-field model (3.38)–(3.43) is derived within a thermodynamic framework, it implies that the first and second thermodynamic laws are naturally underlying the model. However, from the numerical point of view, thermodynamic consistency can further serve as a criterion to design the numerical methods. In our phase-field model, the Navier–Stokes equations are coupled with the Cahn–Hilliard equations and energy balance equation, which leads to a nonlinear system. Moreover, as rapid variations in the solutions of the phase variable occur near the interfacial region, the energy stability of the numerical method is critical. Recently, the preservation of the thermodynamic laws at the discrete level has been reported to play an important role in the design of numerical methods (e.g. Lin & Liu Reference Lin and Liu2006; Lin, Liu & Zhang Reference Lin, Liu and Zhang2007 for liquid crystal models; Hua et al. Reference Hua, Lin, Liu and Wang2011; Guo et al. Reference Guo, Lin and Lowengrub2014a for phase-field models), which not only immediately implies the stability of the numerical scheme, but also ensures the correctness of the solutions. Hence, in contrast to the derivations, we now show that the first and second laws of thermodynamics can be derived from the system of (3.38)–(3.43), which can be further used to design the numerical methods. In addition, important modelling properties, Onsager reciprocal relations and Galilean invariance will be verified as well.

4.1. Derivation of the model

In order to study situations in which the density in each phase is uniform, it is convenient to adopt a thermodynamic formulation which does not employ the density as an independent variable, as in the model of quasi-incompressible flow considered by Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998). Following their work, we choose the pressure and temperature as independent variables, and work with the Gibbs free energy. In addition, for a binary incompressible fluid system, the free energy density can appear as the ‘per unit mass’ quantity or ‘per unit volume’ quantity. In most phase-field models for two-phase flows (e.g. Hohenberg & Halperin Reference Hohenberg and Halperin1977; Liu & Shen Reference Liu and Shen2002), the densities of the two components are assumed to be constant and equal, and the per unit mass and per unit volume specifications of the free energy density are equivalent. However, in the situation we study here, the densities of the two fluids of the mixture may not be matched and thus the per unit mass and per unit volume forms are not equivalent. As mentioned above, several models have been developed for binary incompressible fluids with different densities, in which the per unit volume form of the free energy density was employed by Boyer (Reference Boyer2002), Ding et al. (Reference Ding, Spelt and Shu2007), Shen & Yang (Reference Shen and Yang2010) and Abels et al. (Reference Abels, Garcke and Grün2012) and the per unit mass form by Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998). Here, we concentrate on the Gibbs free energy density in the per unit mass form, and denote it by ${\hat{g}}(T,p,c,\boldsymbol{{\rm\nabla}}c)$ . Again, similarly to the definition of the free energy (3.9) for binary compressible fluid, we introduce the gradient terms (gradient energy) into the Gibbs free energy of our model, which can then be given in the form

where $g$ is the classical part of the Gibbs free energy density and ${\it\lambda}_{f}(T)$ is a temperature-dependent coefficient and will lead to thermocapillary effects along the interface (see § 4.3 for details). For the classical part of the internal energy defined by (3.7), we have the following thermodynamic relation:

Use of the thermodynamic relation (3.10) leads to

where we note the relations

Here, as we notice that the variable density is independent of temperature and pressure (see (2.5)), the condition of incompressibility can then be written in the terms of the Gibbs free energy,

where the second condition in (4.4) is used. Condition (4.5) implies that the Gibbs free energy is a linear function of pressure (e.g. Lowengrub & Truskinovsky Reference Lowengrub and Truskinovsky1998),

We then redefine the classical internal energy as a function of $T$ , $p$ and $c$ ,

where the relations (4.3) and (4.4) still hold. Similarly to the definition of the internal energy (3.7) and entropy (3.8) for the binary compressible fluid model, the specific internal energy $\hat{u}$ and the specific entropy ${\hat{s}}$ for binary incompressible fluids can be redefined in the form

where $\tilde{u}$ and $\tilde{s}$ are the classical parts of the specific internal energy and entropy associated with the Gibbs free energy, and ${\it\lambda}_{u}$ and ${\it\lambda}_{s}$ are constant. In addition to these classical contributions, we assume that the same thermodynamic relations that hold for the classical terms also hold for the total terms, such that

Thus, from (4.1), (4.7)–(4.9), the relation for the coefficients (3.16) holds as well. The specifications of these three coefficients will be discussed in § 4.2. It should be noted that ${\it\lambda}_{u}$ and ${\it\lambda}_{s}$ together with ${\it\lambda}_{f}(T)$ (in (4.1)) can be further used to relate the surface tension of the phase-field model to that of the sharp-interface model when our phase-field model reduces to its sharp-interface limit (see § 5.4 for details).

Now we derive the system of equations for the quasi-incompressible phase-field model. We still use (3.1)–(3.5) to define the total properties, namely mass $M$ , momentum $\boldsymbol{P}$ , energy $E$ , entropy $S$ and mass constituent $C$ in a material control volume $V(t)$ of the domain ${\it\Omega}$ . We further assume that the corresponding general balance laws (3.17)–(3.21) that hold for the binary compressible fluid also hold for the quasi-incompressible fluid, which can then be written as

where, as the pressure $p$ is not defined in the traditional way, the general stress tensor $\unicode[STIX]{x1D662}$ is not defined explicitly.

It should be noted that, in contrast to the case of binary compressible fluid, the classical part of the internal energy $\tilde{u}$ we defined here does not depend on the entropy $\tilde{s}$ directly. However, as the derivations are carried out within the thermodynamic framework which is based on the entropy generation, a thermodynamic relation between the internal energy and the entropy is still needed. Having in mind the definition of the internal energy (4.7) and using the relations (4.3) and (4.4), we obtain the following relation between the classical part of the internal energy $\tilde{u}$ , the Gibbs free energy $g$ and the entropy $\tilde{s}$ :

Then similarly to the method used for the binary compressible fluid model, we use the unknowns, including $\unicode[STIX]{x1D662}$ , $\boldsymbol{q}_{E}$ , $\boldsymbol{q}_{E}^{nc}$ and $\boldsymbol{q}_{S}^{nc}$ , to express the entropy generation in the form

where we have used (4.11)–(4.16). Here, $\tilde{{\it\mu}}_{C}=\partial g_{0}(c)/\partial c-{\it\lambda}_{f}(T){\rm\Delta}c$ is a potential term. As the pressure is no longer defined by the thermodynamic formulae in this model, we now derive the pressure in an alternative way that was used by Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998), where the pressure was obtained from the non-dissipated part of the general stress $\unicode[STIX]{x1D662}$ . Considering a dissipative process, we denote the general stress tensor by $\unicode[STIX]{x1D662}=\unicode[STIX]{x1D662}_{0}+{\bf\tau}$ , in which ${\bf\tau}$ is the deviatoric stress tensor with zero trace, and $\unicode[STIX]{x1D662}_{0}$ is the unknown part to be defined later. We then denote $\boldsymbol{D}\boldsymbol{v}=\boldsymbol{{\rm\nabla}}\boldsymbol{v}-(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v})\unicode[STIX]{x1D644}/3$ as the deviatoric part of $\boldsymbol{{\rm\nabla}}\boldsymbol{v}$ ( $\text{tr}\,\boldsymbol{D}\boldsymbol{v}=0$ ). The entropy expression (4.17) can be rewritten as

where we have used the mass balance (4.11) and the following identity:

Now we assume that the first two terms on the right-hand side of (4.18) are non-dissipative and define the pressure $p$ by

such that

in which the way we use to define the pressure in (4.20) is analogous to the way that defines the kinematic pressure for the classical Navier–Stokes equations (Batchelor Reference Batchelor2000). To ensure that our model is consistent with the second law of thermodynamics ( $S_{gen}\geqslant 0$ ), we specify the unknown terms as follows:

Here, $\unicode[STIX]{x1D662}=\unicode[STIX]{x1D662}_{0}+{\bf\tau}$ . It should be noted that, comparing with the model developed by Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998), an extra term, ${\it\rho}{\it\lambda}_{f}(T)|\boldsymbol{{\rm\nabla}}c|^{2}\unicode[STIX]{x1D644}$ , appears in the stress tensor $\unicode[STIX]{x1D662}_{0}$  (4.23). This term could be absorbed into the pressure for convenience (see Hua et al. Reference Hua, Lin, Liu and Wang2011; Liu & Shen Reference Liu and Shen2002 for examples). However, we still keep this term in order that the surface gradient of the surface tension in the sharp-interface limit can be recovered (see § 5.4 for details). This treatment is similar to that used by Guo et al. (Reference Guo, Lin and Wang2014b ) and Garcke et al. (Reference Garcke, Hinze and Kahle2014), where in Garcke et al. (Reference Garcke, Hinze and Kahle2014) this term is also required to recover the surface gradient term in the asymptotic analysis. Apart from the momentum equation, the pressure appears in the chemical potential equation as well, which is different from the chemical potential (3.43) for the binary compressible fluid model. By substituting the above terms into (4.11)–(4.15), we obtain the system of equations for the phase-field model governing the quasi-incompressible fluid with thermocapillary effects, By multiplying (4.26), (4.27) and (4.29)–(4.31) by $p/{\it\rho}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{v}/2+u$ , $\boldsymbol{v}$ , $T$ , ${\it\mu}_{C}$ and ${\it\rho}\text{D}c/\text{D}t$ , and summing them up, we can obtain the first law of thermodynamics (3.19) that we used to derive the model. By substituting the terms, including $\unicode[STIX]{x1D662}$ , $\boldsymbol{q}_{E}$ , $\boldsymbol{q}_{E}^{nc}$ , $\boldsymbol{q}_{S}^{nc}$ and $\boldsymbol{q}_{C}$ , into the entropy generation (4.18), we obtain the second law of thermodynamics for our phase-field model,

Similarly to the binary compressible model, the choices of the terms ${\bf\tau},\boldsymbol{q}_{E}$ and $\boldsymbol{q}_{C}$ satisfy the linear relation (3.48) and the Onsager reciprocal relations (§ 3.2.2). Moreover, it can be observed that the entropy generation (4.32) and the system equations are Galilean invariant.

As mentioned above, several phase-field models have been developed for two-phase flows with thermocapillary effects. However, in most of these models, the classical energy balance equation,

was incorporated directly into the phase-field model, where thermodynamic consistency can hardly be maintained. Compared with the classical energy balance equation (4.33), several extra terms appear in our energy balance equation (4.28), which guarantee thermodynamic consistency (see § 4.2).

It should be noted that, if we define a new pressure as $\bar{p}=p-{\it\rho}{\it\lambda}_{f}(T)|\boldsymbol{{\rm\nabla}}c|^{2}$ , and substitute it into the system (4.26)–(4.31), our model, in the isothermal case, reduces to the quasi-incompressible NSCH model developed by Lowengrub & Truskinovsky (Reference Lowengrub and Truskinovsky1998).

By using the variable mass density (2.5), the mass balance equation (4.26) can be further rewritten as

where we have used the Cahn–Hilliard equation (4.30) and let ${\it\alpha}=({\it\rho}_{2}-{\it\rho}_{1})/{\it\rho}_{2}{\it\rho}_{1}$ .

The initial conditions are given by

For the velocity, the usual no-slip boundary conditions can be posed on $\partial {\it\Omega}$ ,

For the phase field, it is normal to employ Neumann boundary conditions on $\partial {\it\Omega}$ ,

For the temperature, Dirichlet and Neumann boundary conditions can be posed on $\partial {\it\Omega}$ ,

for the specified temperature and heat flux on the boundary $\partial {\it\Omega}$ respectively, and Robin boundary conditions can be posed as well.

4.2. Specifications of the model

We now specify the properties including the Gibbs free energy, entropy and internal energy for our phase-field model (4.26)–(4.31). In Anderson, McFadden & Wheeler (Reference Anderson, McFadden and Wheeler2000), a phase-field model for the solidification of a pure material that includes convection in the liquid phase was developed, in which the case of quasi-incompressibility (assuming that the density in each phase is uniform) was discussed. In their work, the Gibbs free energy was suggested in the form

which we have adopted for the present work. Here, $c_{hc}$ is the heat capacity, $T_{0}$ is the reference temperature, $\tilde{u} _{0}$ is the reference internal energy corresponding to $T_{0}$ and ${\it\gamma}(T)$ is a temperature-dependent parameter that will be discussed later in this section. The free energy function $h(c)$ is a double-well potential and is given by

where the wells define the phases, and lead to an interfacial layer with large variations for $c$ (e.g. Gurtin et al. Reference Gurtin, Poligone and Vinale1996). It should be noted that the form for ${\hat{g}}$  (4.39) is consistent with the incompressible condition (4.6), which is a linear function of pressure. Moreover, this form for ${\hat{g}}$ is consistent with an internal energy $\hat{u}$ , which is a linear function of temperature and leads to the classical heat equation in the bulk liquid (Wang et al. Reference Wang, Sekerka, Wheeler, Murray, Coriell, Braun and McFadden1993; Anderson, McFadden & Wheeler Reference Anderson, McFadden and Wheeler2001). The corresponding expressions for the entropy and internal energy are assumed in the form

where $\tilde{u} _{0}$ is the reference internal energy corresponding to $T_{0}$ .

We now specify those coefficients, including ${\it\lambda}_{f}(T)$ , ${\it\lambda}_{s}$ , ${\it\lambda}_{u}$ , ${\it\gamma}_{f}(T)$ , ${\it\gamma}_{s}$ and ${\it\gamma}_{u}$ , that are used to define the internal energy, entropy and free energy of the system (4.39)–(4.43). In the sharp-interface model for thermocapillary flow, the interface is usually represented as a surface of zero thickness with the surface tension as its physical property. An equation of state is required to relate the surface tension to the temperature, where, for the sake of simplicity, we only consider a linear relation in this study,

where ${\it\sigma}_{0}$ is the interfacial tension at the reference temperature $T_{0}$ , and ${\it\sigma}_{T}$ is the rate of change of interfacial tension with temperature, defined as ${\it\sigma}_{T}=\partial {\it\sigma}(T)/\partial T$ . In our phase-field model, however, the interface has finite thickness and the extra reactive stress (Ericksen’s stress) $\unicode[STIX]{x1D64F}$  (4.25) appears in the Navier–Stokes equation to mimic the surface tension, where the coefficient of $\unicode[STIX]{x1D64F}$ ,

is a linear function of temperature. We then try to relate ${\it\sigma}(T)$ and ${\it\lambda}(T)$ by introducing two parameters. The first parameter is ${\it\epsilon}$ with respect to the diffuse-interface thickness and the second one is ${\it\eta}$ , a ratio parameter that relates the two surface tensions. As the interface thickness goes to zero, our phase-field model reduces to its sharp-interface limit, and the value of ${\it\eta}$ can then be determined (see § 5.4 for details). The corresponding coefficients can then be given as

Here, the coefficients ${\it\lambda}_{f}(T)$ , ${\it\lambda}_{s}$ and ${\it\lambda}_{u}$ for the gradient terms are of $O({\it\epsilon}^{2})$ of the coefficients ${\it\gamma}_{f}(T)$ , ${\it\gamma}_{s}$ and ${\it\gamma}_{u}$ for the corresponding classical terms, which agrees with the definition of the Cahn–Hilliard free energy (e.g. Lowengrub & Truskinovsky Reference Lowengrub and Truskinovsky1998; Liu & Shen Reference Liu and Shen2002). With the specifications above, the total energy $E$ of our phase-field model can now be written as

4.3. Non-dimensionalization

With the help of the specification in (4.8), we non-dimensionalize the phase-field model (4.26)–(4.28), (4.30) and (4.31) as follows. We let $L^{\star }$ , $V^{\star }$ and $T^{\star }$ denote the characteristic scales of length, velocity and temperature. Then we introduce the dimensionless variables $\bar{\boldsymbol{x}}=\boldsymbol{x}/L^{\star }$ , $\bar{t}=V^{\star }t/L^{\star }$ , and also $\bar{{\it\epsilon}}={\it\epsilon}/L^{\star }$ , $\bar{\boldsymbol{v}}=\boldsymbol{v}/V^{\star }$ , $\bar{p}=p{\it\rho}_{1}{\it\mu}_{C}^{\star }$ , $\bar{{\it\mu}}_{C}={\it\mu}_{C}/{\it\mu}_{C}^{\star }$ . For the variable density ${\it\rho}(c)$ , viscosity ${\it\mu}(c)$ and thermal conductivity $k(c)$ , (2.5) and (3.44), we let ${\it\rho}_{i}$ , ${\it\mu}_{i}$ and $k_{i}\;(i=1,2)$ denote the corresponding properties of the $i$ th fluid, and introduce the dimensionless variables $\bar{{\it\rho}}_{r}={\it\rho}(c)/{\it\rho}_{1}$ , $\bar{{\it\mu}}_{r}={\it\mu}(c)/{\it\mu}_{1}$ and $\bar{k}_{r}=k(c)/k_{1}$ . Moreover, for the temperature field, we introduce a new dimensionless variable $\bar{T}=(T-T_{0})/T^{\star }$ . The surface tension ${\it\sigma}(T)$  (4.44) is scaled by ${\it\sigma}_{0}$ such that $\bar{{\it\sigma}}(T)={\it\sigma}(T)/{\it\sigma}_{0}$ . Then ${\it\sigma}_{T}$ is scaled by ${\it\sigma}_{0}/T^{\star }$ , such that $\bar{{\it\sigma}_{T}}={\it\sigma}_{T}T^{\star }/{\it\sigma}_{0}$ . Omitting the bar notation, our phase-field model can now be rewritten as

where $M=V^{2}/{\it\mu}_{C}$ is an analogue of the Mach number, $Ca={\it\eta}{\it\epsilon}{\it\sigma}_{0}/{\it\mu}_{C}L$ is the capillary number, which measures the thickness of the interface, $Re={\it\mu}_{1}/{\it\rho}_{1}VL$ is the Reynolds number, $Fr=V^{2}/gL$ is the Froude number, $Pe={\it\rho}LV/m_{C}{\it\mu}_{C}$ is the diffusional Peclet number, $Ec=V^{2}/c_{ch}T$ is the Eckert number, which characterizes energy dissipation, and $Ma={\it\rho}c_{ch}VL/k_{1}$ is the Marangoni number. It should be noted that these non-dimensional system equations will be computed to study the effects of the Marangoni number through the example of thermocapillary migration of a drop. See § 6.5 for details.

5.1. Pillbox argument

In this section, we apply a pillbox argument to our phase-field model (4.26)–(4.31). In contrast to the sharp-interface model, the interface of the phase-field model is diffusive with finite thickness $O({\it\epsilon})$ . The phase variable (here the mass concentration $c$ ) is chosen to characterize the different phases, which take distinct values (here $c=0,1$ ) for the different phases, and change rapidly through the interfacial region. Within this interfacial region, we chose a contour line of $c$ (here $c=0.5$ ) to represent the dividing surface ${\it\Gamma}$ for the following derivations (see Gibbs Reference Gibbs1928; Everett Reference Everett1972; Rowlinson & Widom Reference Rowlinson and Widom1982 for details of the dividing surface). Moreover, as the largest variations of the phase variable occur in the direction normal to the interface, the side faces of the pillbox need to be treated carefully. Figure 1 shows the pillbox-shaped control volume designed for our phase-field model, where the surface is divided into three parts, namely the top $S_{top}$ , bottom $S_{bot}$ and side $S_{side}$ surfaces with their unit normal vectors $\hat{\boldsymbol{n}}_{T}$ , $\hat{\boldsymbol{n}}_{B}$ and $\hat{\boldsymbol{n}}_{S}$ respectively. The volume of the pillbox is $V=V_{1}+V_{2}$ , where $V_{i}$ is the volume of a single component. The pillbox has a thickness of $2{\it\delta}$ , where the top of the pillbox is above the dividing surface ${\it\Gamma}$ at a height ${\it\zeta}={\it\delta}$ and the bottom is below ${\it\Gamma}$ at a height ${\it\zeta}=-{\it\delta}$ . Here, ${\it\zeta}$ is a local coordinate normal to the interface ${\it\Gamma}$ . In addition, the pillbox contains a portion of the diffuse interface with thickness $O({\it\epsilon})$ , in which ${\it\Gamma}$ stands for the dividing surface with its normal and tangent unit vectors $\hat{\boldsymbol{n}}_{I}$ and $\hat{\boldsymbol{m}}_{I}$ . The key limit in the pillbox argument is that ${\it\epsilon}\ll {\it\delta}\ll L$ , where $L$ is a length scale associated with the outer flow. In this limit, the volume of the pillbox becomes negligible on the outer scales, but the variations in the concentration variable $c$ that defines the interfacial region occur over a region fully contained within the pillbox. Also in this limit, the top ( $S_{top}$ ) and bottom surfaces ( $S_{bot}$ ) of the pillbox collapse onto the interface ${\it\Gamma}$ , and have normal vectors with opposite directions,

Moreover, we assume that the dividing interface is moving with the velocity $\boldsymbol{v}_{I}$ (Anderson & McFadden Reference Anderson and McFadden1996; Anderson et al. Reference Anderson, McFadden and Wheeler1998).

Figure 1. A schematic diagram showing a diffuse interface between two fluids intersecting with a pillbox-shaped control volume. Here, $\hat{\boldsymbol{n}}_{T}$ , $\hat{\boldsymbol{n}}_{B}$ and $\hat{\boldsymbol{n}}_{S}$ stand for the unit normal vector of the pillbox boundary on its top, bottom and side respectively. The dotted lines represent the diffuse interface with thickness $O({\it\epsilon})$ . The thickness of the pillbox is $2{\it\delta}$ . In the limit ${\it\epsilon}\ll {\it\delta}\ll L$ , the interface thickness goes to zero, and the interface has constant density. Here, $\hat{\boldsymbol{n}}_{I}$ and $\hat{\boldsymbol{m}}_{I}$ stand for the unit normal and tangent vectors of the interface.

5.2. Governing equations in the sharp-interface limit

We first derive the system of equations in bulk regions away from the interfacial region. Here, we only concentrate on the equations of mass, momentum and energy balance. The system of (4.26)–(4.28) reduces to the classical equations appropriate for incompressible flows in bulk regions,

where ${\it\rho}_{i}$ , ${\it\mu}_{i}$ and $k_{i}$ are the corresponding physical properties for the $i$ th fluid. We now seek to derive the jump conditions for (5.2)–(5.4) at the interface from our phase-field model (4.26)–(4.31).

5.3. Jump condition for mass balance

In the limit ${\it\epsilon}\ll {\it\delta}\ll L$ , we have the properties (Anderson & McFadden Reference Anderson and McFadden1996; Anderson et al. Reference Anderson, McFadden and Wheeler1998)

On substituting (5.5) into the integral of equation (4.11) and using the divergence theorem we obtain

According to our pillbox argument, we break up the above surface integral into pieces for the top, bottom and side surfaces to obtain

Here, the surface integral of the side portion is further written in terms of a line integral on the surface and an integral in the normal direction $\hat{\boldsymbol{n}}_{S}$ , where the line is a closed curve at the side of the control volume that parallel to the interface. For viscous fluid under normal operating conditions, it is an experimentally observed fact (like the no-slip boundary conditions at solid walls) that no slip takes place at the interface (Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011). Therefore, in the limit ${\it\epsilon}\ll {\it\delta}\ll L$ , we have

This condition implies that the third term on the left in (5.8) is bounded and does not contribute to the integral. Equation (5.8) can be reduced to

where $[{\it\chi}]={\it\chi}_{2}-{\it\chi}_{1}$ refers to the jump of the quantity ${\it\chi}$ across the singular interface. Since the pillbox control volume $V$ that contains a portion of the diffuse interface is arbitrary, the integrand in (5.10) must be zero. This then yields the mass balance jump condition at the interface in a two-phase fluid system,

Further, if we assume that there is no phase change (i.e. no flux) across the interface, (5.11) reduces to the jump condition that

5.4. Jump condition for momentum balance

By substituting (5.6) into the integral of momentum equation (4.27), we obtain

where we have used the mass balance equation (4.26), such that

Moreover, in the limit ${\it\epsilon}\ll {\it\delta}\ll L$ , we assume that the gravitational term ${\it\rho}g\hat{\boldsymbol{z}}$ is bounded and thus does not contribute to the volume integral. We then break up the above surface integral into pieces for the top, bottom and sides of the pillbox to obtain

We assume that the most rapid variations in the phase field take place across the interfacial region with the direction normal to the interface ${\it\Gamma}$ . In the limit ${\it\epsilon}\ll {\it\delta}\ll L$ , local to the interface we have (Anderson & McFadden Reference Anderson and McFadden1996; Lowengrub & Truskinovsky Reference Lowengrub and Truskinovsky1998)

such that

Condition (5.9) implies that the fluid velocity term ${\it\rho}\boldsymbol{v}\otimes (\boldsymbol{v}-\boldsymbol{v}_{I})\boldsymbol{\cdot }\hat{\boldsymbol{n}}_{S}$ is bounded and does not contribute to the integral over the side surface of the pillbox. The terms $-{\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{v}+\boldsymbol{{\rm\nabla}}\boldsymbol{v}^{\text{T}})\boldsymbol{\cdot }\hat{\boldsymbol{n}}_{S}$ are bounded and do not contribute to the side integral. We argue that the term $2/3{\it\mu}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v})$ is bounded across the interfacial region and thus does not contribute to the side integral. The pressure $p$ is bounded and does not contribute to the side integral. Further, the non-classical stress term $\unicode[STIX]{x1D64F}$ does not contribute to the integral over the top and bottom surfaces. Equation (5.15) reduces to

where the condition (5.1) is used. Here, for our pillbox argument to make sense, we require that within the pillbox the temperature is continuous and the variations are small over a small distance (of the order of the pillbox thickness ${\it\delta}$ ). In the limit ${\it\epsilon}\ll {\it\delta}\ll L$ , the temperature $T$ is approximately uniform along the direction normal to the interface. It should be noted that a similar assumption for the temperature was also suggested by Jasnow & Vinals (Reference Jasnow and Vinals1996), where a surface tension term with thermocapillary effects was identified from a phase-field model in its sharp-interface limit. Denoting the surface tension by

and substituting into (5.18), we obtain

where, in the limit ${\it\epsilon}\ll {\it\delta}\ll L$ , we assume that the tangential unit vector $\hat{\boldsymbol{m}}_{I}$ is independent of ${\it\zeta}$ and thus can be taken out of the line integral. Use of the surface divergence theorem (Weatherburn Reference Weatherburn1939) leads to

By substituting (5.21) into (5.20), we obtain

Here, $\boldsymbol{{\rm\nabla}}_{s}$ is the surface gradient and ${\it\kappa}=-\boldsymbol{{\rm\nabla}}_{s}\boldsymbol{\cdot }\hat{\boldsymbol{n}}_{I}$ is the mean curvature of the surface (e.g. Weatherburn Reference Weatherburn1939). The first term on the right is the tangential thermocapillary (Marangoni) force that accounts for the non-uniform surface tension, while the second is the normal surface tension force. Again, if we assume that there is no phase change (i.e. no flux) across the interface, (5.22) reduces to the jump condition that

which is the classical momentum balance jump condition at the interface for two-phase incompressible fluid with thermocapillary effects.

Figure 2. The stationary solution $c_{0}$ (solid line) for the phase field. Here, $A$ is a point on the dividing surface ${\it\Gamma}$ , and ${\it\delta}$ and $-{\it\delta}$ are the positions of the top and bottom surfaces of the pillbox (dotted lines).

It should be noted that we can relate the surface tension of our phase-field model $\tilde{{\it\sigma}}(T)$ (identified in (5.19)) to that of the sharp-interface model ${\it\sigma}(T)$ (defined in (4.44)) by letting

The value of the ratio parameter ${\it\eta}$ can then be determined through the following equation:

It has been argued by Chella & Vinals (Reference Chella and Vinals1996) that in the limit of a gently curved interface, and when the motion of the interface is slow, the phase variable $c$ can be approximated by its 1D stationary solution $c_{0}$ along the direction normal to the interface. For simplicity, we now assume that the local coordinate ${\it\zeta}$ coincides with the $y$ direction, and the position of the dividing surface is $y_{0}=0$ . In the 1D case, we have the following stationary solution $c_{0}$ near the interfacial region:

which is shown in figure 2. Here, $y={\it\delta}$ and $y=-{\it\delta}$ are the positions of the top and bottom surfaces of the pillbox separately. In the limit ${\it\epsilon}\ll {\it\delta}\ll L$ , we note the conditions

By substituting (5.26) and the variable density (2.5) into (5.25) we obtain

where the condition (5.27) is used. It should be noted that for the density matched case ( ${\it\rho}_{1}={\it\rho}_{2}$ ), (5.25) leads to a simpler expression for ${\it\eta}$ , which is

This agrees with the result obtained by Rowlinson & Widom (Reference Rowlinson and Widom1982), Yue et al. (Reference Yue, Feng, Liu and Shen2004) and Aland (Reference Aland2012). In § 6, we will compute some examples by using our phase-field model for quasi-incompressible fluids (4.26)–(4.31).

5.5. Jump condition for energy balance

To derive the jump condition for the energy balance at the interface, we first substitute the terms $E$ , $\unicode[STIX]{x1D662}$ , $\boldsymbol{q}_{E}$ and $\boldsymbol{q}_{E}^{nc}$  (3.3), (4.22), (4.23) and (4.25) into the energy balance equation (3.19). In the integral form, we obtain

where we have used the identities

and the following properties which are similar to (5.5) and (5.6): We then break up the above integral (5.30) into pieces for the top, bottom and sides of the pillbox to obtain

where we assume that the heat capacity $c_{hc}$ is a constant. In the limit ${\it\epsilon}\ll {\it\delta}\ll L$ , the non-classical terms of the internal energy $\hat{u}$  (4.43), $\unicode[STIX]{x1D64F}$ , ${\it\lambda}_{u}({\it\rho}\boldsymbol{{\rm\nabla}}c(\text{D}c/\text{D}t))$ and $m_{C}{\it\mu}_{C}\boldsymbol{{\rm\nabla}}{\it\mu}_{C}$ , do not contribute to the top and the bottom surface integrals. The non-classical term $m_{C}{\it\mu}_{C}\boldsymbol{{\rm\nabla}}{\it\mu}_{C}$ and the energy term $k\boldsymbol{{\rm\nabla}}T$ are bounded in the tangential direction and do not contribute to the side integral. Equation (5.35) then reduces to

where in the last term of (5.36), we argue that the interface velocity $\boldsymbol{v}_{I}$ is independent of the local coordinate ${\it\zeta}$ and thus can be taken out of the integral in the normal direction. By using (5.19) and the surface divergence theorem, we obtain

where the energy spent by the interface deformation and the effects of the interface curvature are taken into account in our jump condition for the energy balance at the interface. Equation (5.37) agrees with the result obtained by Andrea (Reference Andrea2011), where the energy balance condition at the interface is derived by using a pillbox for the sharp-interface model. Again if we assume that there is no phase change across the interface, (5.37) then reduces to

If we further ignore the energy spent by the interface deformation and the effects of interface curvature, we can obtain the classical jump condition for the energy equation,

which is widely used for the computations of the sharp-interface model (e.g. Tavener & Cliffe Reference Tavener and Cliffe2002).

6.1. Simplified model and the weak form

For the sake of simplicity, we assume that the densities of the two fluids are matched. The system (4.26)–(4.31) can then be simplified in the form

where the variable thermal conductivity (3.44) is employed. Here, we employed the energy balance equation (3.19) instead of (4.28). The reason is that in the weak formulation of (4.28), the second-order derivative is involved, implying that more reductive $C^{1}$ finite elements are needed for the conformity. However, in the weak formulation of (6.3) ((6.8) below) we find that only first-order derivatives of $c$ are involved, so that the $C^{0}$ finite element method may be used for our computations. The benefits of using $C^{0}$ elements are obvious, in that the method can have more choices of elements and many existing codes can be incorporated to reduce various complications. It should be noted that (6.1)–(6.3) of the system will be computed for the example of thermocapillary convection, (6.1)–(6.5) will be computed for the example of thermocapillary migration with zero Marangoni number, and the non-dimensional system (4.50)–(4.54) will be computed for the example of thermocapillary migration with finite Marangoni number. For simplicity, we only present the numerical scheme for the dimensional system (6.1)–(6.5). The numerical method for the non-dimensional system (4.50)–(4.54) can be obtained correspondingly. By multiplying the system (6.1)–(6.5) by the test functions $q$ , $\boldsymbol{u}$ , ${\it\chi}$ , ${\it\phi}$ and ${\it\psi}$ respectively and using integration by parts, the weak form can be derived straightforwardly (where $\boldsymbol{v}$ , $p$ , $\hat{u}$ , $c$ , ${\it\mu}$ and the test functions $\boldsymbol{u}$ , $q$ , ${\it\chi}$ , ${\it\phi}$ and ${\it\psi}$ are in appropriate spaces),

6.2. Temporal schemes and implementation issues

The solution of the weak form (6.6)–(6.10) is approximated by a finite difference scheme in time and a conformal $C^{0}$ finite element method in space. To ensure the stability of our numerical method, we adopt the fully implicit backward Euler scheme to compute the problem.

We let ${\rm\Delta}t>0$ represent a time step size, and ( $\boldsymbol{v}_{h}^{n},p_{h}^{n},\hat{u} _{h}^{n},c_{h}^{n},{{\it\mu}_{C}}_{h}^{n}$ ) (in a finite dimensional space given by a finite element discretization of the computational domain ${\it\Omega}$ ) is an approximation of $(\boldsymbol{v},p,\hat{u} ,c,{\it\mu})$ at time $t^{n}=n{\rm\Delta}t$ , where $\boldsymbol{v}_{h}^{n}=\boldsymbol{v}(n{\rm\Delta}t)$ , $p_{h}^{n}=p(n{\rm\Delta}t)$ , $\hat{u} _{h}^{n}=\hat{u} (n{\rm\Delta}t)$ , $c_{h}^{n}=c(n{\rm\Delta}t)$ and ${{\it\mu}_{C}}_{h}^{n}={\it\mu}_{C}(n{\rm\Delta}t)$ . Then the approximation at time $t^{n+1}$ is denoted as ( $\boldsymbol{v}_{h}^{n+1},p_{h}^{n+1},\hat{u} _{h}^{n+1},c_{h}^{n+1},{{\it\mu}_{C}}_{h}^{n+1}$ ) and computed by the following finite element scheme:

where $\boldsymbol{v}_{\bar{t}}^{n+1}=(\boldsymbol{v}_{h}^{n+1}-\boldsymbol{v}_{h}^{n})/{\rm\Delta}t$ , $\hat{u} _{\bar{t}}^{n+1}=(\hat{u} _{h}^{n+1}-\hat{u} _{h}^{n})/{\rm\Delta}t$ and $c_{\bar{t}}^{n+1}=(c_{h}^{n+1}-c_{h}^{n})/{\rm\Delta}t$ . It should be noted that the divergence-free equation needs to be treated carefully in incompressible flow computations. Here, we rewrite (6.11) in the penalty formulation, where ${\it\delta}$ is a relatively small parameter and is set to be ${\it\delta}=10^{-6}$ for all the computations. It should be noted that for every time step, $T^{n+1}$ can be obtained by using (4.43), such that

Since the scheme is nonlinearly implicit we need to perform the linearization and then solve a linear system iteratively at each time step. We follow the numerical methods designed by Hua et al. (Reference Hua, Lin, Liu and Wang2011), where the linear system is symmetric and does not depend on time. Therefore, we only need to perform the Cholesky factorization for the symmetric linear system at the initial time step. After the initial time we do not need to factorize the linear system again since the coefficient matrix is independent of time.

For a phase-field model, it is sufficient to finely resolve only the interfacial region, and a fixed grid meshing represents a waste of computational resources. Therefore, an efficient adapting mesh that resolves the thin interfacial region is desirable. For the examples of thermocapillary convection, we design a mesh that has relatively high-resolution grids near the flat interface. For the example of thermocapillary migration, since the interface moves as the drop rises, an adaptive mesh is designed, in which there is a smaller frame that moves with the drop. Within the frame, the resolution of grids is much higher than those outside the moving frame, so that the moving interface of the drop can be resolved purposely. Here, only the meshes for the example of thermocapillary migration are shown.

Figure 3. Schematic diagram showing two immiscible fluids in a microchannel. The temperatures of the lower and upper plates are $T^{b}(x,-b)=T_{h}+T_{0}\cos (kx)$ and $T^{a}(x,a)=T_{c}$ respectively, where $T_{h}>T_{c}>T_{0}$ and $k=2{\rm\pi}/l$ is the wavenumber, and $a$ and $b$ are the heights of the fluids A and B respectively.

6.3. Thermocapillary convection in a two-layer fluid system

We now investigate the thermocapillary convection in a heated microchannel with two-layer superimposed fluids with a planar interface (Pendse & Esmaeeli Reference Pendse and Esmaeeli2010). We consider two-layer fluids (figure 3), where the heights of fluid A (upper) and fluid B (lower) are $a$ and $b$ respectively, and the fluids are of infinite extension in the horizontal direction. The physical properties of the fluids are their densities, viscosities and heat conductivities. The temperature variations in the present study are considered to be small enough so that the thermophysical properties of each fluid are assumed to remain constant, with the exception of surface tension. The temperatures of the lower and upper plates are

respectively, where $T_{h}>T_{c}>T_{0}>0$ , and ${\it\omega}=2{\rm\pi}/l$ is a wavenumber with $l$ being the channel length. The above temperature boundary conditions establish a temperature field that is periodic in the horizontal direction with a period of $l$ . Therefore, it is sufficient to only focus on the solution in one period, i.e. $-l/2<x<l/2$ . In the limit of zero Marangoni number and small Reynolds number, it is possible to ignore the convective transport of momentum and energy. In addition, we assume that the interface is to remain flat. By solving the simplified sharp-interface governing equations with the corresponding jump boundary conditions at the interface, Pendse & Esmaeeli (Reference Pendse and Esmaeeli2010) obtained the analytical solutions for the temperature field $\bar{T}(x,y)$ and stream function $\bar{{\it\psi}}(x,y)$ , where for the upper fluid

and for the lower fluid

In the above equations the unknowns are defined by $\tilde{k}=k_{A}/k_{B}$ , ${\it\alpha}=a{\it\omega}$ , ${\it\beta}=b{\it\omega}$ , $f({\it\alpha},{\it\beta},\tilde{k})=1/(\tilde{k}\sinh ({\it\beta})\cosh {\it\alpha}+\sinh ({\it\alpha})\cosh {\it\beta})$ , $g({\it\alpha},{\it\beta},\tilde{k})=\sinh ({\it\alpha})f({\it\alpha},{\it\beta},\tilde{k})$ and Based on their work, the simulations for our phase-field model are carried out in a 2D domain $[-l/2,l/2]\times [-b,a]$ with $l=1.6\times 10^{-4}$ and $a=b=4\times 10^{-5}$ . As the interface between the two fluids is assumed to be flat and rigid, the initial conditions for the phase variable only depend on $y$ , and can be given in the form

The periodic boundary conditions are applied on the left and right sides of the domain. On the top and bottom walls, the no-slip boundary conditions are imposed such that

Equations (6.17a,b ) are used as the boundary conditions for temperature with $T_{h}=20$ , $T_{c}=10$ and $T_{0}=4$ . We let the ratio parameter ${\it\eta}=6\sqrt{2}$  (5.29). Moreover, the fluid properties are shown in table 1. To show the influences of the thermal conductivity ratio on the stream-function and temperature fields, the simulations are carried out for two cases with different values of $\tilde{k}$ , where $\tilde{k}=0.1$ for case 1 and $\tilde{k}=0.5$ for case 2. Here, the variable thermal conductivity $k(c)$  (3.44) is employed, where we fix $k_{B}(=0.2)$ , and change the value of $k_{A}$ for the two cases. The contours of the temperature fields and stream functions for two cases at ${\it\epsilon}=0.002$ are shown in figures 4 and 5 respectively. It can be seen that our numerical results are in good agreement with the analytical solutions. In order to show that our phase-field model approaches the sharp-interface model as the thickness of the diffuse interface goes to zero, the computations are carried out by using five different values of ${\it\epsilon}$ ( $=0.02,0.01,0.005,0.002,0.001$ ). The $L^{2}$ norms of the relative differences between the numerical results and analytical solutions are shown in table 2. We can observe that as the value of ${\it\epsilon}$ decreases, the $L^{2}$ norm of the relative differences decreases for both the temperature field and the stream function. We also note that there are slight differences between our numerical results and the analytical predictions. The reason is twofold. First, and most importantly, the thickness of the interface of our model is finite, and the thermal diffusivity changes across it. Second, the viscous heating term is considered in our energy balance equation (6.3). As can be observed from the isotherms in figure 4, the cosine-like boundary condition for temperature leads to non-uniform distributions of the temperature along the interface. This results in a shear force along the interface that is from the centre to both sides of the domain. The fluids are set in motion by this shear force and move from the middle towards both sides of the domain. These are then replaced by the fluid flowing downwards (upward) from the top (bottom) boundary. Also, as the domain is periodic in the horizontal direction, the velocities of fluids that move towards both sides decrease and the fluids are forced to move upward (downward) to the top (bottom) of the domain. This mechanism results in the formation of the circulation patterns that can be observed in the stream-function fields (figure 5), where the fluid flow consists of four counter-rotating circulations that divide the domain into four parts. Moreover, in the context of the thermal conductivity ratio, we find that the decrease of $\tilde{k}$ leads to a more non-uniform distribution of temperature along the interface, and thus strengthens the thermocapillary convection. This result agrees with the recent result obtained by Liu et al. (Reference Liu, Valocchi, Zhang and Kang2014), where the same thermocapillary convection in a two-layer fluid system was investigated numerically by using a lattice Boltzmann phase-field method.

Figure 4. Isotherms of the numerical results and analytical solutions for the example of thermocapillary convection in a two-layer fluid system with different thermal diffusivity ratios $\tilde{k}=0.1$ and $\tilde{k}=0.5$ .

Figure 5. Streamlines of the numerical results and analytical solutions for the example of thermocapillary convection in a two-layer fluid system with different thermal diffusivity ratios $\tilde{k}=0.1$ and $\tilde{k}=0.5$ . Positive (negative) values of the stream function indicate the clockwise (the counterclockwise) circulation.

6.4. Thermocapillary migration in the limit of zero Marangoni number

The thermocapillary migration of a drop was first examined experimentally by Young et al. (Reference Young, Goldstein and Block1959), who derived an analytical expression for the terminal velocity (also known as the YGB velocity) of the drop in an infinite domain. In this study, both the Marangoni and Reynolds numbers are assumed to be infinitely small, such that the convective transport of momentum and energy is negligible. Instead, the terminal velocity of the drop is derived in an infinite domain with constant temperature gradient fields, and can be given in the form

where $U=-{\it\sigma}_{T}G_{T}R/{\it\mu}_{B}$ is chosen as the velocity scale, $R$ is the radius of the drop, $G_{T}$ stands for the constant temperature gradient, $\tilde{k}=k_{A}/k_{B}$ is the thermal conductivity ratio and $\tilde{{\it\mu}}={\it\mu}_{A}/{\it\mu}_{B}$ is the viscosity ratio between the two fluids. In our simulation, we consider a 2D domain ${\it\Omega}$ of size $[0,7.5R]\times [0,15R]$ where a planar 2D circular drop of fluid A with radius $R=0.1$ is placed inside the medium of fluid B, with the centre of the drop located at the centre of the box $(x_{c},y_{c})=(3.75R,7.5R)$ . We set the initial condition for the phase field as

In figure 6 we present the initial condition (6.27) for the whole domain (left-hand side), and for fixed $x=3.75R$ (right-hand side), where it can be observed that the area with $c=1$ represents the drop (fluid A) and the area with $c=0$ represents the medium (fluid B), between which the value of $c$ varies rapidly, resulting in a diffuse interface with finite thickness. Within this transition layer, the dotted contour line is at the level $c=0.5$ representing the dividing surface ${\it\Gamma}$ . No-slip boundary conditions are imposed on the top and bottom walls, and periodic boundary conditions are imposed in the horizontal direction. A linear temperature field is imposed in the $y$ direction,

with $T_{b}=10$ on the bottom wall and $T_{t}=25$ on the top wall, resulting in a constant temperature gradient, $G_{T}=10$ . Again, we let the ratio parameter ${\it\eta}=6\sqrt{2}$  (5.29). Moreover, the fluid properties are shown in table 3. Using these values, the theoretical terminal velocity of a spherical drop can be given as

Numerically, we use the following equation to calculate the rise velocity $v_{r}$ of the drop for our phase-field model:

where $\hat{\boldsymbol{j}}$ is the component of the unit vector in the $y$ direction.

Figure 7 shows the temporal evolution of the drop velocity normalized by $V_{YGB}$ between two different interface capturing methods, the phase-field method and level-set method (Herrmann et al. Reference Herrmann, Lopez, Brady, Raessi, Moin, Mansour and You2008). Similarly to the previous example in § 6.3, we compute our model by using two different interfacial thicknesses corresponding to ${\it\epsilon}=0.002$ and $0.001$ . Both the phase-field method and the level-set method seem to converge to a value of $v_{r}/V_{YGB}=0.8$ , roughly $20\,\%$ different from the theoretical prediction. The reason for this discrepancy is twofold. First, and most importantly, the theoretical rise velocity is for an axisymmetric sphere, whereas our simulations are carried out for a planar 2D drop. Second, the simulations include small blockage effects from the finite computational domain size as well as minute deformations of the drop, whereas the theoretical formula assumes an infinite domain and a non-deformable drop. As the thickness of the diffuse interface decreases, our results seem to coverage to the that obtained by the level-set method (Herrmann et al. Reference Herrmann, Lopez, Brady, Raessi, Moin, Mansour and You2008). For the case ${\it\epsilon}=0.001$ , we present the streamlines together with the moving interface at $t=1$ and $t=50$ in figure 8, where we observe that the streamlines for both cases exhibit similar patterns, with two asymmetric recirculations around the drop. Figure 8 shows the meshes together with the drop interface at $t=1$ and $t=50$ . Here, the size of the smaller frame is set to be $[3R\times 3R]$ , in which we take the shortest edge of the grids inside the frame as $15R/1000={\it\epsilon}$ , so that at least 7–9 grid cells (corresponding to the definition of the interfacial thickness) are located across the interface to ensure accuracy of our computations. In addition, the moving velocity of the frame is set to be equal to the drop rising velocity $v_{frame}=v_{r}$ , such that, through this relative long-term behaviour, the rising drop is always kept inside the smaller moving frame.

Figure 6. Initial condition of the phase variable $c$ for the example of the thermocapillary migration of a drop. The dotted line stands for the dividing surface.

Figure 7. The time evolution of the normalized migration velocity of a drop. The dashed lines are our numerical results for a 2D planar drop ( $v_{r}$ ), while the solid line represents the numerical results by using the level-set method.

Figure 8. The drop interface (black) and the streamlines (grey-scale lines, left), and the meshes (grey lines, right) at (a)  $t=1$ and (b)  $t=50$ . Positive values of the stream function indicate the clockwise circulation and negative values of the stream function indicate the counterclockwise circulation.

6.5. Thermocapillary migration with finite Marangoni number

We now compute the example of the thermocapillary motion of a drop with finite Marangoni number. Due to the finite Marangoni number, the energy equation (6.3) is coupled with the momentum equation (6.2). This is expected to result in a reduction of the tangential temperature gradients at the drop interface due to the interfacial flow driven by the Marangoni stress, which in turn will also be reduced. In this simulation, we consider a 2D domain ${\it\Omega}$ of size $[0,10R]\times [0,15R]$ , where a planar 2D circular drop of fluid A with radius $R=0.5$ is placed inside the medium of fluid B, with the centre of the drop located at the centre of the box $(x_{c},y_{c})=(0,1.5R)$ . At $t=0$ , (6.27) is employed as the initial condition for the phase variable, and a linear temperature distribution from $T_{b}=0$ at the bottom to $T_{t}=1$ at the top is imposed for the bulk liquid, and we assume that the drop has the same initial linear temperature distribution as the bulk liquid. Again, no-slip boundary conditions are imposed on the top and bottom boundaries, and periodic boundary conditions are imposed in the horizontal direction. The two fluids are assumed to have the same densities and viscosities. We set the thermal conductivity $k_{1}=0.1$ for the drop and $k_{2}=1$ for the bulk fluid. In this section, the non-dimensionalized system (4.50)–(4.54) is computed, where we set the non-dimensional parameters as ${\it\epsilon}=0.002$ , $Re=10$ , $M=1$ , $Pe=100/{\it\epsilon}$ , $Ca=1$ , $Ec=1$ . Five different values of the Marangoni number are employed for the computations, namely $Ma=50,100,500,1000,1500$ .

Figure 9 shows the velocity of the drop versus time for the five cases. As the time processes, the rise velocity reduces in all five cases, and we can observe that an increase in $Ma$ leads to a decrease in the rise velocity, which is consistent with the simulations by Herrmann et al. (Reference Herrmann, Lopez, Brady, Raessi, Moin, Mansour and You2008), Yin et al. (Reference Yin, Gao, Hu and Chang2008) and Zhao et al. (Reference Zhao, Li, Li and Li2010).

Figure 10 shows snapshots of the isotherms at four different times for the corresponding three cases, where the dependence of the migration velocity on the Marangoni number can be easily explained. Obviously, the enhanced convective transport of momentum and heat with increase of the Marangoni number results in more disturbances of the temperature field. Inside the drop, as we increase the Marangoni number, larger variations can be observed, leading to a substantial reduction in the surface temperature gradient and the corresponding rise velocities.

Figure 9. The time evolution of the rise velocity of a drop with different finite Marangoni numbers.

Figure 10. Snapshots of the drop interface (black) and isotherms (grey-scale lines) for different times and different values of $Ma$ as indicated.