%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1 Introduction
Bubbles and droplets are ubiquitous in natural processes and technological applications. Engineers, mathematicians and physicists have attempted to understand these phenomena over the past 200 years (Clift, Grace & Weber Reference Clift, Grace and Weber2005). Studies can be traced back to Young (Reference Young1805) and Laplace (Reference marquis de Laplace1805), who described the nonlinear partial differential equation governing the pressure jump that develops across a curved interface between fluids. This model is known as the Young–Laplace equation. Processes that drive the dynamics of bubbles and droplets include rainfall, air pollution, boiling, flotation, topological phase transition and fermentation (Clift et al. Reference Clift, Grace and Weber2005). Also, there has been intensive study of the coalescence of droplets (Kavehpour Reference Kavehpour2015) as well as droplet impact on solid surfaces (Josserand & Thoroddsen Reference Josserand and Thoroddsen2016).
Several methodologies have been developed to model pressure jumps and sharp gradients. Among them, sharp- and diffuse-interface methods are widely used to model phase segregation problems coupled with fluid dynamics. Sharp-interface models include the level-set (Sethian Reference Sethian1999) and volume-of-fluid (Rider & Kothe Reference Rider and Kothe1998) methods. Diffuse-interface models (DIM) have been used to model topological transitions in compressible (Anderson, McFadden & Wheeler Reference Anderson, McFadden and Wheeler1998) and quasi-incompressible (Lowengrub & Truskinovsky Reference Lowengrub and Truskinovsky1998) flows under consistent thermodynamic premises. A wide range of phenomena is described via diffuse interfaces, ranging from material sciences (Loginova, Amberg & Ågren Reference Loginova, Amberg and Ågren2001) and fracture mechanics (Spatschek et al. Reference Spatschek, Müller-Gugenberger, Brener and Nestler2007) to fluid dynamics (Khatavkar, Anderson & Meijer Reference Khatavkar, Anderson and Meijer2006, Reference Khatavkar, Anderson and Meijer2007; Gomez, Cueto-Felgueroso & Juanes Reference Gomez, Cueto-Felgueroso and Juanes2013). A diffuse-interface model circumvents several numerical difficulties (Yue et al. Reference Yue, Feng, Liu and Shen2004), and its inherent thermodynamic consistency allows the method to incorporate rheology when the phase-field evolution is described by a free energy density (Emmerich Reference Emmerich2003).
In this work we use the Cahn–Hilliard equation with the Ginzburg–Landau free energy defined in Gómez et al. (Reference Gómez, Calo, Bazilevs and Hughes2008). A derivation of the Cahn–Hilliard equation is presented by Gurtin (Reference Gurtin1996) using a balance law for microforces. Gurtin, Polignone & Viñals (Reference Gurtin, Polignone and Viñals1996) extended that work to couple phase segregation with hydrodynamics, providing the first derivation of the Navier–Stokes–Cahn–Hilliard system completely based on continuum mechanics and thermodynamics arguments. This theory uses the classical balance laws of mass and momentum, together with a new balance law for microforces. This microforce balance takes into account the ‘microscopic work’ done by changes in the order parameter, which represents, in the context of our work, the dimensionless concentration of a phase.
The introduction of microforces and microstresses, as well as their respective balance laws, modifies the mechanical version of the second law of thermodynamics (Gurtin et al.
Reference Gurtin, Polignone and Viñals1996, equation (30)), which in turn influences the constitutive relations under consideration. In particular, a byproduct of this framework is that the capillary stress tensor naturally appears as a constraint imposed by the mechanical version of the second law of thermodynamics on the constitutive equations (Gurtin et al.
Reference Gurtin, Polignone and Viñals1996, equation (48)). Liu (Reference Liu2014) extended the framework of Gurtin to encompass the
$n$
-component Navier–Stokes–Cahn–Hilliard (NSCH) multiphase system with a choice of switching to van der Waals theory.
In this work, we use the dimensionless form of the equations of motion to develop the complete energy budget equation of the coupled NSCH system to model the dynamics of liquid droplets in a liquid continuum. In addition, from a physical point of view, we give an interpretation of each quantity involved in the energy exchange and analyse the behaviour of the mass flux across the interface.
From a numerical point of view, we develop a general and robust formulation based on a finite-dimensional, high-order isogeometric analysis approximation. We use divergence-conforming B-spline spaces to obtain a discrete pointwise divergence-free velocity field (Evans & Hughes Reference Evans and Hughes2013a ,Reference Evans and Hughes b ,Reference Evans and Hughes c ; Sarmiento et al. Reference Sarmiento, Cortes, Garcia, Dalcin, Collier and Calo2015). We implement this discretization in PetIGA, a high-performance isogeometric analysis framework (Dalcin et al. Reference Dalcin, Collier, Vignal, Côrtes and Calo2015; Vignal et al. Reference Vignal, Dalcin, Brown, Collier and Calo2015a ). Two- and three-dimensional numerical results highlight the robustness of the framework. Finally, we compare our numerical results against analytical features of droplet dynamics to verify our model.
The outline of this work is as follows. Section 2 presents the general governing equations, followed by § 3, which describes the constitutive relations adopted here. Section 4 introduces the dimensionless set of equations together with the dimensionless constitutive equations. Section 5 details the overall mass flux behaviour characterized by the Cahn–Hilliard equation. In § 6, we present the energy budget of the NSCH equations under the Boussinesq premise. This concludes the theoretical analysis and preliminaries of the paper. In § 7 we describe the numerical scheme we use, whereas in § 8 the set-up of the numerical experiments is detailed. In § 9, we analyse the numerical results. Finally, in § 10, we briefly detail the conclusions. Appendix A presents the first and second laws of thermodynamics, showing that our model is thermodynamically consistent. Appendix B lists the identities obtained to develop the energy budget of the NSCH flows.
2 General governing equations
The dynamics of binary immiscible fluids involves both mass and momentum transfer. In this work, the phase segregation phenomenon is described by the Cahn–Hilliard equation (Cahn & Hilliard Reference Cahn and Hilliard1958, Reference Cahn and Hilliard1959; Cahn Reference Cahn1959), while the hydrodynamics is incorporated through the incompressible Navier–Stokes equations under the Boussinesq assumption to account for buoyancy effects (Meiburg & Kneller Reference Meiburg and Kneller2010; Espath et al. Reference Espath, Pinto, Laizet and Silvestrini2015b ). We assume a small density difference between the fluids. The resulting set of equations in dimensional form is given by
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\partial {\it\phi}}{\partial t}}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\phi}\boldsymbol{v})+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\left[{\displaystyle \frac{\partial \boldsymbol{v}}{\partial t}}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{v}\otimes \boldsymbol{v})\right]-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}-{\it\phi}\unicode[STIX]{x27E6}{\it\rho}\unicode[STIX]{x27E7}\boldsymbol{g}=\mathbf{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D64F}=\unicode[STIX]{x1D64F}^{\text{T}}, & \displaystyle\end{eqnarray}$$
${\it\rho}$
,
$\boldsymbol{v}$
,
$\unicode[STIX]{x1D64F}$
,
${\it\phi}\unicode[STIX]{x27E6}{\it\rho}\unicode[STIX]{x27E7}\boldsymbol{g}$
,
${\it\phi}$
and
$\boldsymbol{j}$
are the density, velocity, stress, buoyancy force, volume concentration and mass flux, respectively. The density is defined as the weighted mean of the fluid densities,
${\it\rho}={\it\rho}_{2}{\it\phi}+{\it\rho}_{1}(1-{\it\phi})$
, where the subscript indicates the corresponding fluid. We only consider the density excess
$\unicode[STIX]{x27E6}{\it\rho}\unicode[STIX]{x27E7}={\it\rho}_{2}-{\it\rho}_{1}$
in the buoyancy force. We denote first- and second-order tensors by bold lower and upper case symbols, respectively. The differential operators
$\boldsymbol{{\rm\nabla}}(\cdot )$
and
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{\cdot })$
represent the gradient and divergence, respectively. The superscript
$(\cdot )^{\text{T}}$
denotes the transpose and
$\otimes$
the tensor (dyadic) product.
3 Constitutive equations
3.1 Mass flux
According to Landau & Lifshitz (Reference Landau and Lifshitz1959), in the absence of heat transfer, the mass flux may be defined by generalizing Fick’s law,
$\boldsymbol{j}=-{\it\alpha}\,\boldsymbol{{\rm\nabla}}{\it\eta}$
, where
${\it\eta}$
is the chemical potential and
${\it\alpha}$
is the chemical mass diffusivity. Once the chemical potential is defined by the Ginzburg–Landau free energy density, i.e.
${\it\eta}\stackrel{\text{def}}{=}{\it\delta}{\it\Psi}/{\it\delta}{\it\phi}$
[energy/volume] (the variational derivative of
${\it\Psi}$
), the Cahn–Hilliard equation is obtained (Emmerich Reference Emmerich2003; Gómez et al.
Reference Gómez, Calo, Bazilevs and Hughes2008). Additionally, the Ginzburg–Landau free energy, in a volume
${\it\Omega}$
, may be written as
$$\begin{eqnarray}{\it\Psi}[{\it\phi}]=\int _{{\it\Omega}}{\it\psi}\,\text{d}{\it\Omega}=\int _{{\it\Omega}}({\it\psi}_{{\it\phi}}+{\it\psi}_{s})\,\text{d}{\it\Omega},\end{eqnarray}$$
where the bulk and interfacial free energy densities
${\it\psi}_{{\it\phi}}$
and
${\it\psi}_{s}$
, respectively, are defined as (Cahn & Hilliard Reference Cahn and Hilliard1958, equations (3.1) and (3.13))
$$\begin{eqnarray}\displaystyle & {\it\psi}_{{\it\phi}}=N_{v}k_{B}{\it\theta}({\it\phi}\ln {\it\phi}+(1-{\it\phi})\ln (1-{\it\phi}))+N_{v}{\it\omega}(1-{\it\phi}){\it\phi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\psi}_{s}={\displaystyle \frac{{\it\gamma}_{{\it\phi}}}{2}}\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}. & \displaystyle\end{eqnarray}$$
$N_{v}$
is the number of molecules per unit volume,
$k_{B}$
is Boltzmann’s constant and
${\it\omega}$
is an interaction energy given by
${\it\omega}=2k_{B}{\it\theta}_{c}$
. The interaction energy is positive and is related to the critical temperature,
${\it\theta}_{c}$
. In the interfacial free energy density term,
${\it\gamma}_{{\it\phi}}={\it\sigma}l\,[\text{force}]$
represents the magnitude of the interfacial energy. The parameters
${\it\sigma}$
and
$l$
are the interface tension
$[\text{force}/\text{length}]$
and the interfacial thickness
$[\text{length}]$
, respectively. Alternatively,
$l$
is a root mean square effective ‘interaction distance’, according to the original work of Cahn & Hilliard (Reference Cahn and Hilliard1959). This force term
${\it\gamma}_{{\it\phi}}$
is defined in Cahn & Hilliard (Reference Cahn and Hilliard1958) (equation (3.12)) and Cahn & Hilliard (Reference Cahn and Hilliard1959) (equation (4.1)) by
${\it\gamma}_{{\it\phi}}=N_{v}{\it\omega}l^{2}$
.
Without loss of generality, we define the interface thickness to have the simplest expression for
${\it\gamma}_{{\it\phi}}$
, i.e.
${\it\sigma}l$
. Yue et al. (Reference Yue, Feng, Liu and Shen2004) opted to define this relation as
${\it\gamma}_{{\it\phi}}=3{\it\sigma}l/(2\sqrt{2})$
, based on the quartic bulk free energy density when the diffuse interface is at equilibrium (
${\it\delta}{\it\Psi}/{\it\delta}{\it\phi}=0$
).
The scalar parameter
${\it\alpha}\;[\text{time}\times \text{volume}/\text{mass}]$
in the mass flux is rewritten as
${\it\alpha}=M({\it\phi}){\it\beta}$
, where
$M({\it\phi})=M_{o}(1-{\it\phi}){\it\phi}\,[\text{length}^{2}/\text{time}]$
represents the mobility (with
$M_{o}$
a positive constant) and
${\it\beta}^{-1}={\it\rho}_{1}u_{m}^{2}\,([\text{mass}/\text{volume}]\times [\text{length}/\text{time}]^{2})$
is related to the kinetic energy at a molecular scale,
${\it\rho}_{1}u_{m}^{2}$
. The mobility is degenerate, i.e. the phase dependence of the mobility confines the molecular movement to the interface region. This effect is detailed in § 5. The mass flux can then be defined as
$$\begin{eqnarray}\boldsymbol{j}=-M({\it\phi}){\it\beta}\boldsymbol{{\rm\nabla}}{\it\eta}.\end{eqnarray}$$
Considering the buoyancy effects, we must include the potential energy in the free energy, i.e.
$$\begin{eqnarray}{\it\Psi}[{\it\phi}]=\int _{{\it\Omega}}{\it\psi}\,\text{d}{\it\Omega}=\int _{{\it\Omega}}({\it\psi}_{{\it\phi}}+{\it\psi}_{s}+e_{p})\,\text{d}{\it\Omega},\end{eqnarray}$$
where
$e_{p}={\it\phi}\unicode[STIX]{x27E6}{\it\rho}\unicode[STIX]{x27E7}gx_{2}$
is the potential energy,
$g$
is the gravity acceleration and
$x_{2}$
is the vertical coordinate. Thus, the total mass flux is now defined as
$$\begin{eqnarray}\boldsymbol{j}=-M({\it\phi}){\it\beta}\,\boldsymbol{{\rm\nabla}}({\it\eta}_{{\it\phi}}+{\it\eta}_{s})-M_{o}{\it\beta}\boldsymbol{{\rm\nabla}}{\it\eta}_{p},\end{eqnarray}$$
where
${\it\eta}_{{\it\phi}}$
,
${\it\eta}_{s}$
and
${\it\eta}_{p}$
are the chemical potentials related to the bulk, interfacial and potential energies, respectively. The mass flux due to the potential energy is defined with a constant mobility
$M_{o}$
, since this is not an interfacial but a bulk mass flux. Thus, for this term the chemical mass diffusivity
${\it\alpha}$
is constant.
3.2 Stress
For an incompressible Newtonian fluid, the viscous stress is given by Gurtin, Fried & Anand (Reference Gurtin, Fried and Anand2010) as
$$\begin{eqnarray}\unicode[STIX]{x1D64F}^{visc}=2{\it\mu}\unicode[STIX]{x1D63F},\end{eqnarray}$$
with
$$\begin{eqnarray}\unicode[STIX]{x1D63F}={\textstyle \frac{1}{2}}((\boldsymbol{{\rm\nabla}}\boldsymbol{v})^{\text{T}}+\boldsymbol{{\rm\nabla}}\boldsymbol{v}),\end{eqnarray}$$
where the dynamic viscosity,
${\it\mu}\stackrel{\text{def}}{=}{\it\mu}({\it\phi})$
relates the strain rate
$\unicode[STIX]{x1D63F}$
to the viscous stress tensor
$\unicode[STIX]{x1D64F}^{visc}$
. The Cauchy stress tensor is defined as
$\unicode[STIX]{x1D64F}^{c}=\unicode[STIX]{x1D64F}^{visc}-p\mathbf{1}$
(Gurtin et al.
Reference Gurtin, Fried and Anand2010), with the last term associating the pressure to the Lagrange multiplier that enforces the incompressibility constraint.
Considering a binary mixture where each phase is endowed with a different viscosity, we assume a smooth transition through the interface which is given by
${\it\mu}({\it\phi})={\it\mu}_{1}\text{e}^{m{\it\phi}}$
where
$m=\ln ({\it\mu}_{2}/{\it\mu}_{1})$
, with
${\it\mu}_{1}$
and
${\it\mu}_{2}$
being the viscosity of each phase. Here, we assume
${\it\mu}_{1}<{\it\mu}_{2}$
. To simplify the exposition in what follows, we define
$c=\text{e}^{m{\it\phi}}$
.
Using the balance of microforces, Gurtin et al. (Reference Gurtin, Polignone and Viñals1996) derived a thermodynamically consistent complement to the constitutive relation for the Cauchy stress via the mechanical version of the second law of thermodynamics, which models capillarity effects. The interfacial capillary stress assumes the form
$$\begin{eqnarray}\unicode[STIX]{x1D64F}^{s}=-{\it\gamma}_{{\it\phi}}\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}.\end{eqnarray}$$
Finally, the total stress associated with the macroscopic motion of the fluid is given by
$$\begin{eqnarray}\unicode[STIX]{x1D64F}=\unicode[STIX]{x1D64F}^{c}+\unicode[STIX]{x1D64F}^{s}=2{\it\mu}\unicode[STIX]{x1D63F}-p\mathbf{1}-{\it\gamma}_{{\it\phi}}\,\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi},\end{eqnarray}$$
where the stress has units of
$[\text{force}/\text{length}^{2}]$
.
Finally, our model differs from Gurtin’s model in the sense that we account for the potential energy in the free energy, and we assume that the density difference is small. This assumption justifies the use of the Boussinesq assumption in the momentum equation to account for buoyancy effects. In addition, we employ a logarithmic function to represent the free energy rather than the more often found polynomial approximation. In this sense, our work also differs from that of Jacqmin (Reference Jacqmin1999) and Jamet et al. (Reference Jamet, Lebaigue, Coutris and Delhaye2001), where the double-well polynomial function is used to describe the free energy, and where the potential energy is not included in the free energy description.
4 Governing equations in dimensionless form
To make the governing equations dimensionless, a length
$b$
and a velocity
$u$
are chosen as the characteristic length and velocity scales, respectively. The reference length and velocity are usually chosen as the droplet diameter or radius and the terminal velocity or buoyancy velocity, respectively. We scale the viscosities by the smallest viscosity
${\it\mu}_{1}$
, the pressure by
${\it\mu}_{1}u/b$
, the time by
$b/u$
and the mobility by
$M_{o}$
. The chemical potential is rendered dimensionless by
${\it\psi}_{c}=2N_{v}k_{B}{\it\theta}_{c}=N_{v}{\it\omega}$
, which is the critical free energy density. The phase field is inherently a dimensionless quantity and is normalized between
$(0,1)$
according to the bulk free energy density. The mixture law assumes the density to be of the form
${\it\rho}={\it\rho}_{1}+{\it\phi}\unicode[STIX]{x27E6}{\it\rho}\unicode[STIX]{x27E7}$
, with
$\unicode[STIX]{x27E6}{\it\rho}\unicode[STIX]{x27E7}={\it\rho}_{2}-{\it\rho}_{1}$
. In our examples,
${\it\rho}_{1}$
is the lighter fluid; nevertheless, the association between lower viscosity and lower density is incidental and can be reversed if the modelling requires this. In addition, due to the Boussinesq assumption we have that
$O({\it\phi}\unicode[STIX]{x27E6}{\it\rho}\unicode[STIX]{x27E7}/{\it\rho}_{1})\ll 1$
; usually of the order of 5 % would be admissible. Table 1 lists the dimensionless numbers we obtain with these scalings.
Dimensionless groups.
To the best of our knowledge, two new dimensionless groups appear in this physical interpretation. We denote these new numbers by
$Lm$
and
$Ln$
;
$Lm$
is the ratio between the critical free energy density and the capillary energy, i.e. the interfacial curvature energy. The capillary energy is closely related to the Laplace pressure obtained from the Young–Laplace equation (Myers Reference Myers1990). This pressure difference (or Laplace pressure)
$\unicode[STIX]{x27E6}p\unicode[STIX]{x27E7}={\it\sigma}(R_{1}^{-1}+R_{2}^{-1})$
relates the surface tension to the principal radii of curvature. In the particular case of a spherical droplet,
$\unicode[STIX]{x27E6}p\unicode[STIX]{x27E7}=2{\it\sigma}R^{-1}$
, where
$R$
is the radius of the sphere. The interfacial curvature energy takes the form
${\it\sigma}/b=\unicode[STIX]{x27E6}p\unicode[STIX]{x27E7}\!/4$
. The number
$Ln$
is the ratio between the kinetic energy at a molecular scale and the critical free energy density. In the particular case of an ideal gas,
$Ln$
can be interpreted by defining the molecular velocity
$u_{m}$
as
$u_{m}^{2}=u_{rms}^{2}=3k_{B}{\it\theta}/m$
, where
$u_{rms}$
is the root-mean-square speed of a single molecule of mass
$m$
. This assumption yields
$Ln={\it\rho}_{1}u_{rms}^{2}/{\it\psi}_{c}=3/2{\it\vartheta}^{-1}$
, where
${\it\vartheta}={\it\theta}_{c}/{\it\theta}$
defines the ratio between the critical and the absolute temperatures. Therefore,
$Ln$
can measure the deviation of the absolute temperature from the critical one in an ideal gas.
Henceforth, all quantities considered are dimensionless according to table 1. The dimensionless forms of the free energy densities are
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\psi}_{{\it\phi}}=\frac{1}{2{\it\vartheta}}({\it\phi}\ln {\it\phi}+(1-{\it\phi})\ln (1-{\it\phi}))+(1-{\it\phi}){\it\phi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\psi}_{s}=\frac{Cn}{2Lm}\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle e_{p}=\frac{Bo}{We}{\it\phi}x_{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}{\it\eta}={\displaystyle \frac{{\it\delta}{\it\Psi}}{{\it\delta}{\it\phi}}}={\displaystyle \frac{\partial {\it\psi}}{\partial {\it\phi}}}-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\displaystyle \frac{\partial {\it\psi}}{\partial \boldsymbol{{\rm\nabla}}{\it\phi}}}={\displaystyle \frac{1}{2{\it\vartheta}}}\ln {\displaystyle \frac{{\it\phi}}{1-{\it\phi}}}+1-2{\it\phi}-{\displaystyle \frac{Cn}{Lm}}{\rm\Delta}{\it\phi}+{\displaystyle \frac{Bo}{We}}x_{2}.\end{eqnarray}$$
The set of dimensionless differential equations is given by
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\partial {\it\phi}}{\partial t}}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\partial \boldsymbol{v}}{\partial t}}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{v}-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}+{\displaystyle \frac{Bo}{We}}{\it\phi}\boldsymbol{e}_{2}=\mathbf{0}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\boldsymbol{j}=-{\displaystyle \frac{1}{PeLn}}\left\{M({\it\phi})\,\boldsymbol{{\rm\nabla}}\left({\displaystyle \frac{1}{2{\it\vartheta}}}\ln {\displaystyle \frac{{\it\phi}}{1-{\it\phi}}}+1-2{\it\phi}-{\displaystyle \frac{Cn}{Lm}}{\rm\Delta}{\it\phi}\right)+{\displaystyle \frac{Bo}{We}}\boldsymbol{e}_{2}\right\},\end{eqnarray}$$
the stress is given by
$$\begin{eqnarray}\unicode[STIX]{x1D64F}={\displaystyle \frac{2c}{Re}}\unicode[STIX]{x1D63F}-{\displaystyle \frac{1}{Re}}p\mathbf{1}-{\displaystyle \frac{Cn}{We}}\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi},\end{eqnarray}$$
and
$c$
scales the kinetic viscosities of the different fluids. The rescaled mobility is equal to
$M({\it\phi})={\it\phi}(1-{\it\phi})$
. For convenience, we assume gravity points in the
$-\boldsymbol{e}_{2}$
direction.
To understand the components of the mass flux, we split the chemical potential into
${\it\eta}={\it\eta}_{{\it\phi}}+{\it\eta}_{s}+{\it\eta}_{p}$
where
${\it\eta}_{{\it\phi}}=1/2{\it\vartheta}\ln ({\it\phi}/(1-{\it\phi}))+1-2{\it\phi}$
,
${\it\eta}_{s}=-(Cn/Lm)\,{\rm\Delta}{\it\phi}$
and
${\it\eta}_{p}=(Bo/We)\,x_{2}$
, which yield the bulk, interfacial and gravitational mass fluxes, respectively.
5 Mass flux behaviour
In the Cahn–Hilliard equation, the mass flux is a phobic phenomenon that segregates the phases. The mass flux exhibits an ‘anomalous’ behaviour across the interface thickness, as it may change its sign (twice) and deviate from the gradient direction of the phase field. In the following, we present the analysis of the orientation of these mass fluxes. We start by analysing the flux structure in a one-dimensional example. Next, we describe some of the restrictive conditions under which the two mass flux terms are parallel. We conclude this section with a detailed description of a general two-dimensional simulation. To later explain the energy exchanges in the free energy densities, we first develop a good understanding of the mass flux behaviour. To do so, we rewrite the mass flux (4.4) as
$$\begin{eqnarray}\displaystyle \boldsymbol{j} & = & \displaystyle -{\displaystyle \frac{1}{PeLn}}M({\it\phi})\boldsymbol{{\rm\nabla}}({\it\eta}_{{\it\phi}}+{\it\eta}_{s})-{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Bo}{We}}\boldsymbol{{\rm\nabla}}{\it\eta}_{p}\nonumber\\ \displaystyle & = & \displaystyle -{\displaystyle \frac{M({\it\phi})}{PeLn}}\left[\left(-2+{\displaystyle \frac{1}{2{\it\vartheta}{\it\phi}(1-{\it\phi})}}\right)\boldsymbol{{\rm\nabla}}{\it\phi}-{\displaystyle \frac{Cn}{Lm}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\right]-{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Bo}{We}}\boldsymbol{e}_{2}.\end{eqnarray}$$
To simplify notation, we define the bulk,
$\boldsymbol{j}_{{\it\phi}}$
, interfacial,
$\boldsymbol{j}_{s}$
, and gravitational,
$\boldsymbol{j}_{p}$
, mass fluxes as
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}_{{\it\phi}}=-{\displaystyle \frac{1}{PeLn}}\left(-2+{\displaystyle \frac{1}{2{\it\vartheta}{\it\phi}(1-{\it\phi})}}\right)M({\it\phi})\boldsymbol{{\rm\nabla}}{\it\phi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}_{s}={\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Cn}{Lm}}M({\it\phi})\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}_{p}=-{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Bo}{We}}\boldsymbol{e}_{2}, & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D643}$
denotes the Hessian of
${\it\phi}$
. In general,
$\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}$
and
$\boldsymbol{{\rm\nabla}}{\it\eta}_{s}$
are not parallel, so the mass flux might not be normal to the interface implicitly defined by isovalues of the phase field.
First, let us observe how the mass flux may change sign across the interface. Figure 1 shows, in panel (a), the phase field and its derivatives up to third order in a one-dimensional problem. In this figure, lower values of
${\it\phi}$
represent the droplet (inclusion). Inward and outward are with respect to this region. In addition, figure 1 shows the mass flux behaviour assuming constant and degenerate mobilities (panels (b) and (c), respectively). In each panel the first row shows the phase field. The second and third rows in panels (b) and (c) show the dependence of the bulk mass flux
$\boldsymbol{j}_{{\it\phi}}$
on the temperature ratio
${\it\vartheta}$
. Finally, the fourth row depicts the interfacial mass flux
$\boldsymbol{j}_{s}$
. Red indicates the regions where the mass flux points inward, whereas blue indicates the regions where the mass flux points outward. The bulk mass flux,
$\boldsymbol{j}_{{\it\phi}}$
, points in the steepest descent direction if the temperature ratio is lower than unity,
${\it\vartheta}<1$
, while this flux changes its sign twice if the temperature ratio is greater than one,
${\it\vartheta}>1$
. On the other hand, the interfacial mass flux
$\boldsymbol{j}_{s}$
always exhibits two changes of sign. Thus, the total mass flux
$\boldsymbol{j}$
may have two changes in sign. These changes in the mass flux direction explain the phobic interactions between the phases. The function
$-2+1/(2{\it\vartheta}{\it\phi}(1-{\it\phi}))$
is positive definite if
${\it\vartheta}<1$
, otherwise if
${\it\vartheta}>1$
it becomes indefinite, with
${\it\vartheta}$
a finite value. Finally, if
${\it\vartheta}\rightarrow \infty$
the function is negative definite.
(a) Phase field and its spatial derivatives in one spatial dimension. (b) Mass fluxes with constant mobility in the cases
${\it\vartheta}<1$
and
${\it\vartheta}>1$
(
${\it\vartheta}={\it\theta}_{c}/{\it\theta}$
). (c) The same mass fluxes with degenerate mobility.
Remark 1 (The overall behaviour of the mass fluxes).
The overall behaviour of
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}$
is such that it has a positive inner product with
$-\boldsymbol{{\rm\nabla}}{\it\phi}$
when Neumann-free boundary conditions are considered. For
$\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{n}=0|_{{\it\Gamma}}$
, with
$\boldsymbol{n}$
the outward unit normal vector to the boundary
${\it\Gamma}$
, we obtain
$\int _{{\it\Omega}}-\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\,\text{d}{\it\Omega}=\int _{{\it\Omega}}({\rm\Delta}{\it\phi})^{2}\,\text{d}{\it\Omega}$
.
Remark 2 (Mass flux in two dimensions).
We consider an idealized two-dimensional case to understand how the interfacial mass flux
$\boldsymbol{j}_{s}$
may deviate the total mass flux from the steepest descent direction of the phase field. Let us define a local coordinate system given by the Frenet basis
$(\boldsymbol{e}_{\bot },\boldsymbol{e}_{\Vert })$
(Kühnel Reference Kühnel2006). The basis is defined over isovalue curves of
${\it\phi}$
, where
$\boldsymbol{e}_{\bot }$
is the unit vector in the direction of the gradient of
${\it\phi}$
,
$n$
, and
$\boldsymbol{e}_{\Vert }$
is the unit vector in the tangential direction,
${\it\tau}$
. We describe the coordinate transformation between the Frenet and Cartesian bases by taking the vector
$\text{d}\boldsymbol{x}=(\text{d}x_{1},\,\text{d}x_{2})$
and establishing that
$$\begin{eqnarray}\text{d}\boldsymbol{x}=\partial _{n}^{}\boldsymbol{x}\,\text{d}n+\partial _{{\it\tau}}^{}\boldsymbol{x}\,\text{d}{\it\tau},\end{eqnarray}$$
where
$\partial _{n}^{}=\partial /\partial n$
and
$\partial _{{\it\tau}}^{}=\partial /\partial {\it\tau}$
. We set the metric coefficients as the normal
$h_{n}\stackrel{\text{def}}{=}|\partial _{n}^{}\boldsymbol{x}|$
and tangential
$h_{{\it\tau}}\stackrel{\text{def}}{=}|\partial _{{\it\tau}}^{}\boldsymbol{x}|$
components. Thus, the gradient operator, the gradient of
${\it\phi}$
, the Laplacian operator and the Laplacian of
${\it\phi}$
are defined in this coordinate system as
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}(\cdot )=h_{n}^{-1}\partial _{n}^{}(\cdot )\boldsymbol{e}_{\bot }+h_{{\it\tau}}^{-1}\partial _{{\it\tau}}^{}(\cdot )\boldsymbol{e}_{\Vert }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}{\it\phi}=h_{n}^{-1}\partial _{n}^{}{\it\phi}\boldsymbol{e}_{\bot }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\rm\Delta}(\cdot )=(h_{n}h_{{\it\tau}})^{-1}[\partial _{n}^{}(h_{{\it\tau}}h_{n}^{-1}\partial _{n}^{}(\cdot ))+\partial _{{\it\tau}}^{}(h_{{\it\tau}}^{-1}h_{n}\partial _{{\it\tau}}^{}(\cdot ))], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\rm\Delta}{\it\phi}=(h_{n}h_{{\it\tau}})^{-1}\partial _{n}^{}(h_{{\it\tau}}h_{n}^{-1}\partial _{n}^{}{\it\phi}). & \displaystyle\end{eqnarray}$$
${\it\phi}$
in this orthogonal curvilinear coordinate system is 
$$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643} & = & \displaystyle h_{n}^{-1}\{\partial _{n}^{}(h_{n}h_{{\it\tau}})^{-1}[\partial _{n}^{}(h_{n}^{-1}h_{{\it\tau}}\partial _{n}^{}{\it\phi})]+(h_{n}h_{{\it\tau}})^{-1}[\partial _{n}^{2}(h_{n}^{-1}h_{{\it\tau}}\partial _{n}^{}{\it\phi})]\}\boldsymbol{e}_{\bot }\nonumber\\ \displaystyle & & \displaystyle +\,h_{{\it\tau}}^{-1}\{\partial _{{\it\tau}}^{}(h_{n}h_{{\it\tau}})^{-1}[\partial _{n}^{}(h_{n}^{-1}h_{{\it\tau}}\partial _{n}^{}{\it\phi})]+(h_{n}h_{{\it\tau}})^{-1}[\partial _{{\it\tau}}^{}\partial _{n}^{}(h_{n}^{-1}h_{{\it\tau}}\partial _{n}^{}{\it\phi})]\}\boldsymbol{e}_{\Vert }.\end{eqnarray}$$
By assuming that isocurves of
${\it\phi}$
are parallel (i.e. constant interface thickness) and the tangential metric along
${\it\tau}$
is constant (i.e. the isocurves of
${\it\phi}$
are circumferences or straight lines), we conclude that (
$\partial _{n}^{}\boldsymbol{e}_{\bot }=\partial _{n}^{}\boldsymbol{e}_{\Vert }=\mathbf{0}$
$\Rightarrow$
$\partial _{n}^{}h_{n}=\partial _{{\it\tau}}^{}h_{n}=0$
) and (
$\partial _{{\it\tau}}^{}h_{{\it\tau}}=0$
), respectively. Thus, under these premises, we write the divergence of the Hessian of
${\it\phi}$
as
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}={\it\varphi}\boldsymbol{e}_{\bot }\end{eqnarray}$$
with
${\it\varphi}=h_{n}^{-1}\partial _{n}^{}h_{{\it\tau}}[\partial _{n}^{}h_{{\it\tau}}\partial _{n}^{}{\it\phi}+h_{{\it\tau}}\partial _{n}^{2}{\it\phi}+(h_{n}^{2}h_{{\it\tau}})^{-1}\partial _{n}^{2}h_{{\it\tau}}\partial _{n}^{}{\it\phi}+h_{n}^{-2}\partial _{n}^{3}{\it\phi}]$
.
The idealized example explained above deals with isocontours of
${\it\phi}$
endowed with parallel and constant curvature to show the simplest behaviour of the mass flux, yielding a mass flux normal to the interfaces. In the general case, the mass flux deviates from the normal direction due to changes in the metric coefficients in both directions.
Remark 3 (Mass flux in steady state solutions).
From Remark 2, we conclude that the steady state solution of the Cahn–Hilliard equation yields
$\boldsymbol{{\rm\nabla}}{\it\phi}\Vert \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}$
. We stress that this conclusion does not consider hydrodynamic effects.
Figure 2 shows three stages in the merging of two rising droplets and how the mass flux behaves during this process. At an early stage, details 1–3 depict the mass flux on three isocontours of the phase field,
${\it\phi}=0.1$
, 0.5 and 0.9. Across the interface, in the leading zone of the top droplet, the mass flux points inward and is almost normal to the interface; cf. detail 1. However, detail 2 shows a gradual change in the mass flux direction. In this region, the mass flux points outward on a isocontour
${\it\phi}=0.5$
, whereas on a isocontour
${\it\phi}=0.9$
the mass flux points inward. In the regions where the phase field isocontours exhibit a large change in curvature, which happens when droplets merge, we observe a significant deviation of the mass flux from the steepest descent direction, as detail 3 depicts. In the trailing zone, the mass flux points outward during the entire merging process. The largest values of the tangential component of the mass flux are located in the highly curved region of the phase field. At the last stage, instant
$t_{3}$
, large mass flux deviations are not observed.
Mass flux over isocontours of
${\it\phi}=0.1$
, 0.5 and 0.9 at three instants
$t_{1}$
,
$t_{2}$
and
$t_{3}$
when two droplets are merging and rising. Three details at
$t_{1}$
show regions where the mass flux changes its sign and deviates from the normal direction.
In figure 3(a) the red (blue) region represents where the mass flux points inward (outward) with respect to the droplets. In the leading (trailing) zone, the mass flux points inward (outward). In the highly curved region, the mass flux pushes the droplet shape to reduce its curvature. In the contact region between the droplets, the mass flux points toward the interiors of the droplets and the mass flux changes its sign across the contact interface between two droplets. In figure 3(b) the regions where the tangential component is large are shown in yellow. Figures 2 and 3 were obtained from the second example (case 2) detailed in § 9. This example involves the merging of two rising liquid droplets embedded in a continuous liquid. Finally, as
$\boldsymbol{j}_{p}$
is a divergence-free field, it does not play any role in the mass transfer equation. However, that is not the case for the energy budget equation, as we will show in the following section.
(a) Red (blue) depicts the region where the mass flux points inward (outward) and (b) yellow depicts the region where the tangential component of the mass flux is largest.
6 Energy budget
In this section, we describe the temporal evolution of the different energy components of the flow and analyse their physical meaning. The energy budget plays a central role in the scientific understanding of physical phenomena. We analyse the exchanges of the kinetic, potential, bulk and interfacial free energies. In particular, the energy budget describes the energy exchange mechanisms that govern the flow. Moreover, according to Liu & Walkington (Reference Liu and Walkington2000) and Guo, Lin & Lowengrub (Reference Guo, Lin and Lowengrub2014), the convergence of finite-dimensional approximations to the solutions of the partial differential equations is also linked to those energy transfers, particularly for non-smooth solutions.
For the sake of clarity, subscripts
$k$
,
${\it\phi}$
,
$s$
and
$p$
will refer to quantities derived from the kinetic energy, bulk free energy density
${\it\psi}_{{\it\phi}}$
, interfacial free energy density
${\it\psi}_{s}$
and potential energy, respectively. Assuming that the control volume
${\it\Omega}$
does not deform, in the following subsections we derive energy exchanges in an Eulerian description. The identities and definitions used here are detailed in appendix B.
6.1 Energies
The different energies relevant to the incompressible NSCH equations under the Boussinesq assumption in the dimensionless form are
$$\begin{eqnarray}\displaystyle & \displaystyle E_{k}=\int _{{\it\Omega}}e_{k}\,\text{d}{\it\Omega}=\int _{{\it\Omega}}{\displaystyle \frac{1}{2}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{v}\,\text{d}{\it\Omega}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle E_{{\it\phi}}=\int _{{\it\Omega}}{\it\psi}_{{\it\phi}}\,\text{d}{\it\Omega}=\int _{{\it\Omega}}\left({\displaystyle \frac{1}{2{\it\vartheta}}}({\it\phi}\ln {\it\phi}+(1-{\it\phi})\ln (1-{\it\phi}))+(1-{\it\phi}){\it\phi}\right)\,\text{d}{\it\Omega}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle E_{s}=\int _{{\it\Omega}}{\it\psi}_{s}\,\text{d}{\it\Omega}=\int _{{\it\Omega}}\left({\displaystyle \frac{Cn}{2Lm}}\,\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}\right)\,\text{d}{\it\Omega}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle E_{p}=\int _{{\it\Omega}}e_{p}\,\text{d}{\it\Omega}=\int _{{\it\Omega}}{\displaystyle \frac{Bo}{We}}{\it\phi}x_{2}\,\text{d}{\it\Omega}. & \displaystyle\end{eqnarray}$$
Consider the material derivative
$\dot{(\cdot )}$
of each energy involved in the problem, i.e.
${\dot{e}}_{\jmath }$
where
$\jmath =\{k,{\it\phi},s,p\}$
. Integrating these material derivatives over the domain results in the total derivative
$\text{d}(\cdot )/\text{d}t$
for each energy in (6.1). In the following, we consider the equations formulated for a given, undeformed control volume
${\it\Omega}$
and boundary
${\it\Gamma}$
,
$$\begin{eqnarray}{\displaystyle \frac{\text{d}E_{\jmath }}{\text{d}t}}=S_{\jmath }+Q_{\jmath }=\int _{{\it\Gamma}}{\it\epsilon}_{\jmath }^{io}\,\text{d}{\it\Gamma}+\int _{{\it\Omega}}{\it\epsilon}_{\jmath }\,\text{d}{\it\Omega},\end{eqnarray}$$
where, as before,
$\jmath =\{k,{\it\phi},s,p\}$
denotes the energy type. The term
$Q_{\jmath }$
incorporates the sources and sinks of the energy
$E_{\jmath }$
in the control volume, and
$S_{\jmath }$
is the total diffusive flux of energy across the boundary of the control volume. Integrating over time the total derivative of each energy, i.e.
$S_{\jmath }+Q_{\jmath }$
, we obtain the quantities that describe energy exchanges. We denote the extensive (intensive) quantities by uppercase (lowercase) letters, such as
$E_{(\cdot )}$
(
$e_{(\cdot )}$
) for the extensive (intensive) energies and
$\mathscr{E}_{(\cdot )}$
(
${\it\epsilon}_{(\cdot )}$
) for the extensive (intensive) energy exchanges, i.e. extensive (intensive) sources/sinks of energy.
6.2 Kinetic energy exchange
The form of the flux and source/sink terms in the case of the kinetic energy stemming from (6.2) can be derived using the momentum equations and are given by
$$\begin{eqnarray}\displaystyle Q_{k} & = & \displaystyle -\int _{{\it\Omega}}\unicode[STIX]{x1D63F}\boldsymbol{ : }[\unicode[STIX]{x1D64F}^{visc}+\unicode[STIX]{x1D64F}^{s}]\,\text{d}{\it\Omega}-\int _{{\it\Omega}}{\displaystyle \frac{Bo}{We}}{\it\phi}v_{2}\,\text{d}{\it\Omega}\nonumber\\ \displaystyle & = & \displaystyle -\int _{{\it\Omega}}\unicode[STIX]{x1D63F}\boldsymbol{ : }\left[{\displaystyle \frac{2c}{Re}}\unicode[STIX]{x1D63F}-{\displaystyle \frac{Cn}{We}}\,\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}\right]\,\text{d}{\it\Omega}-\int _{{\it\Omega}}{\displaystyle \frac{Bo}{We}}{\it\phi}v_{2}\,\text{d}{\it\Omega},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle S_{k} & = & \displaystyle \int _{{\it\Gamma}}(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D64F})\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}\nonumber\\ \displaystyle & = & \displaystyle \int _{{\it\Gamma}}\left[\boldsymbol{v}\boldsymbol{\cdot }\left({\displaystyle \frac{2c}{Re}}\unicode[STIX]{x1D63F}-{\displaystyle \frac{p}{Re}}\mathbf{1}-{\displaystyle \frac{Cn}{We}}\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}\right)\right]\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}.\end{eqnarray}$$
Physical interpretation. Regarding the volumetric term
$Q_{k}$
, the first term,
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}$
, is the rate of work done by the viscous forces of the flow. Given the constitutive model assumed for the viscous stress, (4.5), the system can be shown to lose kinetic energy from the friction between fluid particles, as the inner product of the strain rate and the viscous stress is always positive. The first term
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}$
is a sink, as can be inferred from its negative sign.
The second term
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{s}$
is the rate of work done by the capillary forces in the diffuse interface and can be thought of as follows. The vector
$\boldsymbol{{\rm\nabla}}{\it\phi}$
is parallel to the normal
$\boldsymbol{n}$
of the isosurfaces of
${\it\phi}$
. We define a local coordinate system at a point on one of these surfaces. In this local coordinate system, the capillary tensor is simply
$(\partial {\it\phi}/\partial n)^{2}$
, whereas the strain rate contribution is
$\partial v_{n}/\partial n$
, with
$v_{n}$
denoting the normal velocity. The contribution to the change of kinetic energy density around this point is then
$(\partial v_{n}/\partial n)(\partial {\it\phi}/\partial n)^{2}$
. Depending on the growth direction
$v_{n}$
, the energy exchange can behave as both a source and a sink of energy. Due to incompressibility, the pressure does not play any role in
$Q_{k}$
. Since we use divergence-conforming discretizations, the orthogonality of the pressure and the strain rate is preserved at the discrete level.
The first and second terms in the kinetic energy exchange can also be rewritten as the power expenditures carried out in an infinitesimal increment of the strain rate
$\text{d}\unicode[STIX]{x1D63F}$
done by the stresses
$(\unicode[STIX]{x1D64F}^{visc}+\unicode[STIX]{x1D64F}^{s})\boldsymbol{ : }\text{d}\unicode[STIX]{x1D63F}$
, i.e.
$Q_{k}$
is the rate of work done by the viscous and capillary stresses to achieve an increment
$\text{d}\unicode[STIX]{x1D63F}$
in strain rate. We give a detailed description of the third term when we interpret the potential energy exchanges below. Moreover, the boundary term
$S_{k}$
represents the rate of work on the boundary given by the external power done by the total stress.
6.3 Bulk free energy exchange
The form of the terms in (6.2) for the bulk free energy is
$$\begin{eqnarray}\displaystyle Q_{{\it\phi}} & = & \displaystyle \int _{{\it\Omega}}\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}\boldsymbol{\cdot }\boldsymbol{j}\,\text{d}{\it\Omega}\nonumber\\ \displaystyle & = & \displaystyle -\int _{{\it\Omega}}{\displaystyle \frac{1}{PeLn}}\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}\boldsymbol{\cdot }[M({\it\phi})\boldsymbol{{\rm\nabla}}({\it\eta}_{{\it\phi}}+{\it\eta}_{s})+\boldsymbol{{\rm\nabla}}{\it\eta}_{p}]\,\text{d}{\it\Omega},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle S_{{\it\phi}} & = & \displaystyle -\int _{{\it\Gamma}}{\it\eta}_{{\it\phi}}\,\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}\nonumber\\ \displaystyle & = & \displaystyle \int _{{\it\Gamma}}{\displaystyle \frac{1}{PeLn}}{\it\eta}_{{\it\phi}}[M({\it\phi})\boldsymbol{{\rm\nabla}}({\it\eta}_{{\it\phi}}+{\it\eta}_{s})+\boldsymbol{{\rm\nabla}}{\it\eta}_{p}]\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}.\end{eqnarray}$$
Physical interpretation. In the volumetric term
$Q_{{\it\phi}}$
, the power expenditure carried out in an infinitesimal increment of the bulk free energy gradient
$\text{d}\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}$
done by the mass flux is given by
$\boldsymbol{j}\boldsymbol{\cdot }\text{d}\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}$
. Given an infinitesimal increment
$\text{d}\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}$
, the local exchange of bulk free energy is proportional to the mass flux in the direction of that increment. That is,
$Q_{{\it\phi}}$
is the rate of work done by the mass flux to achieve an increment in bulk free energy gradient. This term may act as either a source or a sink of energy. The boundary term
$S_{{\it\phi}}$
represents the bulk free energy diffusion on the boundaries by the mass flux.
Remark 4. Substituting the mass flux equation (3.3) into (6.4) and using the Gauss divergence theorem, we find where the sources and sinks of bulk free energy come from. The volumetric
$Q_{{\it\phi}}$
and boundary
$S_{{\it\phi}}$
terms take the following form in terms of
${\it\phi}$
:
$$\begin{eqnarray}\displaystyle Q_{{\it\phi}} & = & \displaystyle -\int _{{\it\Omega}}{\displaystyle \frac{1}{PeLn}}\left\{{\it\phi}(1-{\it\phi})\left(-2+{\displaystyle \frac{1}{2{\it\vartheta}{\it\phi}(1-{\it\phi})}}\right)^{2}\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}\right.\nonumber\\ \displaystyle & & \displaystyle +\,{\displaystyle \frac{Cn}{Lm}}\left[\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}-2(1-2{\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}\right]\nonumber\\ \displaystyle & & \displaystyle +\,\left.{\displaystyle \frac{Bo}{We}}\left(-2+{\displaystyle \frac{1}{2{\it\vartheta}{\it\phi}(1-{\it\phi})}}\right)\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{e}_{2}\right\}\,\text{d}{\it\Omega},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle S_{{\it\phi}} & = & \displaystyle \int _{{\it\Gamma}}{\displaystyle \frac{1}{PeLn}}\left\{\left({\displaystyle \frac{1}{2{\it\vartheta}}}\ln {\displaystyle \frac{{\it\phi}}{1-{\it\phi}}}+1-2{\it\phi}\right)\left[\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\boldsymbol{{\rm\nabla}}{\it\phi}\right.\right.\nonumber\\ \displaystyle & & \displaystyle -\!\!\left.\left.{\displaystyle \frac{Cn}{Lm}}{\it\phi}(1-{\it\phi})\,\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}+{\displaystyle \frac{Bo}{We}}\boldsymbol{e}_{2}\right]+{\displaystyle \frac{Cn}{Lm}}\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\right\}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
The first term proportional to
$\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}$
of
$Q_{{\it\phi}}$
in (6.5) is always negative due to its negative sign, i.e. this term is a sink of energy. Although
$\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}$
is always positive, the function that multiplies it is indefinite if
${\it\vartheta}>1$
. Otherwise, if
${\it\vartheta}<1$
this term is a sink of energy. The last term involving
$f({\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}$
is indefinite since both
$f({\it\phi})$
and the Hessian of
${\it\phi}$
are indefinite. Normalizing
$\boldsymbol{{\rm\nabla}}{\it\phi}$
, this quadratic form
$\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}$
geometrically describes how the curvature in the
$\boldsymbol{{\rm\nabla}}{\it\phi}$
direction changes as we move along this direction, i.e. how the curvature changes across the interface.
6.4 Interfacial free energy exchange
The exchanges (6.2) in interfacial free energy can be expressed as
$$\begin{eqnarray}\displaystyle Q_{s} & = & \displaystyle \int _{{\it\Omega}}{\displaystyle \frac{Cn}{Lm}}[\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{j}-\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}]\,\text{d}{\it\Omega}\nonumber\\ \displaystyle & = & \displaystyle -\!\int _{{\it\Omega}}{\displaystyle \frac{Cn}{Lm}}\left[{\displaystyle \frac{1}{PeLn}}\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}(M({\it\phi})\boldsymbol{{\rm\nabla}}({\it\eta}_{{\it\phi}}+{\it\eta}_{s})+\boldsymbol{{\rm\nabla}}{\it\eta}_{p})+\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}\right]\text{d}{\it\Omega},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle S_{s} & = & \displaystyle -\int _{{\it\Gamma}}{\displaystyle \frac{Cn}{Lm}}(\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{j})\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}\nonumber\\ \displaystyle & = & \displaystyle \int _{{\it\Gamma}}{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Cn}{Lm}}[\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(M({\it\phi})\,\boldsymbol{{\rm\nabla}}({\it\eta}_{{\it\phi}}+{\it\eta}_{s})+\boldsymbol{{\rm\nabla}}{\it\eta}_{p})]\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}.\end{eqnarray}$$
Physical interpretation. The first term of the volumetric term
$Q_{s}$
has a similar meaning to that given for the kinetic and bulk free energy. The power expenditure carried out in an infinitesimal increment in the curvature of the phase field,
$\text{d}\unicode[STIX]{x1D643}$
, done by the mass flux gradient is given by
$\boldsymbol{{\rm\nabla}}\boldsymbol{j}\boldsymbol{ : }\text{d}\unicode[STIX]{x1D643}$
. Given an infinitesimal increment
$\text{d}\unicode[STIX]{x1D643}$
, the local exchange of interfacial free energy is proportional to the mass flux gradient. That is,
$Q_{s}$
is the rate of work done by the mass flux gradient to achieve an increment in the curvature of the phase field. This interfacial free energy exchange acts as a source or a sink of energy.
We can interpret the first term in
$Q_{s}$
in a simple way. Suppose that the isosurface of
${\it\phi}$
is locally flat, and that the mass flux varies only along the normal direction to the surface. Thus, the mass flux gradient leads to different numbers of particles diffusing at different points along the normal, therefore displacing the isosurfaces of
${\it\phi}$
. This in turn leads to a change in the magnitude of the gradient of
${\it\phi}$
, represented by the corresponding components of
$\unicode[STIX]{x1D643}$
, thus changing the interfacial free energy. Now suppose that the mass flux varies along the surface of constant
${\it\phi}$
. Thus, the mass flux gradient leads to a different number of diffusing particles along the tangential direction of the isosurface. This in turn induces a change in the curvature of the isosurface, thus changing the magnitude of the gradient of
${\it\phi}$
, which leads to a change in interfacial free energy. This is represented by the components of the Hessian of
${\it\phi}$
describing the curvature of phase isosurfaces.
The second term in
$Q_{s}$
appears in the kinetic energy exchange. Even though it is discussed there, we interpret it here in the context of the interfacial energy. Looking at a small region on a phase field surface, if the normal velocity increases (decreases) in the direction of the isosurface normal (which is parallel to
$\boldsymbol{{\rm\nabla}}{\it\phi}$
), then the interface widens (narrows), thereby decreasing (increasing)
$|\boldsymbol{{\rm\nabla}}{\it\phi}|$
and thus becoming a sink (source) of interfacial energy. This term differs from that found in the kinetic energy exchange in scale and sign. As a consequence of the change of sign, if this term is a source (sink) of interfacial free energy, it will be a sink (source) of kinetic energy. The boundary term
$S_{s}$
represents the interfacial free energy diffusion on the boundaries by the mass flux.
Remark 5. By substituting the mass flux from (3.3) into (6.6) and using the Gauss divergence theorem we find where the source and sink of interfacial free energy come from. The volumetric
$Q_{s}$
and boundary
$S_{s}$
terms take the following forms in terms of
${\it\phi}$
:
$$\begin{eqnarray}\displaystyle Q_{s} & = & \displaystyle -\int _{{\it\Omega}}\left\{{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Cn}{Lm}}\left[-2(1-2{\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}+\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}\right.\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.\left.{\displaystyle \frac{Cn}{Lm}}{\it\phi}(1-{\it\phi})(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\right]+{\displaystyle \frac{Cn}{Lm}}\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}\right\}\,\text{d}{\it\Omega},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle S_{s} & = & \displaystyle \int _{{\it\Gamma}}{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Cn}{Lm}}\left\{\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\left[\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \left(-2(1-2{\it\phi})\boldsymbol{{\rm\nabla}}{\it\phi}-{\displaystyle \frac{Cn}{Lm}}(1-2{\it\phi})\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\right)\right.\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\unicode[STIX]{x1D643}-{\displaystyle \frac{Cn}{Lm}}{\it\phi}(1-{\it\phi})\boldsymbol{{\rm\nabla}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\right]\nonumber\\ \displaystyle & & \displaystyle +\,\left.{\displaystyle \frac{Cn}{Lm}}{\it\phi}(1-{\it\phi})(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\boldsymbol{\cdot }\unicode[STIX]{x1D643}\right\}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}.\end{eqnarray}$$
There is only one term that does not change its sign, that is, the contraction on itself of
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}$
. This term always acts as a sink. All the remaining terms that appear in
$Q_{s}$
in (6.7) are indefinite, i.e. they may be either a source or a sink of interfacial free energy. The term
$\propto f_{1}({\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}+f_{2}({\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}$
also appears in the bulk free energy exchange
$Q_{{\it\phi}}$
, with the same sign and scales. Thus, its physical meaning is exactly the same.
6.5 Potential energy exchange
Finally, the forms of the terms in (6.2) for the potential energy are
$$\begin{eqnarray}\displaystyle Q_{p} & = & \displaystyle \int _{{\it\Omega}}{\displaystyle \frac{Bo}{We}}\boldsymbol{e}_{2}\boldsymbol{\cdot }\boldsymbol{j}\,\text{d}{\it\Omega}+\int _{{\it\Omega}}{\displaystyle \frac{Bo}{We}}{\it\phi}v_{2}\,\text{d}{\it\Omega}\nonumber\\ \displaystyle & = & \displaystyle -\int _{{\it\Omega}}{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Bo}{We}}\boldsymbol{e}_{2}\boldsymbol{\cdot }[M({\it\phi})\boldsymbol{{\rm\nabla}}({\it\eta}_{{\it\phi}}+{\it\eta}_{s})+\boldsymbol{{\rm\nabla}}{\it\eta}_{p}]\,\text{d}{\it\Omega}+\int _{{\it\Omega}}{\displaystyle \frac{Bo}{We}}{\it\phi}v_{2}\,\text{d}{\it\Omega},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle S_{p} & = & \displaystyle -\int _{{\it\Gamma}}{\displaystyle \frac{Bo}{We}}x_{2}\;\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}\nonumber\\ \displaystyle & = & \displaystyle \int _{{\it\Gamma}}{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Bo}{We}}x_{2}[M({\it\phi})\boldsymbol{{\rm\nabla}}({\it\eta}_{{\it\phi}}+{\it\eta}_{s})+\boldsymbol{{\rm\nabla}}{\it\eta}_{p}]\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}.\end{eqnarray}$$
Physical interpretation. The sources and sinks of potential energy described by
$Q_{p}$
are related to relative motion in the
$\boldsymbol{e}_{2}$
direction. The first term describes the diffusive flux in this direction and therefore can be a source or a sink of potential energy. This contribution acts as a source – if the flux points in the same direction as
$\boldsymbol{e}_{2}$
– since it leads to fluid particles diffusing against the pull of gravity (pointing towards
$-\boldsymbol{e}_{2}$
), thereby gaining potential energy. This contribution acts as a sink of energy if the direction of the mass flux is parallel to gravity, thus accepting the pull and losing potential energy. The second term is related to the non-diffusive motion of the fluid particles. If the velocity points towards
$\boldsymbol{e}_{2}$
then the particle moves against gravity, gaining potential energy and vice versa. This term cancels out the third term appearing in the kinetic energy exchanges and therefore it does not play any role in the energy budget equation.
Remark 6. As we did in Remarks 4 and 5, we obtain explicit expressions in terms of
${\it\phi}$
for the volumetric
$Q_{p}$
and boundary
$S_{p}$
terms:
$$\begin{eqnarray}\displaystyle Q_{p} & = & \displaystyle -\int _{{\it\Omega}}{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Bo}{We}}\left[\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{e}_{2}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.{\displaystyle \frac{Cn}{Lm}}(1-2{\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{e}_{2}+{\displaystyle \frac{Bo}{We}}\right]\text{d}{\it\Omega}+\int _{{\it\Omega}}{\displaystyle \frac{Bo}{We}}{\it\phi}v_{2}\,\text{d}{\it\Omega},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle S_{p} & = & \displaystyle \int _{{\it\Gamma}}{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Bo}{We}}\left\{{\displaystyle \frac{Cn}{Lm}}{\it\phi}(1-{\it\phi})\boldsymbol{e}_{2}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.x_{2}\left[\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\boldsymbol{{\rm\nabla}}{\it\phi}-{\displaystyle \frac{Cn}{Lm}}{\it\phi}(1-{\it\phi})\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}+{\displaystyle \frac{Bo}{We}}\boldsymbol{e}_{2}\right]\right\}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Despite identifying at least one negative-definite term (sink of energy) in the energy exchanges discussed above, the definiteness of the potential energy exchanges cannot be determined beforehand. Finally, normalizing
$\boldsymbol{{\rm\nabla}}{\it\phi}$
, the term
$\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{e}_{2}$
geometrically describes how the curvature in the
$\boldsymbol{{\rm\nabla}}{\it\phi}$
direction changes as we move along the
$\boldsymbol{e}_{2}$
direction.
6.6 Total energy budget
After detailing the different energy exchange terms in the system, we address the total energy exchanges in the system. Integrating the source and sink terms along with the fluxes over time leads to the total energy exchange terms, where the volumetric terms are defined as
$$\begin{eqnarray}\mathscr{E}_{\jmath }=\int _{t}Q_{\jmath }\,\text{d}t=\int _{t}\int _{{\it\Omega}}{\it\epsilon}_{\jmath }\,\text{d}{\it\Omega}\,\text{d}t,\quad \jmath =\{k,{\it\phi},s,p\},\end{eqnarray}$$
and the input/output of energy diffused across the boundaries is
$$\begin{eqnarray}\mathscr{E}_{\jmath }^{io}|_{\mathit{diff}}=\int _{t}S_{\jmath }\,\text{d}t=\int _{t}\int _{{\it\Gamma}}{\it\epsilon}_{\jmath }^{io}|_{\mathit{diff}}\,\text{d}{\it\Gamma}\,\text{d}t,\quad \jmath =\{k,{\it\phi},s,p\}.\end{eqnarray}$$
If there exists a momentum flux across the boundaries, the energy must be advected according to
$$\begin{eqnarray}\mathscr{E}_{\jmath }^{io}|_{\mathit{adv}}=-\int _{t}\int _{{\it\Gamma}}e_{\jmath }\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}\,\text{d}t,\quad \jmath =\{k,{\it\phi},s,p\},\end{eqnarray}$$
and the sum
$\mathscr{E}_{\jmath }^{\text{io}}|_{\mathit{diff}}+\mathscr{E}_{\jmath }^{\text{io}}|_{\mathit{adv}}$
yields the total input/output of energy across the boundaries.
Finally, the energy budget is defined as
$$\begin{eqnarray}E_{k}+E_{{\it\phi}}+E_{s}+E_{p}-(\mathscr{E}_{k}+\mathscr{E}_{{\it\phi}}+\mathscr{E}_{s}+\mathscr{E}_{p}+\mathscr{E}_{k}^{\text{io}}+\mathscr{E}_{{\it\phi}}^{\text{io}}+\mathscr{E}_{s}^{\text{io}}+\mathscr{E}_{p}^{\text{io}})=C,\end{eqnarray}$$
for an arbitrary constant
$C$
. We analyse the energy exchange over the whole domain and use no-flux boundary conditions for all fields in the numerical examples discussed in the following. Under such a set-up, the energy budget equation reads
$$\begin{eqnarray}\mathop{\sum }_{\jmath =\{k,{\it\phi},s,p\}}(E_{\jmath }-\mathscr{E}_{\jmath })(t)=E_{k}+E_{{\it\phi}}+E_{s}+E_{p}-(\mathscr{E}_{k}+\mathscr{E}_{{\it\phi}}+\mathscr{E}_{s}+\mathscr{E}_{p})=C.\end{eqnarray}$$
As mentioned before, we relate
$E_{\jmath }$
and
$\mathscr{E}_{\jmath }$
,
$\jmath =\{k,{\it\phi},s,p\}$
, to extensive energies and their energy exchanges, whereas
$e_{\jmath }$
and
${\it\epsilon}_{\jmath }$
are associated with intensive energies and their energy exchanges, respectively.
Remarks 4–6 show which terms work as sinks or as both sources and sinks of energy, as functions of the phase field. Concluding, the sum of energy exchanges in terms of the phase field reads
$$\begin{eqnarray}\displaystyle \mathop{\sum }_{\jmath =\{k,{\it\phi},s,p\}}{\it\epsilon}_{\jmath } & = & \displaystyle {\displaystyle \frac{Cn(Lm-We)}{WeLm}}\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}-{\displaystyle \frac{2c}{Re}}\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D63F}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{{\it\phi}(1-{\it\phi})}{PeLn}}\left(-2+{\displaystyle \frac{1}{2{\it\vartheta}{\it\phi}(1-{\it\phi})}}\right)^{2}\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{{\it\phi}(1-{\it\phi})}{PeLn}}\left({\displaystyle \frac{Cn}{Lm}}\right)^{2}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{2}{PeLn}}{\displaystyle \frac{Cn}{Lm}}\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}\nonumber\\ \displaystyle & & \displaystyle +\,{\displaystyle \frac{4}{PeLn}}{\displaystyle \frac{Cn}{Lm}}(1-2{\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{1+{\it\phi}(1-{\it\phi})}{PeLn}}{\displaystyle \frac{Bo}{We}}\left(-2+{\displaystyle \frac{1}{2{\it\vartheta}{\it\phi}(1-{\it\phi})}}\right)\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{e}_{2}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Bo}{We}}\left[(1-2{\it\phi}){\displaystyle \frac{Cn}{Lm}}\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{e}_{2}+{\displaystyle \frac{Bo}{We}}\right].\end{eqnarray}$$
In the system governed by NSCH equations, there are three volumetric terms that always work as energy sinks. They are the terms proportional to
$f_{1}({\it\phi})\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D63F}$
,
$f_{2}({\it\phi})\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}$
and
$f_{3}({\it\phi})(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})$
. If the critical temperature ratio
${\it\vartheta}$
is less than one, the term proportional to
$f_{4}({\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}$
is also a sink. The functions
$f_{1}$
,
$f_{2}$
,
$f_{3}$
and
$f_{4}$
are stated as
$$\begin{eqnarray}\displaystyle & \displaystyle f_{1}({\it\phi})=-{\displaystyle \frac{2c}{Re}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f_{2}({\it\phi})=-{\displaystyle \frac{{\it\phi}(1-{\it\phi})}{PeLn}}\left(-2+{\displaystyle \frac{1}{2{\it\vartheta}{\it\phi}(1-{\it\phi})}}\right)^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f_{3}({\it\phi})=-{\displaystyle \frac{1}{PeLn}}\left({\displaystyle \frac{Cn}{Lm}}\right)^{2}{\it\phi}(1-{\it\phi}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f_{4}({\it\phi})=-2{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Cn}{Lm}}\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right). & \displaystyle\end{eqnarray}$$
Figure 4 depicts the behaviour of the term
$f({\it\phi}(x),{\it\vartheta})\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}=[-2{\it\phi}(1-{\it\phi})+(1/2{\it\vartheta})]\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}$
that appears in both the bulk and the interfacial free energies, along the
${\it\phi}$
axis for
${\it\vartheta}$
from 0.5 to
$\infty$
. The phase field is defined by a hyperbolic tangent,
${\it\phi}(x)=0.95[0.5\tanh (10x)+0.5]+0.025$
. The grey region represents the interface length, whereas the red region depicts the region where this function lives when phase segregation takes place, i.e.
${\it\vartheta}>1$
. The term
$f({\it\phi}(x),{\it\vartheta})\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}$
is always positive if
${\it\vartheta}<1$
. After scaling by
$-2(1/PeLn)(Cn/Lm)$
, the energy exchange term becomes negative definite, thus it is a sink of energy. However, if
${\it\vartheta}>1$
(red region depicted in figure 4), the energy exchange term is indefinite, being a source of energy in the middle region of the interface and becoming a sink of energy as we move away from the interface. When the phase is segregated, far away from the interfaces, this term does not play any role. In the extreme case,
${\it\vartheta}\rightarrow \infty$
, this term acts as a source of energy.
The plot depicts the function that appears in both the bulk and the interfacial free energies,
$f({\it\phi},{\it\vartheta})\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}=[-2{\it\phi}(1-{\it\phi})+1/2{\it\vartheta}]\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}$
, along the
${\it\phi}$
axis for
${\it\vartheta}$
from 0.5 to
$\infty$
. The phase field is defined by a hyperbolic tangent,
${\it\phi}(x)=0.95[0.5\,\text{tanh}(10x)+0.5]+0.025$
. The grey region represents the interface length
$l$
, whereas the red region depicts the region where this function lives when the phase segregation is assumed, i.e.
${\it\vartheta}>1$
.
Finally, the sum of boundary sources and sinks of energies, as functions of the phase field, reads
$$\begin{eqnarray}\displaystyle \mathop{\sum }_{\jmath =\{k,{\it\phi},s,p\}}{\it\epsilon}_{\jmath }^{io}|_{\mathit{diff}}\! & = & \displaystyle \!\left\{\boldsymbol{v}\boldsymbol{\cdot }\left({\displaystyle \frac{2c}{Re}}\unicode[STIX]{x1D63F}-{\displaystyle \frac{p}{Re}}\mathbf{1}-{\displaystyle \frac{Cn}{We}}\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\,{\displaystyle \frac{1}{PeLn}}\left({\displaystyle \frac{1}{2{\it\vartheta}}}\ln {\displaystyle \frac{{\it\phi}}{1-{\it\phi}}}+1-2{\it\phi}\right)\left[\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\boldsymbol{{\rm\nabla}}{\it\phi}\right.\nonumber\\ \displaystyle & & \displaystyle -\!\left.{\displaystyle \frac{Cn}{Lm}}{\it\phi}(1-{\it\phi})\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}+{\displaystyle \frac{Bo}{We}}\boldsymbol{e}_{2}\right]\!+{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Cn}{Lm}}\!\left(\!-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\!\right)\!\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{1}{PeLn}}{\displaystyle \frac{Cn}{Lm}}\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\left[\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \left(2(1-2{\it\phi})\boldsymbol{{\rm\nabla}}{\it\phi}+{\displaystyle \frac{Cn}{Lm}}(1-2{\it\phi})\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\right)\right.\nonumber\\ \displaystyle & & \displaystyle -\,\left.\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\unicode[STIX]{x1D643}+{\displaystyle \frac{Cn}{Lm}}{\it\phi}(1-{\it\phi})\boldsymbol{{\rm\nabla}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\right]\nonumber\\ \displaystyle & & \displaystyle +\,{\displaystyle \frac{{\it\phi}(1-{\it\phi})}{PeLn}}\left({\displaystyle \frac{Cn}{Lm}}\right)^{2}(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\boldsymbol{\cdot }\unicode[STIX]{x1D643}+{\displaystyle \frac{{\it\phi}(1-{\it\phi})}{PeLn}}{\displaystyle \frac{Bo}{We}}{\displaystyle \frac{Cn}{Lm}}\boldsymbol{e}_{2}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\nonumber\\ \displaystyle & & \displaystyle +\,{\displaystyle \frac{x_{2}}{PeLn}}{\displaystyle \frac{Bo}{We}}\left[\left(-2{\it\phi}(1-{\it\phi})+{\displaystyle \frac{1}{2{\it\vartheta}}}\right)\boldsymbol{{\rm\nabla}}{\it\phi}\right.\nonumber\\ \displaystyle & & \displaystyle -\,{\it\phi}(1-{\it\phi})\left.\left.{\displaystyle \frac{Cn}{Lm}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}+{\displaystyle \frac{Bo}{We}}\boldsymbol{e}_{2}\right]\right\}\boldsymbol{\cdot }\,\boldsymbol{n}.\end{eqnarray}$$
In our numerical experiments, we present the energy terms using (6.1a–d ) and their related energy exchange terms using (6.3), (6.4), (6.6) and (6.8). We have opted to compute energy exchanges using the expressions containing the mass flux as a variable, in those terms related to bulk, interfacial and potential energies, i.e. (6.4), (6.6) and (6.8) instead of (6.5), (6.7) and (6.9) presented in Remarks 4–6. The reasons to opt for this approach are twofold. First, we take advantage of the mixed formulation of Cahn–Hilliard equation, since the auxiliary variable is the chemical potential. Second, at most second-order derivatives are required for the phase field, whereas (6.5), (6.7) and (6.9) require fourth-order derivatives of the phase field, which leads to lower-order approximations of the fields of interest.
To simplify the exposition of results, we split the volumetric energy exchange terms. Three terms appear in the kinetic energy exchanges: the power expenditure done by the viscous stress
${\it\epsilon}_{k}^{visc}=-\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}$
, the power expenditure done by the capillary stress
${\it\epsilon}_{k}^{s}=-\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{s}$
and a term related to the Boussinesq approximation
${\it\epsilon}_{k}^{buoy}=-(Bo/We){\it\phi}v_{2}$
. The bulk free energy exchange is described by only one term, that is, the power expenditure done by the mass flux on the bulk free energy gradient,
${\it\epsilon}_{{\it\phi}}=\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}\boldsymbol{\cdot }\boldsymbol{j}$
. The interfacial free energy exchange is split into a term related to changes in the curvature done by the mass flux gradient,
${\it\epsilon}_{s}^{curv}=(Cn/Lm)\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{j}$
, and a term that relates the strain rate and capillary effects,
${\it\epsilon}_{s}^{s}=-(Cn/Lm)\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}$
, similarly to
${\it\epsilon}_{k}^{s}$
. Finally, in the potential energy exchange, the first term describes the mass flux in the vertical direction
${\it\epsilon}_{p}^{mass}=(Bo/We)\boldsymbol{e}_{2}\boldsymbol{\cdot }\boldsymbol{j}$
, while the second term
${\it\epsilon}_{p}^{buoy}=-{\it\epsilon}_{k}^{buoy}=(Bo/We){\it\phi}\,v_{2}$
, which also appears in the kinetic energy exchange, is related to the Boussinesq approximation. The same notation is employed for extensive quantities. Finally, we employ the trapezoidal rule to integrate over time the energy exchanges.
7 Numerical scheme: divergence-conforming B-spline spaces
Isogeometric analysis has been used successfully to solve high-order phase-field models, including the Cahn–Hilliard equation (Gómez et al. Reference Gómez, Calo, Bazilevs and Hughes2008; Vignal, Collier & Calo Reference Vignal, Collier and Calo2013), the Navier–Stokes–Korteweg equations (Gomez et al. Reference Gomez, Hughes, Nogueira and Calo2010), the Swift–Hohenberg equation (Gomez & Nogueira Reference Gomez and Nogueira2012a ) and the phase-field crystal equation (Gomez & Nogueira Reference Gomez and Nogueira2012b ; Vignal et al. Reference Vignal, Collier and Calo2013, Reference Vignal, Dalcin, Brown, Collier and Calo2015a ). Stability conditions in models for electromagnetism such as the Maxwell equations (Buffa, Sangalli & Vázquez Reference Buffa, Sangalli and Vázquez2010) and flow models including the Stokes (Buffa, De Falco & Sangalli Reference Buffa, De Falco and Sangalli2011a ; Evans & Hughes Reference Evans and Hughes2013a ) and Navier–Stokes equations (Evans & Hughes Reference Evans and Hughes2013b ) have been solved exactly using curl- and divergence-conforming spaces, respectively. Solving incompressible flow models using divergence-conforming spaces produces discrete pointwise divergence-free velocity fields. The advantages of such fields in the conservation of kinetic energy, vorticity, enstrophy and helicity are discussed by Evans & Hughes (Reference Evans and Hughes2013c ) for the Navier–Stokes equations and are extensively exploited herein.
The discrete model is solved using the PetIGA-MF (Sarmiento et al. Reference Sarmiento, Cortes, Garcia, Dalcin, Collier and Calo2015; Vignal et al. Reference Vignal, Sarmiento, Côrtes, Dalcin and Calo2015b ), a high-performance framework built on top of PetIGA (Côrtes et al. Reference Côrtes, Vignal, Sarmiento, García, Collier, Dalcin and Calo2014; Dalcin et al. Reference Dalcin, Collier, Vignal, Côrtes and Calo2015), which uses structure-preserving B-spline basis functions. This framework simplifies the solution of systems of high-order partial differential equations, where multifield strategies can provide a high order of approximation and smoothness in the basis functions, as well as structure-preserving discretizations that allow exact fulfilment of discrete stability conditions (Buffa et al. Reference Buffa, Rivas, Sangalli and Vázquez2011b ; Côrtes et al. Reference Côrtes, Coutinho, Dalcin and Calo2015).
The idea of structure-preserving spaces is based on satisfying the exact sequence given by the discrete version of the de Rham diagram (Buffa et al. Reference Buffa, Rivas, Sangalli and Vázquez2011b ). Here, two spaces are said to conform to an operator if they satisfy a step in the sequence that corresponds to that particular operator. PetIGA-MF admits multifield discretizations and provides gradient-, curl-, divergence- and integral-conforming discrete spaces, allowing the user to discretize a specific problem in a stable manner.
We use a mixed formulation of the Cahn–Hilliard equation, taking the chemical potential
${\it\eta}$
as an auxiliary variable. We do this to reduce the computational cost of using high-order and high-continuity basis functions (Collier et al.
Reference Collier, Pardo, Dalcin, Paszynski and Calo2012, Reference Collier, Dalcin, Pardo and Calo2013; Collier, Dalcin & Calo Reference Collier, Dalcin and Calo2014), and to avoid complications with the imposition of nonlinear boundary conditions that arise from the discretization of the mass flux in the primal form. Thus, we have four variables, the velocity, pressure, phase and chemical potential fields. We discretize the velocity and pressure variables using a divergence- and integral-conforming conjugated pair of spaces that satisfy the inf-sup stability condition exactly and render a pointwise divergence-free discrete velocity field. The phase field and chemical potential variables are discretized using
$H^{1}$
spaces. No penetration and free-slip boundary conditions are imposed on the velocity since we are not seeking to study boundary-layer effects in this work. In addition, homogeneous Neumann boundary conditions are applied to all the other variables.
We employ the generalized-
${\it\alpha}$
time integrator presented by Chung & Hulbert (Reference Chung and Hulbert1993), Jansen, Whiting & Hulbert (Reference Jansen, Whiting and Hulbert2000) and Espath et al. (Reference Espath, Braun, Awruch and Dalcin2015a
) to advance all fields in time. In all examples presented here, we employ a spectral radius,
${\it\rho}_{\infty }$
, equal to 0.9. This means that we add a small amount of numerical dissipation in the system to avoid numerical instabilities. We also use time step adaptivity based on keeping the local truncation error of the method under a prescribed tolerance.
8 Problem set-up
We perform four numerical experiments to investigate the energy exchanges in the NSCH system. In addition, we analyse the main features of this kind of flow as well as the time evolution of the half-distance between the meniscus, i.e. the bridge length as defined in Eggers, Lister & Stone (Reference Eggers, Lister and Stone1999). The first two-dimensional (2-D) simulation (case 1) is performed without buoyancy effects, whereas the second and third 2-D simulations (cases 2 and 3) take into account buoyancy effects to analyse the merging of droplets when they are rising. The fourth simulation (case 4), which is three-dimensional (3-D), includes buoyancy effects and deals with a rising droplet that coalesces with a film of fluid.
Table 2 lists (in the first column, case 1) the physical parameters employed in the simulation without buoyancy and (in the second and third columns, cases 2 and 3, respectively) the physical parameters employed in each simulation that takes into account buoyancy. Cases 2 and 3 differ only in the dimensionless
$Cn/We$
ratio. We employ
$(Cn/We)^{-1}=10^{2}$
and
$10^{3}$
for those simulations, respectively. The fourth column of table 2 lists the parameters used in case 4.
The domains we consider are shown in figure 5, where the vertical axis is the
$x_{2}$
direction. Figure 5(a) shows the 2-D domain and the initial location of the droplets employed in case 1. Figure 5(b) depicts the domain and the initial location of the droplets employed in cases 2 and 3. Figure 5(c) shows the domain with the initial configuration for the 3-D case 4. Due to symmetry, only a quarter of the domain is considered in case 4.
Table 3 lists (from the first to third columns) the number of mesh nodes used in the 2-D cases 1, 2 and 3, respectively. The last column lists the number of mesh nodes employed in the 3-D case 4. Table 4 shows the polynomial degree and continuity employed. All meshes are uniform.
Physical parameters.
Meshes.
In the first simulation, case 1, the droplets and the interstitial fluid are endowed with the same density, but the droplet is less viscous than the interstitial fluid. In the remaining cases, 2–4, droplets are lighter and less viscous (endowed with density
${\it\rho}_{1}$
and viscosity
${\it\mu}_{1}$
) than the interstitial fluid (endowed with density
${\it\rho}_{2}$
and viscosity
${\it\mu}_{2}$
). To represent the phobic interactions between the phases,
${\it\vartheta}=3/2>1$
is employed.
Finally, the initial condition in case 1 is given by
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}h=0.01;\quad b_{1}=1.0;\quad b_{2}=0.75;\\ a_{1}=\tanh \left\{\!\left[\!\left({\displaystyle \frac{x_{1}-0.375}{b_{1}}}\right)^{2}\!+\left({\displaystyle \frac{x_{2}-0.375-0.15}{b_{2}}}\right)^{2}\!-0.15^{2}\right]\!{\displaystyle \frac{h^{-1}}{0.3+h}}\!\right\}0.499-0.5;\\ a_{2}=\tanh \left\{\!\left[\!\left({\displaystyle \frac{x_{1}-0.375}{b_{1}}}\right)^{2}\!+\left({\displaystyle \frac{x_{2}-0.375+0.15}{b_{2}}}\right)^{2}\!-0.15^{2}\right]\!{\displaystyle \frac{h^{-1}}{0.3+h}}\!\right\}0.499-0.5;\\ {\it\phi}(x_{1},x_{2},t=0)=a_{1}(x_{1},x_{2})+a_{2}(x_{1},x_{2})+1;\end{array}\right\}\qquad\end{eqnarray}$$
while in cases 2 and 3 they are
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}h=0.01;\\ a_{1}=\tanh \left\{[(x_{1}-0.375)^{2}+(x_{2}-0.65)^{2}-0.15^{2}]{\displaystyle \frac{h^{-1}}{0.3+h}}\right\}0.499-0.5;\\ a_{2}=\tanh \left\{[(x_{1}-0.375)^{2}+(x_{2}-0.30)^{2}-0.10^{2}]{\displaystyle \frac{h^{-1}}{0.2+h}}\right\}0.499-0.5;\\ {\it\phi}(x_{1},x_{2},t=0)=a_{1}(x_{1},x_{2})+a_{2}(x_{1},x_{2})+1.\end{array}\right\}\end{eqnarray}$$
In case 4, the initial condition is defined as
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}h=0.02;\\ a_{1}=\tanh \left\{\vphantom{{\displaystyle \frac{h^{-1}}{0.3+h}}}[(x_{1}-0.375)^{2}+(x_{2}-0.5)^{2}+(x_{3}-0.375)^{2}-0.15^{2}]\right.\\ \times \,\left.{\displaystyle \frac{h^{-1}}{0.3+h}}\right\}0.499-0.5;\hspace{60.0pt}\\ a_{2}=-\!\tanh \left({\displaystyle \frac{x_{2}-1.25}{h}}\right)0.499+0.5;\\ {\it\phi}(x_{1},x_{2},x_{3},t=0)=a_{1}(x_{1},x_{2},x_{3})+a_{2}(x_{2})+1.\end{array}\right\}\qquad\end{eqnarray}$$
Initial conditions in: (a) the 2-D domain without buoyancy effects, case 1; (b) the 2-D domain with buoyancy effects, cases 2 and 3; and (c) the 3-D domain with buoyancy effects, case 4.
Discrete spaces.
9 Numerical investigations
In this section, we present the four numerical experiments described in the previous section. These use simple domains and focus on the physical aspects of the energy exchanges. The first simulation deals with the merging of droplets without gravitational effects in two dimensions. The second and third simulations deal with rising droplets that merge as they evolve. The fourth and final simulation is 3-D and deals with a single droplet rising to merge with a thin film of fluid. In the first simulation, the droplets are less viscous than the interstitial fluid and both fluids have the same density. In the remaining simulations, the droplets are lighter and less viscous than the interstitial fluid.
9.1 Two-dimensional investigation: case 1
In the absence of gravity, the flow is driven by surface tension and liquid droplets in close proximity immersed in a liquid continuum merge into a single bigger droplet.
9.1.1 General flow features in droplet dynamics: case 1
Let us define two characteristic regions in the merging of droplets. The first one is the region where the interface of the droplets breaks down to combine them into a single entity. This region does not contain the interface, only the region between the broken interfaces. We call it the confluence region. The second one is the broken interface, which has the highest curvature. We call it the meniscus region. These regions are depicted of figure 6(b). In figure 6(a), we depict the bridge length and the meniscus curvature, following the usual definition by Eggers et al. (Reference Eggers, Lister and Stone1999).
A schematic representation of the meniscus and confluence regions defined over the isocurve
${\it\phi}=0.5$
during the merging of droplets.
Initially, both droplets are elliptical. To achieve a state that requires minimum energy, they merge into a bigger droplet, which minimizes the surface area. This results in a circular droplet. Figure 7 shows isocurves of the phase-field evolution, including the initial state, two elliptical droplets and the final state, a steady bigger circular droplet.
Isocurves for
${\it\phi}=0.1$
, 0.5 and 0.9.
9.1.2 Meniscus and bridge evolution: case 1
According to the theory presented by Eggers et al. (Reference Eggers, Lister and Stone1999), at early stages of the merging process, the bridge length (half the distance between the menisci) scales with
$t\ln t$
, i.e.
$r(t)\propto t\ln t$
, and this scaling holds up to
$r(t)\lesssim 0.03R$
, where
$t$
is the time from the initial contact between the droplets and
$R$
is the droplet radius. In Eggers et al. (Reference Eggers, Lister and Stone1999), the authors used the Stokes system together with the sharp-interface method to track the interface and model the viscous motion. If
$r\gtrsim 0.03R$
the scaling law changes and the bridge length scales with
$\sqrt{t}$
, i.e.
$r(t)\propto \sqrt{t}$
. Although the theoretical scaling law at early stages of the merging is
${\sim}t\ln t$
, a linear growth has been reported by Aarts et al. (Reference Aarts, Lekkerkerker, Guo, Wegdam and Bonn2005) and Thoroddsen, Takehara & Etoh (Reference Thoroddsen, Takehara and Etoh2005) in experimental analyses.
Figure 8 depicts the meniscus evolution defined on the isocurve
${\it\phi}=0.5$
during the merging. Figure 9 shows the bridge length evolution, measured from the symmetry axis to the meniscus point. Due to the fact that the scaling law for the early stage,
$r\lesssim 0.03R$
, holds for a very tiny period of time, we restrict ourselves to analysing the bridge length for
$r\gtrsim 0.03R$
. Thus, we fit the bridge length obtained numerically, depicted by the red line in figure 9, with a fitted function proportional to
$\sqrt{t}$
, depicted by the green line in the same figure.
Visualization of the meniscus and bridge evolution on
${\it\phi}=0.5$
, for case 1, from
$t_{i}=0.2698$
to
$t_{f}=t_{i}+0.055$
using equally spaced isocurves over time.
Bridge length evolution
$r(t)$
in case 1. Red line is the numerical result obtained herein,
$r(t)$
. Green line is the fitted square root curve,
$\bar{r}(t)\propto \sqrt{t}$
.
The simulation results, obtained from the NSCH equations, fit with the scaling law
${\sim}\sqrt{t}$
, also verified by the Stokes equation coupled to a sharp-interface method (Eggers et al.
Reference Eggers, Lister and Stone1999) and by experimental results (Aarts et al.
Reference Aarts, Lekkerkerker, Guo, Wegdam and Bonn2005; Thoroddsen et al.
Reference Thoroddsen, Takehara and Etoh2005). To further verify our model, we also fit the bridge length
$r(t)$
with a function proportional to
$t^{b}$
,
$b$
being a coefficient to be determined. The resulting coefficient is
$b=0.504$
. As mentioned before, the theoretical prediction is
$b=0.5$
to yield the function
${\sim}\sqrt{t}$
. This allows us to conclude that there is excellent agreement between our model and the theoretical/experimental models in the available literature.
9.1.3 Energy budget of the flow: case 1
Figure 10 presents the extensive energies and their respective energy exchanges versus time, for case 1. In figure 10(a) the kinetic energy and its energy exchanges are shown. During the merging, the kinetic energy (solid blue line,
$E_{k}$
) shows a peak which is quickly damped. The source of kinetic energy yielding that peak is provided by capillary effects and damped by the viscous ones. Energy exchange terms are depicted using dashed lines. The sink of the energy done by the viscous stress is depicted by the dashed red line,
$-\mathscr{E}_{k}^{visc}$
, related to
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}$
, while the source of energy is provided by the capillary stress and depicted by the dashed green line,
$-\mathscr{E}_{k}^{s}$
, related to
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{s}$
. The term
$-\mathscr{E}_{k}^{visc}$
is always a sink of energy, thus its behaviour is monotonic. In this experiment, the term
$-\mathscr{E}_{k}^{s}$
acts as a source of energy. Nevertheless, its behaviour is non-monotonic, meaning that the overall behaviour of the capillarity is not only a source but also a sink of energy.
Energy budget for case 1 2-D simulation. (a) Kinetic energy and its energy exchanges. (b) Bulk free energy and its energy exchange. (c) Interfacial free energy and its energy exchanges. Energies are depicted in solid lines and their energy exchanges in dashed lines.
In figure 10(b) we depict the bulk free energy (solid blue line,
$E_{{\it\phi}}$
) with its energy exchange term (dashed red line,
$-\mathscr{E}_{{\it\phi}}$
) related to
$\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}\boldsymbol{\cdot }\boldsymbol{j}$
. The bulk free energy is decreasing and its energy exchange shows a monotonic increase. The overall behaviour of this energy exchange is dissipative, i.e. it acts as a sink of energy. This means that the overall behaviour of the mass flux has a positive inner product with the steepest descent direction of the phase field almost everywhere. This confirms our theoretical prediction (see Remark 1).
In figure 10(c) we depict the interfacial free energy (solid blue line,
$E_{s}$
) with its energy exchange terms (dashed red line,
$-\mathscr{E}_{s}^{curv}$
) related to
$\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{j}$
and (dashed green line,
$-\mathscr{E}_{s}^{s}$
) related to
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}$
. The interfacial free energy decreases in a non-monotonic manner due to hydrodynamic effects, since the minimum surface area is not achieved monotonically. That is, the droplet oscillates around the circular configuration before settling at the steady configuration. A similar phenomenon is observed when an under-damped system oscillates around equilibrium. Both energy exchange terms appear as energy sinks. Regarding
$-\mathscr{E}_{s}^{curv}$
, this means that the Hessian of
${\it\phi}$
,
$\unicode[STIX]{x1D643}$
, and the mass flux gradient,
$\boldsymbol{{\rm\nabla}}\boldsymbol{j}$
, have a positive inner product. Finally, the term
$-\mathscr{E}_{s}^{s}$
has the same meaning as
$-(-\mathscr{E}_{k}^{s})$
.
The total energy should be constant. However, due to numerical errors and dissipation introduced by the time integrator, we observe a decrease of approximately 0.07 % in the total energy. Mass though is preserved exactly to machine precision.
9.2 Two-dimensional investigation: cases 2 and 3
In these examples, a gravitational field induces a density-driven flow and thus the droplets we consider in a liquid continuum rise and merge into a single bigger droplet.
9.2.1 Meniscus and bridge evolution: cases 2 and 3
Since this experimental set-up is beyond the assumptions of Eggers et al. (Reference Eggers, Lister and Stone1999), we fit the evolution of the bridge length with a function proportional to
$t^{b}$
. In figure 11 we depict the meniscus evolution defined on the isocurve
${\it\phi}=0.5$
during the rising and merging. Figure 12 shows the bridge length evolution, measured from the symmetry axis to the meniscus point. The scaling law obtained by Eggers et al. (Reference Eggers, Lister and Stone1999) (
${\sim}\sqrt{t}$
) does not hold for cases 2 and 3 due to buoyancy effects and different droplet radii. However, by using a least-squares fitting of the function
${\sim}t^{b}$
for the exponent
$b$
, we find that the bridge length can be approximated by
${\sim}t^{0.57}$
and
${\sim}t^{0.62}$
for cases 2 and 3, respectively. This suggests that the exponent
$b$
departs from 0.5 once buoyancy and non-uniform radii are considered. In case 1 the flow is driven by capillary effects, whereas in cases 2 and 3 the motion is mainly driven by buoyancy effects. This fact suggests that the flow may have different behaviours, i.e. different scaling laws.
Visualization of the meniscus and bridge evolution on
${\it\phi}=0.5$
. (a) Case 2, from
$t_{i}=4.12142$
to
$t_{f}=t_{i}+0.04$
. (b) Case 3, from
$t_{i}=4.11107$
to
$t_{f}=t_{i}+0.11$
. The isocurves are equally spaced over time.
Bridge length evolution
$r(t)$
(red lines), in (a) case 2 and (b) case 3. Curves are fitted with a polynomial function
$\bar{r}(t)\propto t^{b}$
.
9.2.2 General flow features in droplet dynamics: cases 2 and 3
The temporal evolution of the phase field depicted in figures 13 and 14 shows the phase-field dynamics during the merging and rising of droplets. In case 2, hydrodynamic effects lead to a variable interface thickness and a deformed shape as the droplets rise (cf. figure 13). The leading zone shows a narrower interface, whereas the trailing zone shows a thicker one. The droplet obtained after the merging shows a double-elliptic shape, the lower ellipse being flatter than the upper one. Additionally, in case 3, as the
$(Cn/We)^{-1}$
ratio increases (cf. figure 14), the surface tension dominates over the inertial effects, thus yielding a droplet which is almost spherical with constant interface thickness.
Phase-field time evolution: case 2. (a–e) Snapshots related to times
$t=0$
,
${\approx}4.09$
,
${\approx}4.19$
,
${\approx}4.83$
and
${\approx}8.70$
.
Phase-field time evolution: case 3. (a–e) Snapshots related to times
$t=0$
,
${\approx}4.07$
,
${\approx}4.10$
,
${\approx}4.25$
and
${\approx}9.31$
.
To characterize the hydrodynamics, the phase, vorticity, pressure and velocities are presented in figures 15 and 16 for cases 2 and 3, respectively. Relevant flow features are shown when the merging of the droplets starts. The highest vorticity occurs in the meniscus zone. Case 3, with the highest
$(Cn/We)^{-1}$
ratio, shows a vorticity that is seven times larger than that of case 2. The pressure jump across the interface of the biggest droplet is an order of magnitude higher as
$(Cn/We)^{-1}$
grows by an order of magnitude. This jump grows from
$10^{3}$
to
$10^{4}$
as the
$(Cn/We)^{-1}$
ratio grows from
$10^{2}$
to
$10^{3}$
. Likewise, the velocities increase as the
$(Cn/We)^{-1}$
ratio increases. The horizontal velocity,
$v_{1}$
, increases by a factor of five, while the vertical velocity,
$v_{2}$
, increases modestly by a factor of at most two.
Features of the flow: case 2. (a–e) Snapshots of the phase field
${\it\phi}$
, vorticity
$(1/2)\boldsymbol{{\rm\nabla}}\times \boldsymbol{v}$
, pressure
$p$
, horizontal
$v_{1}$
and vertical
$v_{2}$
velocities at
$t\approx 4.09$
.
Features of the flow: case 3. (a–e) Snapshots of the phase field
${\it\phi}$
, vorticity
$(1/2)\boldsymbol{{\rm\nabla}}\times \boldsymbol{v}$
, pressure
$p$
, horizontal
$v_{1}$
and vertical
$v_{2}$
velocities at
$t\approx 4.07$
.
The droplet merging process is essentially modified (from case 2 to case 3) as the
$(Cn/We)^{-1}$
ratio increases by an order of magnitude. In case 2, the smaller (bottom) droplet pushes the larger (top) droplet, which accompanies the motion. In case 3, the droplets push against each other at the time of merging. The vertical velocities evidence these phenomena. In case 2, the vertical velocities inside the droplets are positive, whereas in case 3, the highest positive vertical velocity occurs in the northern hemisphere of the smaller (bottom) droplet and the highest negative vertical velocity is encountered at the southern hemisphere of the larger (top) one.
Stress profiles depicted in the last two snapshots of figures 17 and 18 for cases 2 and 3, respectively, show the cross-section integrals
$\int T_{22}^{visc}\,\text{d}A_{2}$
and
$\int T_{22}^{s}\,\text{d}A_{2}$
along the vertical axis
$x_{2}$
, with
$\text{d}A_{2}=\text{d}x_{1}\text{d}x_{3}$
. The capillary stress is higher at the interface, having its peak value at the leading zone during the merging and rising. The magnitude of the capillary stress is roughly five times higher for the highest
$(Cn/We)^{-1}$
ratio, in case 3, also showing a smaller difference in magnitude between the leading and trailing zones after merging (this difference decreases as the droplet takes a constant interface thickness). During the merging, the viscous stress is higher at the merging region for both cases. As the
$(Cn/We)^{-1}$
ratio increases, in case 3, the viscous stress increases by a factor of two, whereas during the rising process this difference is five times. During the rising, the resultant viscous stress takes positive values in the leading zone and negative ones in the trailing zone of the droplet. The local viscous stress is negative in the northern hemisphere of the droplet, but the viscous stress that lies outside of the droplet is positive and higher in magnitude. This is due to the higher viscosity of the interstitial fluid, which results in a positive viscous stress profile. The opposite behaviour is observed in the southern hemisphere, which has a negative resultant viscous stress profile.
9.2.3 Energy exchange in the merging and rising of droplets: cases 2 and 3
During the merging of droplets, the highest intensive kinetic energy is located where the curvature of
${\it\phi}$
is highest, i.e. in the meniscus region (see figures 17 and 18 for cases 2 and 3, respectively). The source of this kinetic energy is given by the capillary stress. Evidence for this is taken from the second snapshot depicted in figures 19 and 20 for cases 2 and 3, respectively, which correspond to the energy exchange
${\it\epsilon}_{k}^{s}$
. It is the source/sink given by the capillary stress
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{s}$
acting as a source of kinetic energy in the meniscus zone. Although the highest
$(Cn/We)^{-1}$
ratio, in case 3, shows the highest kinetic energy (approximately 20 times greater than that observed for a lower
$(Cn/We)^{-1}$
ratio, case 2) after the merging, the resultant droplets in both cases have similar values of intensive kinetic energy. This intensity decreases by one and two orders of magnitude after the merging phenomenon, for cases 2 and 3, respectively (cf. the second snapshot in figures 17 and 18).
Intensive energies and stress profiles at (a–e)
$t\approx 4.09$
and (f–j)
$t\approx 8.70$
for case 2. From left to right, the panels in each row show the intensive kinetic energy
$e_{k}$
, bulk free energy
$e_{{\it\phi}}$
, interfacial free energy
$e_{s}$
, capillary stress profile
$T_{22}^{visc}$
and viscous stress profile
$T_{22}^{s}$
(
$\text{d}A_{2}=\text{d}x_{1}\text{d}x_{3}$
).
Intensive energies and stress profiles at (a–e)
$t\approx 4.07$
and (f–j)
$t\approx 9.31$
for case 3. From left to right, the panels in each row show the intensive kinetic energy
$e_{k}$
, bulk free energy
$e_{{\it\phi}}$
, interfacial free energy
$e_{s}$
, capillary stress profile
$T_{22}^{visc}$
and viscous stress profile
$T_{22}^{s}$
(
$\text{d}A_{2}=\text{d}x_{1}\text{d}x_{3}$
).
The viscous sink of kinetic energy
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}$
, i.e. the energy exchange
${\it\epsilon}_{k}^{visc}$
(first snapshot depicted in figure 19 and figure 20 for cases 2 and 3, respectively) is also high in the confluence and meniscus zones. Increasing the
$(Cn/We)^{-1}$
ratio by an order of magnitude (case 2
$\Rightarrow$
case 3) yields a viscous sink greater by two orders of magnitude. Nonetheless, just before and after the merging, the region of highest viscous dissipation is located outside the droplet, in the interstitial fluid, near the equatorial region of the droplet (cf. figure 21), as is usual in a flow around a sphere, with the same intensity in both cases. The capillary source/sink of kinetic energy
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{s}$
(second snapshot depicted in figures 19 and 20 for cases 2 and 3, respectively) acts as a sink of energy in the confluence region while in the meniscus region it acts as a source of energy. During the rising of the resultant droplet in the case of the lowest
$(Cn/We)^{-1}$
ratio, the capillary source of kinetic energy is located along the leading zone, whereas capillarity acts as a sink of kinetic energy on the trailing zone. However, in the case of the highest
$(Cn/We)^{-1}$
ratio, at the leading zone the capillarity acts as a source (sink) of kinetic energy just inside (outside) the droplet. At the trailing zone, the opposite behaviour is observed.
In the merging of droplets, the intensity of the bulk free energy is greater by an order of magnitude for the lowest
$(Cn/We)^{-1}$
ratio, case 2. After the merging, the maximum and minimum values are equal for both cases. Meanwhile, the bulk free energy experiences a reduction by an order of magnitude for case 2 (lowest
$(Cn/We)^{-1}$
ratio), whereas for case 3 (highest
$(Cn/We)^{-1}$
ratio) the maximum and minimum values remain the same during the merging and rising. The bulk free energy density shows its highest value at the interface, whereas its lowest values define a hydrodynamic wake below the droplets as they rise (see the third snapshot in figures 17 and 18 for cases 2 and 3, respectively).
The source/sink of bulk free energy occurs mainly at the interface (see the third snapshot in figures 19 and 20 for cases 2 and 3, respectively). However, there is a region of energy production contouring the trailing zone that goes down towards the wake. This region arises due to hydrodynamic effects. In case 2 (lowest
$(Cn/We)^{-1}$
ratio) the bulk free energy exchange is a source and a sink at the trailing and leading zones, respectively. Nevertheless, in case 3 (highest
$(Cn/We)^{-1}$
ratio) the bulk free energy exchange is a sink at the interface.
As expected, the interfacial free energy is concentrated along the interface. During the merging, when the interfaces are coalescing, the interfacial energy decreases. In addition, at that time and place, the highest source of energy is provided by the first term in the interfacial energy exchange,
${\it\epsilon}_{s}^{curv}\propto \unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{j}$
. The second term has the same effect as the capillary power for kinetic energy exchanges with the opposite sign.
To understand the overall instantaneous behaviour of the energy exchanges, figure 21 uses blue to indicate regions of energy sources and red to indicate regions of energy sinks.
Intensive energy exchanges at (a–e)
$t\approx 4.09$
and (f–j)
$t\approx 8.70$
for case 2. From left to right, the panels in each row show the intensive energy exchanges
${\it\epsilon}_{k}^{visc}$
,
${\it\epsilon}_{k}^{s}$
,
${\it\epsilon}_{{\it\phi}}$
,
${\it\epsilon}_{s}^{curv}$
and
${\it\epsilon}_{s}^{s}$
.
Intensive energy exchanges at (a–e)
$t\approx 4.07$
and (f–j)
$t\approx 9.31$
for case 3. From left to right, the panels in each row show the intensive energy exchanges
${\it\epsilon}_{k}^{visc}$
,
${\it\epsilon}_{k}^{s}$
,
${\it\epsilon}_{{\it\phi}}$
,
${\it\epsilon}_{s}^{curv}$
and
${\it\epsilon}_{s}^{s}$
.
Intensive energy exchanges, source (blue regions) and sink (red regions) of energy: (a) case 2; (b) case 3.
Energy budget: (a,e) potential energy and its energy exchange. (b,f) Kinetic energy and its energy exchanges. (c,g) Bulk free energy and its energy exchange. (d,h) Interfacial free energy and its energy exchanges, for cases 2 and 3, respectively. Energies are depicted in solid lines and their energy exchanges in dashed lines.
9.2.4 Energy budget of the flow: cases 2 and 3
Figure 22 presents the extensive energies and their related energy exchanges over time for cases 2 and 3. In figure 22(a,e) the potential energy and its energy exchanges are shown. The potential energy (solid blue line,
$E_{p}$
) shows a linear decay. Its first term (dashed red line,
$-\mathscr{E}_{p}^{mass}$
) related to
$\boldsymbol{e}_{2}\boldsymbol{\cdot }\boldsymbol{j}$
in the energy exchange is almost zero. The second term (dashed green line,
$-\mathscr{E}_{p}^{buoy}$
) related to
${\it\phi}\,v_{2}$
in the potential energy exchange works as a sink of energy, but it cancels out with the third term of the kinetic energy exchange,
$-\mathscr{E}_{k}^{buoy}$
.
In figure 22(b,f) the kinetic energy and its energy exchanges are shown. During the merging, the kinetic energy (solid blue line,
$E_{k}$
) shows a smooth (sharp) peak in case 2 (case 3) which is quickly damped. The energy exchange terms are depicted using dashed lines. The sink of energy done by the viscous stress is depicted by the dashed red line,
$-\mathscr{E}_{k}^{visc}$
, related to
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}$
, while the sink (source) of energy in case 2 (case 3) is provided by the capillary stress and depicted by the dashed green line,
$-\mathscr{E}_{k}^{s}$
, related to
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{s}$
. The overall contribution
$-\mathscr{E}_{k}^{s}$
changes its behaviour from a sink to a source of energy when the
$(Cn/We)^{-1}$
ratio increases. Nevertheless, in both cases the kinetic energy peak is provided by the energy exchange
$-\mathscr{E}_{k}^{s}$
.
In figure 22(c,g), the bulk free energy (solid blue line,
$E_{{\it\phi}}$
) is depicted with its energy exchange term (dashed red line,
$-\mathscr{E}_{{\it\phi}}$
) related to
$\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}\boldsymbol{\cdot }\boldsymbol{j}$
. In all examples presented here, the overall behaviour of these energies is quite similar: the bulk free energy has a monotonic decay and its energy exchange acts as a sink of energy.
In figure 22(d,h), the interfacial free energy (solid blue line,
$E_{s}$
) is depicted with its energy exchange terms (dashed red line,
$-\mathscr{E}_{s}^{curv}$
) related to
$\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{j}$
and (dashed green line,
$-\mathscr{E}_{s}^{s}$
) related to
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}$
. The interfacial free energy experiences a sudden drop during the merging. This drop increases as the
$(Cn/We)^{-1}$
ratio increases. Regarding
$-\mathscr{E}_{s}^{curv}$
, this means that the Hessian of
${\it\phi}$
,
$\unicode[STIX]{x1D643}$
, and the mass flux gradient,
$\boldsymbol{{\rm\nabla}}\boldsymbol{j}$
, have a positive inner product considering the overall behaviour. Finally, the term
$-\mathscr{E}_{s}^{s}$
has the same meaning as
$-(-\mathscr{E}_{k}^{s})$
.
(a)
$Q$
-criterion evolution. (b) Phase-field evolution.
Energy budget: case 4 3-D simulation. (a) Potential energy and its energy exchange. (b) Kinetic energy and its energy exchanges. (c) Bulk free energy and its energy exchange. (d) Interfacial free energy and its energy exchanges. Energies are depicted in solid lines and their energy exchanges in dashed lines.
Due to numerical errors and dissipation introduced by the time integrator, we observe a decrease of less than 3.8 % and 3.6 % in the total energy, for cases 2 and 3, respectively. As in case 1, mass is preserved exactly to machine precision.
9.3 Three-dimensional investigation: case 4
In the 3-D case, the flow is also driven by buoyancy effects. The droplet is lighter and less viscous than the interstitial fluid.
9.3.1 General flow features in droplet dynamics: case 4
The
$Q$
-criterion is defined by the scalar
$Q=1/2(\unicode[STIX]{x1D652}\boldsymbol{ : }\unicode[STIX]{x1D652}-\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D63F})$
, where
$\unicode[STIX]{x1D63F}$
is the symmetric part of the velocity gradient, i.e. the strain rate tensor, and
$\unicode[STIX]{x1D652}$
the anti-symmetric part of the velocity gradient, i.e. the spin rate tensor. Figure 23 depicts at the top two isosurfaces of the
$Q$
-criterion, i.e.
$Q=-0.05$
and 0.05. Positive (negative) values show regions where the spin (strain) rate overcomes the strain (spin) rate.
9.3.2 Energy budget of the flow: case 4
Figure 24 presents the extensive energies and their related energy exchanges over time for this 3-D simulation, case 4. In figure 24(a) the potential energy and its energy exchanges are shown. The potential energy (solid blue line,
$E_{p}$
) shows a linear decay until the droplet merges in a thin film of fluid above the interstitial one, while its first term (dashed red line,
$-\mathscr{E}_{p}^{mass}$
), related to
$\boldsymbol{e}_{2}\boldsymbol{\cdot }\boldsymbol{j}$
in the energy exchange, is almost zero. The second term (dashed green line,
$-\mathscr{E}_{p}^{buoy}$
), related to
${\it\phi}v_{2}$
, works in the potential energy exchange as a sink of energy, but it cancels out with the third term of the kinetic energy exchange
$-\mathscr{E}_{k}^{buoy}$
, as mentioned in cases 2 and 3.
In figure 24(b) the kinetic energy and its energy exchanges are shown. Before the merging between the droplet and the thin film of fluid, the kinetic energy (solid blue line,
$E_{k}$
) grows. When the droplet is close to the thin film the kinetic energy decreases. During the merging, the kinetic energy shows a smooth peak which is quickly damped by viscous effects. The energy exchange terms are depicted using dashed lines. The sink of energy contribution done by the viscous stress is depicted by the dashed red line,
$-\mathscr{E}_{k}^{visc}$
, related to
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}$
, while the source/sink of energy provided by the capillary stress is depicted by the dashed green line,
$-\mathscr{E}_{k}^{s}$
, related to
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{s}$
. Before the merging process, capillary acts as a sink of energy. During the merging, capillary changes behaviour to act as a source of energy.
In figure 24(c) the bulk free energy (solid blue line,
$E_{{\it\phi}}$
) is depicted with its energy exchange term (dashed red line,
$-\mathscr{E}_{{\it\phi}}$
) related to
$\boldsymbol{{\rm\nabla}}{\it\eta}_{{\it\phi}}\boldsymbol{\cdot }\boldsymbol{j}$
. In all examples presented here, the overall behaviour of these energies is quite similar.
In figure 24(d) the interfacial free energy (solid blue line,
$E_{s}$
) is depicted with its energy exchange terms (dashed red line,
$-\mathscr{E}_{s}^{curv}$
) related to
$\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{j}$
and (dashed green line,
$-\mathscr{E}_{s}^{s}$
) related to
$\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}$
. Before the merging, the interfacial free energy grows to a stable baseline. However, the interfacial free energy experiences a sudden drop during the merging. The energy exchange term
$-\mathscr{E}_{s}^{curv}$
acts as a source of energy up to the merging, and changes its behaviour during the merging to act as a sink of energy. A similar behaviour is observed in the second term of the energy exchange,
$-\mathscr{E}_{s}^{s}$
. Finally, the term
$-\mathscr{E}_{s}^{s}$
has the same meaning as
$-(-\mathscr{E}_{k}^{s})$
.
Due to numerical errors and dissipation introduced by the time integrator, we observe a decrease of less than 0.7 % in the total energy. The error in the mass conservation is 0.002 %.
10 Conclusions
We develop the energy budget equation of the coupled NSCH flow, the chief result of this work. We derive two new dimensionless numbers and denote these by
$Lm$
and
$Ln$
, which relate the critical free energy density to the capillary energy and the kinetic energy at the molecular scale to the critical free energy, respectively. In addition, we describe how the mass flux behaves in simple configurations, in one and two dimensions. We perform highly resolved simulations, showing the energy exchanges and the energy budget of the flow. As a particular result, we show how the merging phenomenon of droplets is essentially modified when the
$(Cn/We)^{-1}$
ratio increases. To achieve such a detailed reproduction of the physical features of these multiphase flows we relied on robust and efficient (in-house and open-source) software tools, such as PetIGA and PetIGA-MF, to implement a robust discretization based on isogeometric analysis. We describe and analyse all energetic interactions at the discrete level as direct counterparts of those derived at the continuous level. We are able to achieve this by using divergence-conforming discretizations for the velocity–pressure pair and a mixed form for the phase field and chemical potential. Finally, we show that our numerical simulations are in good agreement with the available analytical predictions and experimental results, concluding that employing NSCH equations to model droplet dynamics is a sensible approach worthy of further research.
Acknowledgements
This work is part of the European Union’s Horizon 2020 research and innovation programme of the Marie Skłodowska-Curie grant agreement no. 644602.2014-0191. This publication also was made possible by a National Priorities Research Program grant NPRP 7-1482-1-278 from the Qatar National Research Fund (a member of The Qatar Foundation) and Agencia Nacional de Promoción Científica y Tecnológica (Argentina) grants PICT 2014–2660 and PICT-E 2014–0191.
Supplementary movies
Supplementary movies are available at http://dx.doi.org/10.1017/jfm.2016.277.
Appendix A. First and second laws of thermodynamics
A.1 First law
Here, we develop the first law in the classical manner, i.e. directly from the governing equations of momenta and mass transfer, cf. Guo & Lin (Reference Guo and Lin2015). The inner product between the material derivative of the velocity and the velocity yields the balance between the external and internal mechanical rates of work done by the stress tensor. The product between the material derivative of the phase field and the chemical potential yields the balance between the external and internal chemical rates of work done by the mass flux, including buoyancy effects.
We look for the following form:
$$\begin{eqnarray}\displaystyle & {\dot{e}}_{k}+{\dot{w}}_{m}^{int}={\dot{w}}_{m}^{ext}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \partial _{{\it\phi}}^{}{\it\psi}\dot{{\it\phi}}+\partial _{\boldsymbol{{\rm\nabla}}{\it\phi}}^{}{\it\psi}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}{\it\phi})^{\boldsymbol{\cdot }}+{\dot{w}}_{c}^{int}={\dot{w}}_{c}^{ext}, & \displaystyle\end{eqnarray}$$
${\dot{w}}_{m}^{int}$
and
${\dot{w}}_{m}^{ext}$
(
${\dot{w}}_{c}^{int}$
and
${\dot{w}}_{c}^{ext}$
) account for the internal and external mechanical power (internal and external chemical power), respectively. Here,
$\partial _{{\it\phi}}^{}$
is the partial derivative with respect to
${\it\phi}$
and
$\partial _{\boldsymbol{{\rm\nabla}}{\it\phi}}^{}$
is the partial derivative with respect to
$\boldsymbol{{\rm\nabla}}{\it\phi}$
. Consider the free energy in the form (3.4), i.e.
${\it\psi}={\it\psi}_{{\it\phi}}+{\it\psi}_{s}+e_{p}$
. In the general case
$\dot{{\it\psi}}=\partial _{{\it\phi}}^{}{\it\psi}\,\dot{{\it\phi}}+\partial _{\boldsymbol{{\rm\nabla}}{\it\phi}}^{}{\it\psi}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}{\it\phi})^{\boldsymbol{\cdot }}+\partial _{{\it\theta}}^{}{\it\psi}\,\dot{{\it\theta}}$
, where
$\partial _{{\it\theta}}^{}$
is the partial derivative with respect to
${\it\theta}$
, whereas in the absence of heat transfer
$\dot{{\it\psi}}=\partial _{{\it\phi}}^{}{\it\psi}\,\dot{{\it\phi}}+\partial _{\boldsymbol{{\rm\nabla}}{\it\phi}}^{}{\it\psi}\cdot (\boldsymbol{{\rm\nabla}}{\it\phi})^{\boldsymbol{\cdot }}$
. Thus, (A 1b
) is
$\dot{{\it\psi}}+{\dot{w}}_{c}^{int}={\dot{w}}_{c}^{ext}$
.
In the dimensional form, the internal and external powers are
$$\begin{eqnarray}\displaystyle & {\dot{w}}_{m}^{int}=\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\dot{w}}_{m}^{ext}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D64F})-{\it\phi}\unicode[STIX]{x27E6}{\it\rho}\unicode[STIX]{x27E7}gv_{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\dot{w}}_{c}^{int}=-\boldsymbol{{\rm\nabla}}{\it\eta}\boldsymbol{\cdot }\boldsymbol{j}+{\it\gamma}_{{\it\phi}}\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\dot{w}}_{c}^{ext}=-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\eta}\,\boldsymbol{j}-{\it\gamma}_{{\it\phi}}\dot{{\it\phi}}\boldsymbol{{\rm\nabla}}{\it\phi})+{\it\phi}\unicode[STIX]{x27E6}{\it\rho}\unicode[STIX]{x27E7}gv_{2}; & \displaystyle\end{eqnarray}$$
$j$
in the right-hand side in both (2.1b
) and (A 2d
).
The first law of thermodynamics represents an energy balance and states the interplay between the kinetic energy
$e_{k}$
, the internal energy
$e_{i}$
, the rate at which (mechanical and chemical) power is expended and the rate at which energy in the form of heat is transferred, i.e.
$$\begin{eqnarray}\dot{e_{T}}={\dot{e}}_{k}+{\dot{e}}_{i}={\dot{w}}_{m}^{ext}+{\dot{w}}_{c}^{ext}-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{q}+q,\end{eqnarray}$$
where
$\boldsymbol{q}$
is the heat flux and
$q$
is a heat sink/source.
Finally, using (A 2), we obtain the first law of thermodynamics
$$\begin{eqnarray}{\dot{e}}_{i}=\partial _{{\it\phi}}^{}{\it\psi}\dot{{\it\phi}}+\partial _{\boldsymbol{{\rm\nabla}}{\it\phi}}^{}{\it\psi}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}{\it\phi})^{\boldsymbol{\cdot }}+\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}-\boldsymbol{{\rm\nabla}}{\it\eta}\boldsymbol{\cdot }\boldsymbol{j}-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{q}+q.\end{eqnarray}$$
In the absence of heat transfer the first thermodynamic law is
$$\begin{eqnarray}{\dot{e}}_{i}=\dot{{\it\psi}}+\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}-\boldsymbol{{\rm\nabla}}{\it\eta}\boldsymbol{\cdot }\boldsymbol{j}.\end{eqnarray}$$
A.2 Second law
The second law of thermodynamics (in the form of the Clausius–Duhem inequality or entropy imbalance) states that the entropy
$s$
should grow at least with a rate given by the entropy flux
$\boldsymbol{q}/{\it\theta}$
added to the entropy supply
$q/{\it\theta}$
, i.e.
$$\begin{eqnarray}{\dot{s}}\geqslant -\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left({\displaystyle \frac{\boldsymbol{q}}{{\it\theta}}}\right)+{\displaystyle \frac{q}{{\it\theta}}}={\displaystyle \frac{1}{{\it\theta}}}\left(-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{q}+{\displaystyle \frac{1}{{\it\theta}}}\boldsymbol{{\rm\nabla}}{\it\theta}\boldsymbol{\cdot }\boldsymbol{q}+q\right),\end{eqnarray}$$
where
$s$
is the entropy. By definition, the free energy is
$$\begin{eqnarray}{\it\psi}=e_{i}-{\it\theta}s.\end{eqnarray}$$
Taking the material derivative, we obtain
$$\begin{eqnarray}\dot{{\it\psi}}={\dot{e}}_{i}-\dot{{\it\theta}}s-{\it\theta}{\dot{s}}.\end{eqnarray}$$
Substituting the first law (A 4) and (A 8) into (A 6), and considering that the entropy is
$s=-\partial _{{\it\theta}}^{}{\it\psi}$
, the second law of thermodynamics is obtained in the form of the entropy imbalance, i.e.
$$\begin{eqnarray}{\dot{s}}={\displaystyle \frac{1}{{\it\theta}}}\left(\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D64F}^{visc}-\boldsymbol{{\rm\nabla}}{\it\eta}\boldsymbol{\cdot }\boldsymbol{j}-{\displaystyle \frac{1}{{\it\theta}}}\boldsymbol{{\rm\nabla}}{\it\theta}\boldsymbol{\cdot }\boldsymbol{q}\right)\geqslant 0.\end{eqnarray}$$
In the absence of heat transfer, the entropy imbalance yields
$$\begin{eqnarray}2{\it\mu}({\it\phi})\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D63F}+{\it\alpha}({\it\phi})\boldsymbol{{\rm\nabla}}{\it\eta}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\eta}\geqslant 0.\end{eqnarray}$$
Finally, this shows that our model guarantees the entropy production if
${\it\mu}({\it\phi}),{\it\alpha}({\it\phi})\geqslant 0$
.
Appendix B. Identities used in the energy budget
To obtain the energy dissipation related to each energy, we employ some identities coupled with the governing equations as well as the constitutive relations. We include the ones used in our derivations here. It is worth noting that all identities have been derived considering that both second-order tensors
$\unicode[STIX]{x1D63F}$
and
$\unicode[STIX]{x1D643}$
are symmetric.
Finally, in the kinetic energy
$$\begin{eqnarray}\displaystyle & \left({\textstyle \frac{1}{2}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{v}\right)^{\boldsymbol{\cdot }}=\boldsymbol{v}\boldsymbol{\cdot }\dot{\boldsymbol{v}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\boldsymbol{v})=(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D63F})\boldsymbol{\cdot }\boldsymbol{v}+\unicode[STIX]{x1D63F}\boldsymbol{ : }\unicode[STIX]{x1D63F}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{v})=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi})\boldsymbol{\cdot }\boldsymbol{v}+\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \dot{{\it\psi}}_{{\it\phi}}({\it\phi})={\it\psi}_{{\it\phi}}^{\prime }\dot{{\it\phi}}={\it\eta}_{{\it\phi}}\dot{{\it\phi}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(f({\it\phi})\unicode[STIX]{x1D643}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi})=f^{\prime }({\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}+f({\it\phi})(\unicode[STIX]{x1D643}\boldsymbol{ : }\unicode[STIX]{x1D643}+(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}). & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \left({\textstyle \frac{1}{2}}\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}\right)^{\boldsymbol{\cdot }}=\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\dot{{\it\phi}}-\unicode[STIX]{x1D63F}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}{\it\phi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\dot{{\it\phi}}=(\boldsymbol{{\rm\nabla}}{\it\phi})^{\boldsymbol{\cdot }}+\boldsymbol{{\rm\nabla}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{j})=\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{j}+\boldsymbol{{\rm\nabla}}{\it\phi}\boldsymbol{\cdot }{\rm\Delta}\boldsymbol{j}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(f({\it\phi})\unicode[STIX]{x1D643}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})) & = & \displaystyle f^{\prime }({\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}\nonumber\\ \displaystyle & & \displaystyle +\;\,f({\it\phi})[(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})+\unicode[STIX]{x1D643}\boldsymbol{ : }(\boldsymbol{{\rm\nabla}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643})].\end{eqnarray}$$
In the potential energy
$$\begin{eqnarray}\displaystyle & ({\it\phi}x_{2})^{\boldsymbol{\cdot }}=x_{2}\dot{{\it\phi}}+{\it\phi}v_{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(x_{2}\,\boldsymbol{j})=(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j})x_{2}+\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{e}_{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(f({\it\phi})\boldsymbol{e}_{2}\boldsymbol{\cdot }\unicode[STIX]{x1D643})=f^{\prime }({\it\phi})\unicode[STIX]{x1D643}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\it\phi}\otimes \boldsymbol{e}_{2}+f({\it\phi})\boldsymbol{e}_{2}\boldsymbol{\cdot }(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\unicode[STIX]{x1D643}). & \displaystyle\end{eqnarray}$$
'/%3e%3cpath%20d='M17%207.971L15.711%206.677%2010.912%2011.495%2010.912%200%209.088%200%209.088%2011.495%204.289%206.677%203%207.971%2010%2015z'%20transform='translate(8%208)'/%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
