On a thermodynamic framework for developing boundary conditions for Korteweg fluids

TL;DR

This paper develops thermodynamically consistent boundary conditions for Korteweg fluids using a membrane-based approach, generalizing slip and contact angle conditions, with numerical validation for contact angle hysteresis.

Contribution

It introduces a new thermodynamic framework for deriving boundary conditions for Korteweg fluids, ensuring compatibility with the second law and encompassing contact angle phenomena.

Findings

Derived boundary conditions compatible with thermodynamics

Generalized Navier slip and contact angle conditions

Numerical demonstration of contact angle hysteresis

Abstract

We provide a derivation of several classes of boundary conditions for fluids of Korteweg-type using a simple and transparent thermodynamic approach that automatically guarentees that the derived boundary conditions are compatible with the second law of thermodynamics. The starting assumption of our approach is to describe the boundary of the domain as the membrane separating two different continua, one inside the domain, and the other outside the domain. With this viewpoint one may employ the framework of continuum thermodynamics involving singular surfaces. This approach allows us to identify, for various classes of surface Helmholtz free energies, the corresponding surface entropy production mechanisms. By establishing the constitutive relations that guarantee that the surface entropy production is non-negative, we identify a new class of boundary conditions, which on one hand generalizes in a nontrivial manner the Navier's slip boundary conditions, and on the other hand describes dynamic and static contact angle conditions. We explore the general model in detail for a particular case of Korteweg fluid where the Helmholtz free energy in the bulk is that of a van der Waals fluid. We perform a series of numerical experiments to document the basic qualitative features of the novel boundary conditions and their practical applicability to model phenomena such as the contact angle hysteresis.

Equations

$$\displaystyle\frac{d}{dt}\mbox{\boldmath$\Psi$}(V)=-\mathcal{F}^{\mbox{\tiny{$\Psi$}}}(V)+\mathcal{P}^{\mbox{\tiny{$\Psi$}}}(V)+\mathcal{S}^{\mbox{\tiny{$\Psi$}}}(V)\ .$$

$$\displaystyle\mbox{\boldmath$\Psi$}(V)=\int_{V^{+}\cup V^{-}}\Psi_{V}\,dx+\int_{\Gamma}\Psi_{\Gamma}\,dS\ ,$$

$$\displaystyle\mathcal{F}^{\mbox{\tiny{$\Psi$}}}(V)=\int_{S^{+}\cup S^{-}}\mbox{\boldmath$\Phi$}_{V}^{\mbox{\tiny{$\Psi$}}}\cdot\mathrm{{\bf n}}\,dS+\int_{\partial\Gamma}\mbox{\boldmath$\Phi$}_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\Psi$}}}\cdot{\mbox{\boldmath$\nu$}}\,dl\ ,$$

$$\displaystyle\mathcal{P}^{\mbox{\tiny{$\Psi$}}}(V)=\int_{V^{+}\cup V^{-}}\Pi_{V}^{\mbox{\tiny{$\Psi$}}}\,dx+\int_{\mbox{\scriptsize{$\Gamma$}}}\Pi_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\Psi$}}}\,dS\ ,$$

$$\displaystyle\mathcal{S}^{\mbox{\tiny{$\Psi$}}}(V)=\int_{V^{+}\cup V^{-}}\Sigma_{V}^{\mbox{\tiny{$\Psi$}}}dx+\int_{\Gamma}\Sigma_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\Psi$}}}\,dS\ .$$

$$\frac{\partial\Psi_{V}}{\partial t}+\hskip 0.85358pt\mathrm{div}\left(\mbox{\boldmath$\Phi$}_{V}^{\mbox{\tiny{$\Psi$}}}+\Psi_{V}{\bf v}\right)-\Pi_{V}^{\mbox{\tiny{$\Psi$}}}-\Sigma_{V}^{\mbox{\tiny{$\Psi$}}}=0\ ,$$

$$\frac{D_{\mbox{\scriptsize{$\Gamma$}}}\Psi_{\mbox{\scriptsize{$\Gamma$}}}}{Dt}+\Psi_{\mbox{\scriptsize{$\Gamma$}}}\left(\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\tau}-2K_{M}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\mathrm{n}}\right)+\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}\mbox{\boldmath$\Phi$}_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\Psi$}}}-\Pi_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\Psi$}}}-\Sigma_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\Psi$}}}=-\llbracket\mbox{\boldmath$\Phi$}_{V}^{\mbox{\tiny{$\Psi$}}}+\Psi_{V}({\bf v}-{\bf v}_{\mbox{\scriptsize{$\Gamma$}}})\rrbracket\cdot\mathrm{{\bf n}}_{\mbox{\scriptsize{$\Gamma$}}}\ ,$$

$$\displaystyle\frac{D_{\mbox{\scriptsize{$\Gamma$}}}A}{Dt}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left.\frac{\partial A}{\partial t}\right|_{{\bf X}_{\mbox{\scriptsize{$\Gamma$}}}}\ ,$$

$$\displaystyle{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\tau}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(\mathbb{I}{-}\mathrm{{\bf n}}_{\mbox{\scriptsize{$\Gamma$}}}{\otimes}\mathrm{{\bf n}}_{\mbox{\scriptsize{$\Gamma$}}}){\bf v}_{\mbox{\scriptsize{$\Gamma$}}}\ ,\hskip 28.45274pt{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\mathrm{n}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}}}{\cdot}\mathrm{{\bf n}}_{\mbox{\scriptsize{$\Gamma$}}}\ ,$$

$$\frac{\partial \rho }{\partial t}+\mathrm{div}(\rho \mathbf{v})=0,$$

$$\frac{\partial (\rho \mathbf{v})}{\partial t}+\mathrm{div}(\rho \mathbf{v}\otimes \mathbf{v})=\mathrm{div}\mathbb{T}+\rho \mathbf{b},$$

$$\mathbb{T}={\mathbb{T}}^{\mathrm{T}},$$

$$\frac{\partial \left(\rho (e+\frac{1}{2}|\mathbf{v}{|}^2)\right)}{\partial t}+\mathrm{div}\left(\rho (e+\frac{1}{2}|\mathbf{v}{|}^2)\mathbf{v}\right)=-\mathrm{div}\mathbf{q}+\mathrm{div}\left(\mathbb{T}\mathbf{v}\right)+\rho r,$$

$$\rho \frac{De}{Dt}=-\mathrm{div}\mathbf{q}+\mathbb{T}:\mathbb{D}+\rho r,$$

$$\displaystyle\frac{\partial(\rho\eta)}{\partial t}+\hskip 0.85358pt\mathrm{div}\left(\rho\eta{\bf v}+{\mbox{\boldmath$\Phi$}^{\mbox{\tiny{$\eta$}}}}\right)-{\Sigma^{\mbox{\tiny{$\eta$}}}}={\Pi^{\mbox{\tiny{$\eta$}}}}\hskip 28.45274pt\text{with}\hskip 28.45274pt{\Pi^{\mbox{\tiny{$\eta$}}}}\geq 0\ ,$$

$$\frac{D_{\mbox{\scriptsize{$\Gamma$}}}\rho_{\mbox{\scriptsize{$\Gamma$}}}}{Dt}+\rho_{\mbox{\scriptsize{$\Gamma$}}}(\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\tau}-2K_{M}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\mathrm{n}})=-\llbracket\rho({\bf v}-{\bf v}_{\mbox{\scriptsize{$\Gamma$}}})\rrbracket\cdot\mathrm{{\bf n}}_{\mbox{\scriptsize{$\Gamma$}}}\ ,$$

$$\frac{D_{\mbox{\scriptsize{$\Gamma$}}}(\rho_{\mbox{\scriptsize{$\Gamma$}}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}}})}{Dt}+\rho_{\mbox{\scriptsize{$\Gamma$}}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}}}(\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\tau}-2K_{M}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\mathrm{n}})-\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}\mathbb{T}_{\mbox{\scriptsize{$\Gamma$}}}-\rho_{\mbox{\scriptsize{$\Gamma$}}}{\bf b}_{\mbox{\scriptsize{$\Gamma$}}}=-\llbracket\rho{\bf v}\otimes({\bf v}{-}{\bf v}_{\mbox{\scriptsize{$\Gamma$}}})-\mathbb{T}\rrbracket\cdot\mathrm{{\bf n}}_{\mbox{\scriptsize{$\Gamma$}}}\ ,$$

$$\mathbb{T}_{\mbox{\scriptsize{$\Gamma$}}}=\mathbb{T}_{\mbox{\scriptsize{$\Gamma$}}}^{\mathrm{T}}\ ,$$

$$\displaystyle\frac{D_{\mbox{\scriptsize{$\Gamma$}}}\left(\rho_{\mbox{\scriptsize{$\Gamma$}}}(e_{\mbox{\scriptsize{$\Gamma$}}}+\frac{1}{2}|{\bf v}_{\mbox{\scriptsize{$\Gamma$}}}|^{2})\right)}{Dt}+\rho_{\mbox{\scriptsize{$\Gamma$}}}(e_{\mbox{\scriptsize{$\Gamma$}}}+\frac{1}{2}|{\bf v}_{\mbox{\scriptsize{$\Gamma$}}}|^{2})(\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\tau}-2K_{M}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\mathrm{n}})=-\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}{\bf q}_{\mbox{\scriptsize{$\Gamma$}}}+\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}\left(\mathbb{T}_{\mbox{\scriptsize{$\Gamma$}}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}}}\right)$$

$$\displaystyle+\rho_{\mbox{\scriptsize{$\Gamma$}}}r_{\mbox{\scriptsize{$\Gamma$}}}+\rho_{\mbox{\scriptsize{$\Gamma$}}}{\bf b}_{\mbox{\scriptsize{$\Gamma$}}}\cdot{\bf v}_{\mbox{\scriptsize{$\Gamma$}}}+\left\llbracket-\rho\left(e+\frac{1}{2}|{\bf v}|^{2}\right)({\bf v}{-}{\bf v}_{\mbox{\scriptsize{$\Gamma$}}})+\mathbb{T}{\bf v}-{\bf q}\right\rrbracket\cdot\mathrm{{\bf n}}_{\mbox{\scriptsize{$\Gamma$}}}\ ,$$

$$\displaystyle\frac{D_{\mbox{\scriptsize{$\Gamma$}}}(\rho_{\mbox{\scriptsize{$\Gamma$}}}e_{\mbox{\scriptsize{$\Gamma$}}})}{Dt}+\rho_{\mbox{\scriptsize{$\Gamma$}}}e_{\mbox{\scriptsize{$\Gamma$}}}(\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\tau}-2K_{M}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\mathrm{n}})$$

$$\displaystyle+\left\llbracket-\rho\left(e+\frac{1}{2}|{\bf v}_{\mbox{\scriptsize{$\Gamma$}}}-{\bf v}|^{2}\right)({\bf v}{-}{\bf v}_{\mbox{\scriptsize{$\Gamma$}}})+\mathbb{T}({\bf v}{-}{\bf v}_{\mbox{\scriptsize{$\Gamma$}}})-{\bf q}\right\rrbracket\cdot\mathrm{{\bf n}}_{\mbox{\scriptsize{$\Gamma$}}}\ ,$$

$$\frac{D_{\mbox{\scriptsize{$\Gamma$}}}(\rho_{\mbox{\scriptsize{$\Gamma$}}}\eta_{\mbox{\scriptsize{$\Gamma$}}})}{Dt}+\rho_{\mbox{\scriptsize{$\Gamma$}}}\eta_{\mbox{\scriptsize{$\Gamma$}}}\left(\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\tau}{-}2K_{M}{\bf v}_{\mbox{\scriptsize{$\Gamma$}},\mathrm{n}}\right)+\hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{$\Gamma$}}}\mbox{\boldmath$\Phi$}_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\eta$}}}+\left\llbracket{\mbox{\boldmath$\Phi$}^{\mbox{\tiny{$\eta$}}}}{+}\rho\eta({\bf v}{-}{\bf v}_{\mbox{\scriptsize{$\Gamma$}}})\right\rrbracket\cdot\mathrm{{\bf n}}_{\mbox{\scriptsize{$\Gamma$}}}{-}\Sigma_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\eta$}}}=\Pi_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\eta$}}}\ ,$$

$$\Pi_{\mbox{\scriptsize{$\Gamma$}}}^{\mbox{\tiny{$\eta$}}}\ \geq 0\ .$$

$$e=\hat{e}(\eta ,\rho ,\nabla \rho ).$$

$$\vartheta \stackrel{\mathrm{def}}{=}\frac{\partial \hat{e}}{\partial \eta }.$$

$$\displaystyle\psi=e-\vartheta\eta\ ,\hskip 28.45274pt\ \ \widehat{\psi}(\vartheta,\rho,\nabla\rho)\$$

$$\eta =-\frac{\partial \hat{\psi }}{\partial \vartheta }.$$

$$\frac{\partial \hat{\psi }}{\partial \rho }\frac{D\rho }{Dt}+\frac{\partial \hat{\psi }}{\partial \nabla \rho }\cdot \frac{D\nabla \rho }{Dt}=\frac{De}{Dt}-\vartheta \frac{D\eta }{Dt}.$$

$$\rho \left(\vartheta \frac{D\eta }{Dt}+\frac{p}{{\rho }^2}\frac{D\rho }{Dt}+\frac{\partial \hat{\psi }}{\partial \nabla \rho }\cdot \frac{D\nabla \rho }{Dt}\right)=-\mathrm{div}\mathbf{q}+\mathbb{T}:\mathbb{D}+\rho r,$$

$$p\stackrel{\mathrm{def}}{=}{\rho }^2\frac{\partial \hat{\psi }}{\partial \rho }.$$

Full text

††journal: Int J Eng Sci.

1 Introduction

The seminal papers by Dutch scientists Johannes Diederik van der Waals and Diederik Johannes Korteweg at the turn of the 19th century (van der Waals, 1893; Korteweg, 1901) provided the first thermodynamic insight into the physics of capilarity. In their theory, interaction phenomena at the interfaces between liquid and vapor phases of one substance are described in terms of properties of an interfacial zone of finite thickness where density changes continuously albeit with a very steep gradient. A cornerstone of their theory can be formulated as the assumption that the Helmholtz free energy of such a two-phase system is composed of two contributions - a (local) double well part with two minima related to the two coexisting phases and a gradient term penalizing the volume of the interfacial regions, the latter term being related to the notion of surface energy and surface tension. A considerable effort has been spent in an attempt to incorporate these ideas consistently into the framework of continuum mechanics and thermodynamics and to couple these models of capillarity with equations of flow (e.g. Dunn and Serrin, 1986; Anderson et al., 1998; Heida and Málek, 2010).

Korteweg-type models have gained great popularity in the modeling of granular materials and also in the modeling of two-phase flows; see survey papers Hutter and Rajagopal (1994); Rohde (2018). A key feature of Korteweg-type models is their ability to naturally deal with complex changes of domain topology in contrast with the sharp interface counterparts of these models. On the other hand, their apparent disadvantage is due to the presence of Korteweg stress in the balance of linear momentum that calls for additional boundary conditions which are very difficult to specify in an ad-hoc manner. This issue is clearly not just a mathematical subtlety. In the discussed class of models the boundary conditions describe real physical phenomena such as the motion of the contact line, i.e., dynamics of advancement or retreat of the vapor-fluid interface attached to the solid surface, see for instance Heida (2013) and references therein, in particular Bonn et al. (2009). Another observable real-world phenomenon most likely related to boundary conditions is the so-called contact angle hysteresis, that is, the difference in the measured contact angles of sliding droplets on the advancing and receding parts of the contact line (see e.g. Bormashenko, 2013). Furthermore, one expects that a formulation of the boundary conditions based on solid physical grounds would result in formulations of the problems that might be robust from the point of view of computer simulations and amenable from the point of view of mathematical analysis.

The general aim of this paper is to address the question of the identification of appropriate boundary conditions for problems in continuum thermodynamics. Towards this goal, we use a transparent thermodynamic approach that has been successful in identification of the constitutive equations in the bulk for various complex materials and that stems from specification of the energy storage and dissipation mechanisms. Here, we follow a similar methodology, but we extend it also to surface phenomena. A crucial viewpoint adopted here is that the outer boundary of a liquid-vapor body may be viewed as an interface between this body and its exterior. This viewpoint provides a framework for considering a rather general class of boundary processes and admits a natural coupling between the processes on the surface and in the bulk. This in turn leads to a relatively straightforward procedure for deriving the constitutive relations on the surface delimiting the boundary, i.e., the boundary conditions. This approach is illustrated on the derivation of boundary conditions for a Korteweg-type fluid, for which we can explicitly characterize both the bulk and the surface Helmholtz free energies - in the bulk using the standard thermodynamic relations for van der Waals fluid and at the surface by exploiting the idea of wall-interaction energy for diffuse-interface models (Jacqmin, 2000). Let us, however, note that the methodology developed in this paper can be extended in a relatively straightforward manner to other diffuse interface (or order parameter) models, such as Cahn-Hilliard or Allen-Cahn models; see the concluding remarks in the final section.

The structure of the paper is as follows. In Section 2, we first formulate a general integral form of the balance equation for a quantity comprising bulk and interfacial contributions and provide a corresponding local form of the balances in the bulk and at the interfaces. We explicitly list the local forms of balance equations for mass, linear momenta, energy and entropy for a single-component body. In Section 3, we recapitulate the thermodynamic derivation of a constitutive model for Korteweg-type fluid in the bulk following and slightly modifying the approach of Heida and Málek (2010). In Section 4, we extend this approach to surface phenomena, and by mutual coupling between the bulk and the interface processes, we identify the surface entropy production and the surface entropy flux. The surface entropy production is then rearranged into the form of a sum of the products of mutually related quantities (sometimes called thermodynamic fluxes and thermodynamic affinities) where the individual terms represent different physical mechanisms. Requiring that these mutually related quantities are linearly related111To our understanding, it means that we have employed the framework of linear irreversible thermodynamics (de Groot and Mazur, 1984). (with positive coefficient of proportionality), we not only specify the convex quadratic form for the entropy production, but we also obtain linear constitutive relations on the boundary that automatically comply with the second law of thermodynamics. These constitutive relations (i.e., the boundary conditions) involve a novel type of static and dynamic contact angle boundary conditions as well as a non trivial generalization of the Navier slip boundary condition for the Korteweg model. In Section 5, we discuss a particular variant of the model obtained by considering the bulk Helmholtz free-energy of a van der Waals fluid under isothermal conditions. In Section 6, we present several numerical experiments which demonstrate in a simplified two-dimensional setting the effects of the obtained contact angle and generalized Navier slip boundary conditions for the Korteweg - van der Waals model and show its potential to model dynamic contact angle phenomena and in particular the contact angle hysteresis.

2 General local form of the balance equations in a body with singular surface

Let us consider a material body $\mathcal{B}$ in the current configuration which contains a singular surface $\Gamma$. The singular surface is understood as a mathematical model for a thin wall or membrane which separates one part of the body from another. On the surface $\Gamma$, counterparts of bulk properties and processes may take place. Let us consider arbitrary control volume $V$ and let \mbox{\boldmath\Psi}(V) denote a generic additive quantity (such as mass, momentum, energy, etc.) contained in $V$. Let us consider an integral form of a general balance equation for such a quantity, evaluating the rate of change of \mbox{\boldmath\Psi}(V), as a result of three independent processes: (i) a flux \mathcal{F}^{\mbox{\tiny{\Psi}}} of the quantity $\Psi$ through the boundaries $\partial V$ of the control volume $V$, (ii) an internal production \mathcal{P}^{\mbox{\tiny{\Psi}}} of the quantity $\Psi$ within the control volume $V$ and (iii) an outer supply \mathcal{S}^{\mbox{\tiny{\Psi}}} of the quantity $\Psi$ to the control volume $V$. The general balance equation is thus postulated in the form

[TABLE]

Although this general description is rather formal, we feel it is useful to see all the balance equations of continuum thermodynamics under a unifying frame. The quantity $\Psi$ is assumed to be composed of a bulk contribution and a surface contribution localized at a singular surface $\Gamma$ (see Fig. 1) and both the bulk and the surface contributions are assumed to be representable by corresponding densities $\Psi_{V}$ and $\Psi_{\Gamma}$, respectively. Here, we implicitly follow the standard notion of Gibbs’ surface excess when discussing the surface quantities (Gibbs, 1928).

Considering a control volume $V$ as in Fig. 1, we thereby assume the following representation:

Quantity $\Psi$:

[TABLE]

• 2.

Flux \mathcal{F}^{\mbox{\tiny{\Psi}}}:

[TABLE]

where $\mathrm{{\bf n}}$ is the outer unit normal to the boundary $\partial V$ and $\nu$ is the outer unit normal to the line $\partial\Gamma$ (lying in $\Gamma$), and $S^{+}$, $S^{-}$ are the “outer” parts of $\partial V^{+}$, $\partial V^{-}$ in the sense that $S^{+}{\stackrel{{\scriptstyle\mathrm{def}}}{{=}}}\partial V^{+}{\cap}\partial V$ and $S^{-}{\stackrel{{\scriptstyle\mathrm{def}}}{{=}}}\partial V^{-}{\cap}\partial V$.

• 3.

Production \mathcal{P}^{\mbox{\tiny{\Psi}}}:

[TABLE]

• 4.

Supply \mathcal{S}^{\mbox{\tiny{\Psi}}}:

[TABLE]

We will distinguish between the volume (bulk) contribution $\Psi_{V}$ and its surface counterpart \Psi_{\mbox{\scriptsize{\Gamma}}}, in the sense that in general222The symbol $A|_{S}$ stands for the restriction of the quantity $A$ to a set $S$. \Psi{{}_{V}}|_{{}_{\mbox{\scriptsize{\Gamma}}}}\neq\Psi_{\mbox{\scriptsize{\Gamma}}}, that is, the restriction of the bulk quantity to the surface need not coincide with the corresponding surface quantity. This assumption corresponds to the fact that interfaces are in general zones where material properties may undergo abrupt changes and the transition layers are typically very thin. It is reasonable to treat them as $n{-}1$ dimensional manifolds, $n$ being the dimension of the “bulk” space, and the corresponding averaged (over the thickness of the layer) bulk quantities are then taken as the independent surface counterparts.

By taking arbitrary control volume $V$, using the generalized Reynolds’ transport theorem, Gauss’ theorem, and tools of the differential geometry, under an additional assumption of sufficient smoothness of all the involved quantities we can derive the following local form of the integral balance equations (Slattery, 1990, section 1.3.2):

• [TABLE]

In the bulk $\mathcal{B}^{+}{\cup}\mathcal{B}^{-}$:

[TABLE]

where ${\bf v}$ is the material velocity.

• 2.

At the interface $\Gamma$:

[TABLE]

where \frac{D_{\mbox{\scriptsize{\Gamma}}}}{Dt} denotes the surface material time derivative defined as

[TABLE]

where {\bf X}_{\mbox{\scriptsize{\Gamma}}} denotes the surface material point, (see Slattery, 1990, section 1.2.5). Next, {\bf v}_{\mbox{\scriptsize{\Gamma}}} is the surface velocity and {\bf v}_{\mbox{\scriptsize{\Gamma}},\tau} is its projection to the surface $\Gamma$ and {\bf v}_{\mbox{\scriptsize{\Gamma}},\mathrm{n}} is the normal component:

[TABLE]

where $\mathbb{I}$ is the identity tensor, $K_{M}$ is the mean curvature of the surface, and $\llbracket z\rrbracket$ denotes the jump of the bulk quantity $z$ across the surface $\Gamma$ defined by $\llbracket z\rrbracket:=z^{+}-z^{-}$ where $z^{+}$ and $z^{-}$ are the restrictions of $z|_{\mathcal{B}^{+}}$ and $z|_{\mathcal{B}^{-}}$, respectively, to $\Gamma$, see Fig. 1. Note that in (1b), we keep the surface velocity {\bf v}_{\mbox{\scriptsize{\Gamma}}} inside the “jump” brackets $\llbracket\rrbracket$. This is standard notation in the literature understood in the sense that all surface quantities are tacitly taken as continuous, i.e. {\bf v}_{\mbox{\scriptsize{\Gamma}}}^{+}={\bf v}_{\mbox{\scriptsize{\Gamma}}}^{-} and thus the term on the right hand side in (1b) is interpreted as \llbracket\Psi_{V}({\bf v}-{\bf v}_{\mbox{\scriptsize{\Gamma}}})\rrbracket\cdot\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}}=\llbracket\Psi_{V}{\bf v}\rrbracket\cdot\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}}-\llbracket\Psi_{V}\rrbracket{\bf v}_{\mbox{\scriptsize{\Gamma}}}\cdot\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}}. Finally, \hskip 0.85358pt\mathrm{div}_{\mbox{\scriptsize{\Gamma}}} denotes the surface divergence operator. For the definition of kinematic quantities and operators on the surfaces, see Slattery (1990)(Appendix A) or the coordinate-free exposition by Buscaglia and Ausas (2011).

2.1 Local forms of the balance equations in the bulk

The local forms of the balance equations in the bulk $\mathcal{B}^{+}\cup\mathcal{B}^{-}$ (i.e. outside the iterface $\Gamma$) for a single-component non-polar material read as follows (we omit the subscript V for bulk quantities for brevity):

• [TABLE]

Balance of mass:

[TABLE]

where $\rho$ is the density and ${\bf v}$ is the material velocity.

• 2.

Balance of linear momentum:

[TABLE]

where $\mathbb{T}$ is the Cauchy stress tensor and ${\bf b}$ is the specific body force. Here and in what follows the dyadic product $\mathbf{a}\otimes\mathbf{d}$ of two vectors $\mathbf{a}$ and $\mathbf{d}$ is the second order tensor with the components $(\mathbf{a}\otimes\mathbf{d})_{ij}=a_{i}d_{j}$.

• 3.

Balance of angular momentum is reduced to the statement that the Cauchy stress tensor $\mathbb{T}$ is symmetric, i.e.,

[TABLE]

where T denotes the transposition of a tensor.

• 4.

Balance of energy:

[TABLE]

where $e$ is the specific internal energy, ${\bf q}$ denotes the energy flux, and $r$ is the specific energy supply. Employing the mass and momentum balances, we obtain the balance equation for the energy in the form:

[TABLE]

where $\mathbb{D}=\frac{1}{2}(\nabla{\bf v}+(\nabla{\bf v})^{\mathrm{T}})$ is the symmetric part of the velocity gradient and $\frac{D}{Dt}=\frac{\partial}{\partial t}+{\bf v}\cdot\nabla$ denotes the material time derivative.

• 5.

The formulation of the second law of thermodynamics:

[TABLE]

where $\eta$ is the specific entropy, \mbox{\boldmath\Phi}^{\mbox{\tiny{\eta}}} is the bulk entropy flux, \Sigma^{\mbox{\tiny{\eta}}} is the entropy supply and \Pi^{\mbox{\tiny{\eta}}} is the entropy production, which must be non-negative according to the second law of thermodynamics.

2.2 Local form of the balance quations at the interface

The local forms of the balance equations at the interface $\Gamma$ for a single-component non-polar material read as follows:

• [TABLE]

Balance of mass:

[TABLE]

where \rho_{\mbox{\scriptsize{\Gamma}}} is the surface mass density.

• 2.

Balance of linear momentum:

[TABLE]

where \mathbb{T}_{\mbox{\scriptsize{\Gamma}}} denotes the surface Cauchy stress tensor and {\bf b}_{\mbox{\scriptsize{\Gamma}}} is the specific surface force.

• 3.

Balance of angular momentum (for a non-polar material):

[TABLE]

i.e., symmetry of the surface Cauchy stress tensor.

• 4.

Balance of energy:

[TABLE]

where e_{\mbox{\scriptsize{\Gamma}}} denotes the specific surface internal energy and {\bf q}_{\mbox{\scriptsize{\Gamma}}} is the surface energy flux. Multiplying the surface momentum balance (5b) by {\bf v}_{\mbox{\scriptsize{\Gamma}}}, one can obtain the surface balance equation for kinetic energy. This can be subtracted from (5d), which yields the balance of surface energy in the reduced form

[TABLE]

where \mathbb{T}_{\Gamma}{:}\nabla_{\Gamma}{\bf v}_{\Gamma}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathrm{tr}\left(\mathbb{T}_{\mbox{\scriptsize{\Gamma}}}\ \nabla_{\mbox{\scriptsize{\Gamma}}}{\bf v}_{\mbox{\scriptsize{\Gamma}}}\right).

• 5.

The formulation of the second law of thermodynamics:

[TABLE]

with

[TABLE]

Here \eta_{\mbox{\scriptsize{\Gamma}}} denotes the surface specific entropy, \mbox{\boldmath\Phi}_{\mbox{\scriptsize{\Gamma}}}^{\mbox{\tiny{\eta}}} is the surface entropy flux, \Sigma_{\mbox{\scriptsize{\Gamma}}}^{\mbox{\tiny{\eta}}} is the surface entropy supply and \Pi_{\mbox{\scriptsize{\Gamma}}}^{\mbox{\tiny{\eta}}} is the surface entropy production, which must be non-negative in order to comply with the second law of thermodynamics.

3 Derivation of constitutive equations for Korteweg-type fluids in the bulk

In this section, we recall and slightly modify the derivation of a constitutive model for a Korteweg-type fluid developed in Heida and Málek (2010). The derivation is based on imposing the following constitutive ansatz for the (specific) internal energy in the bulk:

[TABLE]

Assuming that $\hat{e}$ is differentiable, the thermodynamic temperature $\vartheta$ is introduced through

[TABLE]

The corresponding Helmholtz free energy $\psi$ is then obtained via the Legendre transform giving

[TABLE]

where, in the last equality, we assume the invertibility of (7) with respect to $\eta$. As a consequence of (8), we obtain the standard thermodynamic relation

[TABLE]

Taking the material time derivative of (8), we obtain, after using (9), that

[TABLE]

Multiplying (10) by $\rho$ and applying the energy balance (4e) we obtain

[TABLE]

where $p$ denotes the thermodynamic pressure defined through

[TABLE]

Taking the gradient of (4a) we obtain

[TABLE]

Using (4a) and (13) in (11), we arrive at

[TABLE]

where $m$ denotes the mean normal stress $m\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{3}\mathrm{tr}(\mathbb{T})$ and $()^{d}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}()-\frac{1}{3}\mathrm{tr}(){\mathbb{I}}$ denotes the deviatoric part of a tensor. By the principle of material frame indifference (see, e.g., Truesdell and Noll, 1965), the internal energy $\hat{e}$ can only depend on the magnitude of $\nabla\rho$ which immediately implies symmetry of the tensor $\rho\frac{\partial\widehat{\psi}}{\partial\nabla\rho}\otimes\nabla\rho$. Consequently, we may replace $\nabla{\bf v}$ by its symmetric part $\mathbb{D}$ in the fourth term in the right-hand side of (14). Using the identity

[TABLE]

dividing (14) by $\vartheta$ we obtain, after suitable rearrangements, the following local form of the balance equation for the bulk entropy:

[TABLE]

Recalling (4f) and assuming that the entropy supply is given only by the corresponding energy supply term, that is, postulating that333It is also possible to split $\rho r{=}\rho r_{A}{+}\rho r_{B}$ and postulate \Sigma^{\mbox{\tiny{\eta}}}{=}\frac{\rho r_{A}}{\vartheta}, and incorporate $\frac{\rho r_{B}}{\vartheta}$ among the entropy producing mechanisms. For simplicity, we do not consider this possibility here.

[TABLE]

we can identify the entropy flux in (16) as

[TABLE]

The second and last terms in the right-hand side of (16) represent the entropy production. In accordance with the usual approach in the constitutive theory within linear irreversible thermodynamics (de Groot and Mazur, 1984), we want to express this term as a sum of binary products between the thermodynamic “affinities” (forces) and the corresponding thermodynamic “fluxes”. Even if we pick as the set of affinities $(\hskip 0.85358pt\mathrm{div}{\bf v},\mathbb{D}^{d},\nabla(\frac{1}{\vartheta}))$, the splitting into two groups is still not unique. Note, in particular, that the term $\rho^{2}\frac{\partial\widehat{\psi}}{\partial\nabla\rho}\cdot\nabla\left(\frac{1}{\vartheta}\right)\hskip 0.85358pt\mathrm{div}{\bf v}$ can contribute to both products with affinities $\nabla(\frac{1}{\vartheta})$ and $\hskip 0.85358pt\mathrm{div}{\bf v}$. Without knowing a-priori which choice is preferable, we split this term via a convex combination governed by a free parameter $\alpha\in\langle 0,1\rangle$ between both these terms. With such a choice, we finally arrive at the formula for the rate of entropy production,

[TABLE]

where we set

[TABLE]

We recognize the right-hand side of (19) as three entropy-producing mechanisms, each in the form of a product of two terms (sometimes called thermodynamic “flux” and thermodynamic “affinity”), each couple describing a different physical process. Restricting ourselves here to linear relationships among the two types of terms, we arrive at the constitutive equations

[TABLE]

If these relationships are inserted back into (19), we obtain the rate of entropy production expressed as a piece-wise quadratic function in terms of the “affinities”,

[TABLE]

(here and in what follows $|\mathbb{A}|\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sqrt{\mathbb{A}_{ij}\mathbb{A}_{ij}}$), or equivalently, in terms of “fluxes” as follows

[TABLE]

In both cases we see that the positivity of the coefficients together with the piece-wise quadratic form of (22) and (23) ensure that the second law of thermodynamics, i.e., the non-negativity of the rate of entropy production, is guaranteed. Inserting the formulas (20) into (21), we obtain the following expressions for the Cauchy stress $\mathbb{T}$, energy flux ${\bf q}$, and entropy flux \mbox{\boldmath\Phi}^{\mbox{\tiny{\eta}}} in the bulk:

[TABLE]

where $\alpha\in\langle 0,1\rangle$, $2\mu+3\lambda>0$, $\mu>0$, $\kappa>0$.

4 Derivation of boundary conditions for Korteweg-type fluids

In this section, we extend the constitutive theory for Korteweg-type fluids presented in the previous section to dissipative processes at the boundary which will allow us to formulate thermodynamically based constitutive relations (compatible with the (local form) of the second law of thermodynamics) in the form of boundary conditions. Let us now consider a domain $\Omega$ which contains Korteweg (two-phase) fluid and let us denote its boundary $\partial\Omega$. We will look at the boundary $\partial\Omega$ as an interface between the domain $\Omega$ and its exterior; see Fig. 2. Adopting this viewpoint, we can employ the framework of continuum mechanics with singular surfaces introduced in Section 2.

For Korteweg-type fluids we assume, in accordance with the physical theories of capillarity (e.g. Rowlinson and Widom, 1989), that the outer boundary is in fact a boundary layer with certain specific properties. This layer will be treated as infinitely thin and all the corresponding bulk quantities in this layer will be described by their surface (boundary) counterparts, obtained by space averaging over the thickness of the boundary layer. Unlike Navier-Stokes fluids, Korteweg fluids naturally incorporate the notion of surface tension as the interfacial energy in transition regions separating the phases. With the goal of formulating boundary conditions for Korteweg-type fluids and describing phenomena such as wetting (i.e., contact angles), it appears reasonable to extend the notion of interfacial interaction and to include also interaction of the fluid with the boundary walls. To this end, we postulate the existence of boundary surface energy e_{\mbox{\scriptsize{\Gamma}}}, and boundary surface entropy \eta_{\mbox{\scriptsize{\Gamma}}} expressing this fluid-boundary interaction. By postulating the different constitutive equations for e_{\mbox{\scriptsize{\Gamma}}}, we will obtain a hierarchy of models of various complexity, which will result in a corresponding hierarchy of classes of boundary conditions.

For simplicity, we will investigate a model in which we ignore convective mechanisms on the boundary and we shall thus consider the boundary $\partial\Omega$ to be static by postulating zero surface velocity, i.e.,

[TABLE]

Next, we will also need to specify conditions on the exterior side of the boundary (denoted by a $+$ sign with the convention of exterior unit normal pointing from $-$ to $+$). We will assume that the material outside $\Omega$ is at rest, i.e.,

[TABLE]

implying that any mass exchange between the exterior and interior is excluded. Concerning momentum exchange, the exterior may exert force on the interior domain, but due to (25b), this force does not produce any mechanical power and as such does not contribute to the energy balance. We will, however, assume that the exterior may facilitate (non-convective) energy and entropy transfer and we relate them through the standard relation of thermodynamics (Coleman and Noll, 1963):

[TABLE]

Furthermore, we assume that from the “inside” the boundary is impermeable, i.e.,

[TABLE]

allowing for slip of the Korteweg fluid along the boundary but no penetration.

In the following, it will be convenient to not explicitly invoke the concept of surface mass density, as we will only require the notions of surface energy and surface entropy. Thus instead of using the specific (i.e. related to unit of mass) surface energy and surface entropy, we formulate the equations directly for the products \rho_{\mbox{\scriptsize{\Gamma}}}e_{\mbox{\scriptsize{\Gamma}}} and \rho_{\mbox{\scriptsize{\Gamma}}}\eta_{\mbox{\scriptsize{\Gamma}}}, using the following notation:

[TABLE]

Similarly as in the bulk, the key constitutive relation representing the assumption of a local thermodynamic equlibrium takes the form

[TABLE]

where the dots stand for other state variables. We then define the surface thermodynamic temperature \vartheta_{\mbox{\scriptsize{\Gamma}}} through

[TABLE]

We also introduce the corresponding surface Helmholtz free energy \widetilde{\psi}_{\mbox{\scriptsize{\Gamma}}} via the Legendre transform obtaining

[TABLE]

where in the last equality we assume invertibility of (28) with respect to \widetilde{\eta}_{\mbox{\scriptsize{\Gamma}}}. It then follows from (29) that

[TABLE]

Concerning the structure of the constitutive equation for the surface Helmholtz free energy, we will consider the following three situations:

• [TABLE]

Model A:

Motivated by statistical physics description of the surface free energy (see Rowlinson and Widom (1989)(eq. 4.114)) we assume that

[TABLE]

i.e., the surface Helmholtz free energy depends not only on the surface temperature, but also on the density of the neighboring bulk. Here only dependence on the interior density $\rho^{-}$ is considered.

• 2.

Model B:

[TABLE]

i.e., dependence of the surface Helmholtz free energy only on surface temperature.

• 3.

Model C:

[TABLE]

i.e., the “trivial” model without any surface Helmholtz free energy, which however still provides non-trivial boundary conditions for the bulk terms.

We will now inspect these three cases in detail.

4.1 Model A

Applying the surface material time derivative to (29) with \widetilde{\psi}_{\mbox{\scriptsize{\Gamma}}} given by (31a), using (30), and employing the surface energy balance in the form (5e), we obtain the identity

[TABLE]

The surface being static (see the assumption (25a), all terms containing {\bf v}_{\mbox{\scriptsize{\Gamma}}} vanish and the surface material time derivative in (32) becomes just the partial time derivative, i.e. \frac{D_{\mbox{\scriptsize{\Gamma}}}}{Dt}=\frac{\partial}{\partial t}. We then use (4a) and (25d) and obtain

[TABLE]

With the assumptions (25a–25d) and neglecting for simplicity the surface “body” forces, i.e.,

[TABLE]

the balance of linear momentum on the surface (5b) reads

[TABLE]

We shall consider a membrane model444A generalization that would involve more complex structure of the surface stress tensor or non-constant surface tension is possible, but not straightforward. We will therefore not pursue this possiblity here; see also a comment in the conclusion., where only the surface tension $\sigma$ (which is assumed to be constant here for simplicity) constitutes the surface stress tensor, i.e., we consider

[TABLE]

Then it holds (Slattery, 1990, Appendix A) that

[TABLE]

where $K_{M}$ is the mean curvature of the surface, and thus by (35), we get

[TABLE]

where $()_{\tau}$ denotes the projection of a vector to the tangent plane. Hence

[TABLE]

Due to (25b) and (25d) and by virtue of continuity of tangent traction (39) and the symmetry of $\mathbb{T}$, it holds that

[TABLE]

Assuming Korteweg fluid inside $\Omega$, we have, by (20a) and (21c), the following expression for the energy flux

[TABLE]

Inserting (40) and (41) into (33), we finally obtain

[TABLE]

where we set

[TABLE]

In (43a), we used the fact that \nabla\rho^{-}\cdot{\bf v}^{-}=\nabla_{\mbox{\scriptsize{\Gamma}}}\rho^{-}\cdot{\bf v}^{-}_{\tau}, due to (25d).

The next step consists of transforming (42) into the form (5f) where we must also incorporate all the simplifying assumptions used above, such as (25a)-(25d). Then (5f) takes the form

[TABLE]

which, upon inserting (24c) and (25c), leads to

[TABLE]

We proceed in two different ways. These ways differ in the manner in which the term $\hskip 0.85358pt\mathrm{div}{\bf v}^{-}$ in (42) and (45) is treated. In the first procedure, called Model A1, $\hskip 0.85358pt\mathrm{div}{\bf v}^{-}$ is kept unaltered, while in the second procedure, called Model A2, we will split $\hskip 0.85358pt\mathrm{div}{\bf v}^{-}$ into the surface divergence and the normal derivative.

4.1.1 Model A1

We first observe that the equation (42) can be rewritten in the following form of entropy balance:

[TABLE]

Subtracting (46) from (45) yields

[TABLE]

Since \widetilde{r}_{\mbox{\scriptsize{\Gamma}}} is the surface energy supply, it is reasonable to postulate the surface entropy supply to be

[TABLE]

Since {\bf q}_{\mbox{\scriptsize{\Gamma}}} is the surface energy flux, classical thermodynamics together with (47) suggest setting

[TABLE]

Consequently, (47) reduces to the equation

[TABLE]

which identifies the entropy-producing mechanisms and which has the usual structure of a scalar product

[TABLE]

where we choose

[TABLE]

We propose linear constitutive relations between the “fluxes” $\mathbf{J}$ and “affinities” $\mathbf{A}$, following thus the framework of linear irreversible thermodynamics (de Groot and Mazur, 1984). We shall also consider possible cross-coupling among the vectorial quantities. The constitutive relations in such case take the form

[TABLE]

Since the two affinities -\frac{{\bf v}_{\tau}^{-}}{\vartheta_{\mbox{\scriptsize{\Gamma}}}}, and \nabla_{\mbox{\scriptsize{\Gamma}}}\left(\frac{1}{\vartheta_{\mbox{\scriptsize{\Gamma}}}}\right), for which cross-effect is assumed, have opposite behavior with respect to time reversal (the first one changes the sign, the other does not), the Onsager-Casimir relations (see, e.g., (de Groot and Mazur, 1984)) imply anti-symmetry of the cross-coupling coefficient, i.e., $L_{12}=-L_{21}$. The coefficients must fulfill $L_{ii}\geq 0$, for all $i{=}1,{\dots},5$ and $L_{11}L_{22}+(L_{12})^{2}\geq 0$, in order to ensure non-negativity of the rate of entropy production. We shall impose the stronger yet natural assumption that $L_{ii}{>}0$ for all $i{=}1,{\dots},5$ in order to avoid degeneracy of the system.

Let us interpret the derived constitutive relations (53). The first relation (53a) represents the in-surface heat conduction (Fourier law). The last two relations (53d) and (53e) represent heat transfer across the interface, the so-called Kapitza resistance (Kapitza, 1941). Condition (53b) represents a generalized Navier-slip condition. Inspecting (43a), we can see that it is a relation among the surface traction (\mathbb{T}\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}})^{-}, a cross-coupling term involving the surface gradient of temperature, a term proportional to slip velocity ${\bf v}_{\tau}$, i.e., traditional Navier-slip term, and finally, a term involving the surface Helmholtz free energy and tangent derivative of density. Perhaps the most interesting is the relation (53c), which we will later interpret as the static and dynamic contact angle condition - a condition characterizing the angle between the liquid-vapor interface and the boundary. This interpretation will be made explicit in Sections 5 and 6, where we will consider a particular type of bulk and surface Helmholtz free energy functions and support our arguments with numerical experiments. This condition in that case will relate the normal derivative of the density \frac{\partial\rho^{-}}{\partial\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}}} with two effects - one due to surface tension and the other due to motion of the fluid in the vicinity of the interface. Since the latter effect vanishes when the body is in equilibrium, we will interpret that part as the dynamic contact angle condition, while the first effect persists in equilibrium and will be called the static contact angle condition. Let us note here, that in most of the literature related to Korteweg-type models, the dynamic contact angle condition is completely ignored, and the static one is simplified dramatically to \frac{\partial\rho^{-}}{\partial\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}}}{=}0, which corresponds to the contact angle $\frac{\pi}{2}$.

4.1.2 Model A2

The second approach is based on the decomposition of the term $\hskip 0.85358pt\mathrm{div}{\bf v}^{-}$ appearing in (42) and (45) by means of the following identity from differential geometry (Slattery, 1990, Appendix A):

[TABLE]

where

[TABLE]

where we used (25d) in the last equality in (54). We will also employ the identity

[TABLE]

Incorporating (54) and (56) into (42), we obtain

[TABLE]

which leads to the following balance equation for the entropy:

[TABLE]

Similarly, applying (54) together with the identity

[TABLE]

to (45), we arrive, after rearranging the terms, at the following form of the entropy balance:

[TABLE]

Subtracting (58) from (60), we obtain

[TABLE]

As in Model A1 we postulate the surface entropy supply to be

[TABLE]

and (61) suggests setting

[TABLE]

Consequently, (61) gives

[TABLE]

or, written again as a scalar product of two vectors,

[TABLE]

where

[TABLE]

As in Subsection 4.1.1, restricting ourselves to the linear constitutive relations between the “fluxes” ${\bf J}$ and “affinities” ${\bf A}$ (with cross-coupling only among the vectorial quantities), we end up with the following set of constitutive relations

[TABLE]

Since the two “affinities” -\frac{{\bf v}_{\tau}^{-}}{\vartheta_{\mbox{\scriptsize{\Gamma}}}}, and \nabla_{\mbox{\scriptsize{\Gamma}}}\left(\frac{1}{\vartheta_{\mbox{\scriptsize{\Gamma}}}}\right), for which cross-effect is assumed, have opposite behavior with respect to time reversal, the Onsager-Casimir relations suggest the requirement that the cross-coupling coefficients are anti-symmetric, i.e., $L_{12}{=}-L_{21}$. The coefficients are assumed to satisfy $L_{ii}{>}0$, $i{=}1,{\dots},5$, and $L_{11}L_{22}+(L_{12})^{2}\geq 0$, in order to ensure both non-negativity of the rate of entropy production and non-degeneracy of the system.

The interpretation of the constitutive relations for Model A2 is analogous to Model A1, namely the first two relations (67a) and (67b) represent the (generalized) in-surface heat conduction (Fourier law) and generalized Navier-slip condition, respectively, together with a possible cross-coupling of the two mechanisms. Relations (67d) and (67e) represent heat transfer across the interface (the so-called Kapitza resistance) and relation (67c) is again the static and dynamic contact angle condition, as will become apparent in Sections 5 and 6.

4.2 Model B

The derivation of boundary conditions for Model B proceeds in an analogous way as for Model A. The only difference between these models is the absence of the term $s_{\rho}$, which is now identically zero. Consequently, we obtain the following two sets of boundary conditions, which again differ in the manner how the terms involving $\hskip 0.85358pt\mathrm{div}{\bf v}^{-}$ are treated.

4.2.1 Model B1

[TABLE]

where $L_{12}{=}-L_{21}$ and where $L_{ii}{>}0$, $i{=}1,{\dots},5$ and $L_{11}L_{22}{+}(L_{12})^{2}\geq 0$, in order to ensure both non-negativity of the rate of entropy production and non-degeneracy of the system.

4.2.2 Model B2

[TABLE]

where $L_{12}{=}-L_{21}$ and where $L_{ii}{>}0$, $i{=}1,{\dots},5$ and $L_{11}L_{22}{+}(L_{12})^{2}\geq 0$, in order to ensure both non-negativity of the rate of entropy production and non-degeneracy of the system.

The only difference between Models A and B is due to the absence of terms $s_{\rho}$, its main implication being that conditions (68c) and (69c) represent solely dynamic angle conditions, with static (equilibrium) contact angle (equal to $\frac{\pi}{2}$ for the Korteweg - van der Waals fluid studied in Section 5).

4.3 Model C

Assuming that the surface Helmholtz free energy, and consequently also both the surface internal energy and surface entropy are identically equal to zero, it also makes sense to assume the same for the corresponding energy and entropy surface fluxes. Therefore we set

[TABLE]

Employing also the assumptions on the velocity field (25a)–(25d), the absence of the surface body forces (34), and the character of the surface stress tensor (36), the surface energy balance (5e) reduces to the following form of a jump condition:

[TABLE]

Following the same arguments, the surface entropy balance (5f) reduces, with the use of (24b), (24c), (70), and (25a)–(25d), to

[TABLE]

Employing the identity for the jump of a product of two fields $a,b$

[TABLE]

where $\langle a\rangle$ denotes the average of $a^{+}$ and $a^{-}$, i.e., $\langle a\rangle=\frac{1}{2}\left(a^{+}+a^{-}\right)$, applying this identity to \left\llbracket\frac{{\bf q}}{\vartheta}\right\rrbracket in (72), and inserting its result into (71) instead of $\llbracket{\bf q}\rrbracket$, we obtain the equation for the rate of surface entropy production

[TABLE]

This expression again takes the form of a scalar product

[TABLE]

where we choose

[TABLE]

The linear relations between the “fluxes” ${\bf J}$ and “affinities” A yield the following relations:

[TABLE]

Assuming that $L_{11}{>}0$, $L_{22}{>}0$, $L_{33}{>}0$, the rate of entropy production is non-negative. The boundary condition (77a) represents the Navier-slip, (77b) describes the dynamic contact angle condition, and (77c) stands for the heat transfer across the interface (Kapitza resistance), respectively, see Sections 5 and 6 for further details.

**Isothermal process

**In the following, we will study a variant of the above models in which the temperature is continuous across the interface. This can in particular be achieved if we consider an isothermal process and assume that

[TABLE]

If we also ignore the cross-coupling effects for simplicity and absorb the (constant) temperature into the coefficients, we obtain the following sets of reduced boundary conditions (recall that ${\bf t}_{\rho}$, $s_{\rho}$ and $s_{\nabla\rho}$ are defined in (43)):

Model A1:

[TABLE]

• 2.

Model A2:

[TABLE]

• 3.

Model B1 $\&$ Model C:

[TABLE]

• 4.

Model B2

[TABLE]

where $L_{22}{>}0$ $L_{33}{>}0$ for all sets. Each set consists of a (generalized) Navier-slip condition and a contact angle condition, which is either solely dynamic or combines static and dynamic terms, as clarified in the following sections.

5 Particular example of Korteweg fluid model and boundary conditions

Let us now consider a particular Helmholtz free energy $\widehat{\psi}$ of the form

[TABLE]

where the term $\widehat{\psi}_{vdW}$ corresponds to van der Waals fluid (van der Waals, 1893; Landau and Lifshitz, 1980; Diehl, 2007) and takes the following form:

[TABLE]

where $a^{\prime}$, $b^{\prime}$, $c$, $d$, $e$, $\ell$ and $\sigma$ are constant parameters and $\vartheta_{0}$ is some reference temperature. The correspondence of (84) to the van der Waals model is revealed by identifying the equation of state for thermodynamic pressure associated with $\widehat{\psi_{vdW}}$, which we show next.

Using the standard thermodynamic definition of the thermodynamic pressure (12), we obtain the expression

[TABLE]

Defining $a\stackrel{{\scriptstyle\mathrm{def}}}{{=}}a^{\prime}M_{m}^{2}$, $b\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{M_{m}}{b^{\prime}}$ and $R\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\ell M_{m}$, where $M_{m}$ is the molar mass of the molecules of the considered gas-liquid system, and considering a homogeneous system with $n$ moles in volume $V$, which means $\rho=\frac{nM_{m}}{V}$, we can recast (85) into the standard form

[TABLE]

which is the traditional van der Waals equation of state (e.g. Callen, 1985), provided we suitably interpret the parameters $a$, $b$, and $R$.

The chemical potential for a single-component fluid is simply the Gibbs free energy $g$ as follows from the Euler relation (e.g. Callen, 1985). We thus obtain

[TABLE]

The pressure $p_{vdW}$ and the chemical potential $\mu_{vdW}$ are sketched in Fig. 3. For temperatures above a critical temperature $\vartheta_{c}$, both $p_{vdW}$ and $\mu_{vdW}$ are increasing functions of density, but for a temperature below the critical temperature $\vartheta_{c}$, both functions have two increasing branches separated by a region where these functions are decreasing. The critical temperature as well as the critical pressure are found by identifying the inflection point of the $\widehat{p_{vdW}}$ equation of state, i.e., by finding $(\vartheta_{c},\rho_{c})$ such that

[TABLE]

which for (84) yields

[TABLE]

In the subcritical region, i.e., for each $\vartheta{<}\vartheta_{c}$, there are two states (called Maxwell states) $\rho^{M}_{v}$ (vapor) and $\rho^{M}_{l}$ (liquid) defined by two phase-coexistence equilibrium relations expressing the equality of pressures and chemical potentials (see Fig. 3),

[TABLE]

Let us now consider an isothermal setting below the critical temperature (meaning the temperature $\vartheta$ is uniform and equals the constant $\vartheta_{0}$, satisfying $\vartheta_{0}{<}\vartheta_{c}$) and let us set up the system of governing equations and boundary conditions for such a Korteweg - van der Waals fluid in the bulk. Using the definition of the Helmholtz free energy (83), the Cauchy stress $\mathbb{T}$ reads according to (24a) as

[TABLE]

Using the identity

[TABLE]

the governing equations (balances of mass and momentum) in the bulk for the Korteweg - van der Waals fluids read

[TABLE]

Boundary conditions corresponding to Models A, B, and C, (see (79)–(82)) can now be expressed in more explicit forms since for the Korteweg - van der Waals model we can explicitly evaluate the term

[TABLE]

Employing the definition of ${\bf t}_{\rho}$ (43a), we obtain the following set of boundary conditions (depending on the form \widehat{\psi}_{\mbox{\scriptsize{\Gamma}}} and the way how $\hskip 0.85358pt\mathrm{div}{\bf v}^{-}$ is treated on $\Gamma$):

Model A1:

[TABLE]

• 2.

Model A2:

[TABLE]

• 3.

Model B1 $\&$ Model C:

[TABLE]

• 4.

Model B2:

[TABLE]

where $\alpha$, $\beta$ are some non-negative parameters ($\alpha{=}L_{22}$, $\beta{=}\frac{L_{33}}{\rho^{-}}$), possibly depending on $\vartheta^{-}$, \vartheta_{\mbox{\scriptsize{\Gamma}}}, $\rho^{-}$.

It will be convenient to incorporate the constitutive relation for the Cauchy stress (91) in the above conditions. Since, by (91),

[TABLE]

we obtain by simple manipulation (in particular substituting from the second equation into the first) the following conditions:

Model A1:

[TABLE]

• 2.

Model A2:

[TABLE]

• 3.

Model B1 $\&$ Model C:

[TABLE]

• 4.

Model B2:

[TABLE]

where $\alpha$, $\beta$ are non-negative parameters, possibly depending on $\vartheta^{-}$, \vartheta_{\mbox{\scriptsize{\Gamma}}}, $\rho^{-}$.

In the above sets of boundary conditions, the right-hand side of the second equation represents the static contact angle condition (static in the sense that it does not depend explicitly on ${\bf v}$). Let us note that in the class of diffuse interface methods, to which the Korteweg model presented here belongs (as well as other models including Cahn-Hilliard and Allen-Cahn models or numerous variants of the level set method), a standard way to impose a given static contact angle $\varphi$ is expressed through the formula (imposed on the boundary)

[TABLE]

see, e.g., Brackbill et al. (1992). Such a formula, despite its apparent simplicity, is problematic from several points of view. First, it cannot be incorporated into the framework developed above as this would require that the surface Helmholtz free energy \widehat{\psi}_{\mbox{\scriptsize{\Gamma}}} depends also on $|\nabla\rho^{-}|$ instead of just $\rho^{-}$ (and temperature). Second, the term \frac{\partial\rho^{-}}{\partial\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}}} appears naturally in the weak formulation of associated initial boundary value problems (see Section 6.2, eq. (117c)) and from (104)2, we can see that it is a non-linear function of the density gradient. This, however, represents quite a severe constraint on the regularity of the density field in terms of mathematical well-posedness. Last, but not least, in order to apply the formula (104) in numerical calculations, a sufficiently accurate numerical approximation of the term $|\nabla\rho^{-}|$ on the boundary is required. Interestingly, all these problems can be circumvented in the framework developed here. In particular, it is possible to replace (104) by a relation (not involving $\nabla\rho$ at all) of the form

[TABLE]

where $P(\rho^{-})$ is a low-order polynomial and $\gamma(\varphi)$ is a function depending only on the imposed contact angle $\varphi$. Towards this goal, let us follow the so-called energy-based approach (see e.g. Jacqmin, 2000) and postulate the surface Helmholtz free energy as follows:

[TABLE]

where $\sigma_{{vw}}$ and $\sigma_{{lw}}$ are the vapor-wall and liquid-wall surface tensions, respectively, and \psi_{\mbox{\scriptsize{\Gamma}}}^{0} is a constant. In this form, the surface Helmholtz free energy \widehat{\psi}_{\mbox{\scriptsize{\Gamma}}} is constant in the boundary regions that are in contact with the pure bulk phases characterized by the Maxwell states $\rho^{M}_{v}$ or $\rho^{M}_{l}$. The value of the constant differs for the two cases by $(\sigma_{{lw}}{-}\sigma_{{vw}})$, and this jump takes place across the boundary “contact line”, i.e., across the region on the boundary where the phases change from one to another. Let us substitute into (106) the standard contact-angle formula (Young equation)

[TABLE]

relating the liquid-vapor surface tension $\sigma_{{lv}}$ with $\sigma_{{lw}}$ and $\sigma_{{vw}}$ through the cosine of the wetting angle $\varphi$. Here the wetting angle $\varphi$ denotes the contact angle of the liquid-vapor interface with respect to the wall measured inside the liquid domain. Employing the ansatz for the surface Helmholtz free energy (106), the static part (i.e., corresponding to ${\bf v}{=}{\bf 0}$) of the boundary conditions (100b) and (101b) becomes

[TABLE]

which is of the desired form (105). We will test this formula in the numerical simulations in Section 6, where we also provide the specific value for the parameter $\gamma_{0}$. Let us only note here that due to the temperature dependence of the Maxwell states (see (90)), $\gamma_{0}$ depends on temperature even in the current setting with constant surface tensions.

Let us summarize the conditions (100) – (103) corresponding to the ansatz for the surface Helmholtz free energy \widehat{\psi}_{\mbox{\scriptsize{\Gamma}}} of the form (106):

Model A1:

[TABLE]

• 2.

Model A2:

[TABLE]

• 3.

Models B1 and C:

[TABLE]

• 4.

Model B2:

[TABLE]

As will be documented in the following numerical simulations, by explicitly evaluating the parameter $\gamma_{0}$ from (108) in Models A1 and A2, the value $\varphi$ equals the equilibrium (static) contact angle for the Korteweg - van der Waals fluid while the (positive) value of parameter $\beta$ governs the dynamic relaxation to this equilibrium. Clearly, Models B and C admit only homogeneous Neumann boundary conditions in equilibrium \frac{\partial\rho^{-}}{\partial\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}}}{=}0, i.e., the static contact angle is $\frac{\pi}{2}$.

6 Numerical experiments focused on contact angle phenomena

Numerical experiments presented below are focused on the qualitative understanding of phenomena connected with the novel boundary conditions (109). We validate our interpretation of the static and dynamic parts of the contact-angle conditions by demonstrating that the first term on the right-hand side of (109b) and (110b) determines the equilibrium contact angle, while the remaining terms on the right-hand side cause a dynamic delay in the attainment of this equilibrium contact angle value, see Experiments 1 and 2 below. Finally, we show that the dynamic terms have the potential to describe the phenomenon called dynamic contact angle hysteresis, see Experiment 3 below.

This section is structured in the following way. We first provide a dimensionless form of the governing equations. Then we proceed with identifying the corresponding continuous weak form and its discrete counterpart based on the Galerkin discretization. Finally, we briefly describe the numerical method and show the results of the three numerical experiments.

6.1 Dimensionless formulation

We introduce the same scaling as in Gomez et al. (2010). Each field quantity is expressed as $\varphi{=}[\varphi]\tilde{\varphi}$, where $[\varphi]$ denotes the scale of the quantity and $\tilde{\varphi}$ denotes the dimensionless counterpart. We introduce a spatial scale $[{\bf x}]{=}L_{0}$, and consider the scaling of spatial differential operators555This spatial scaling is clearly not optimal in the interfacial regions where another length scale corresponding to the thickness of the interfacial zone should probably be introduced. However, since we do not perform any scaling-based simplifications and the scaling only serves to provide dimensionless formulation, this issue can be ignored. $[\nabla]{=}L_{0}^{-1}$ and $[\hskip 0.85358pt\mathrm{div}]{=}L_{0}^{-1}$. We scale the density by $[\rho]{=}b^{\prime}$, where $b^{\prime}$ occurs as a parameter in the van der Waals model (89). The temperature is scaled by the critical temperature, i.e., $[\vartheta]{=}\vartheta_{c}$ (see (84)) and the pressure $p^{vdW}$ is scaled by $[p^{vdW}]{=}a^{\prime}(b^{\prime})^{2}$. For time, we pick the scale $[t]{=}L_{0}/\sqrt{a^{\prime}b^{\prime}}$ and the velocity is thus scaled by $[{\bf v}]{=}[{\bf x}]/[t]=\sqrt{a^{\prime}b^{\prime}}$. We introduce the dimensionless numbers

[TABLE]

and consequently, we can rewrite the system of balance equations (93) as

[TABLE]

Next, we introduce the dimensionless numbers

[TABLE]

and a dimensionless function $\mathcal{C}(\tilde{\vartheta})$ defined in (127) in the Appendix. Then the boundary conditions (109) read as follows:

[TABLE]

Model A1:

[TABLE]

• 2.

Model A2:

[TABLE]

• 3.

Model B1 and C:

[TABLE]

• 4.

Model B2:

[TABLE]

6.2 Weak formulations of the initial and boundary value problems and their numerical discretization

In this subsection, we first introduce the weak formulations to the system of governing equations (113)–(114) and then we present its discretization. For simplicity, we avoid using tildes in the dimensionless formulations (113) and (114). Since with respect to density, the strong form of the momentum balance (113b) involves the third derivative, we employ a mixed formulation by introducing

[TABLE]

as a new variable. We also assume that both the bulk and the shear viscosities are constant, meaning that $\mu{=}[\mu]$, $\lambda{=}[\lambda]$ (implying that $\tilde{\mu}{=}\tilde{\lambda}{=}1$).

In order to specify a weak solution to (113)–(115), we first introduce several standard function spaces: the Lebesgue space $(L^{2}(\Omega),(\cdot,\cdot)_{\Omega})$, the Sobolev space $W^{1,2}(\Omega)$, and its dual $\left(W^{1,2}(\Omega)\right)^{*}$ with the corresponding duality pairing $\langle\cdot,\cdot\rangle_{\Omega}$. We also set the space

[TABLE]

and introduce the spaces

[TABLE]

We say that ${\bf u}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(\rho,{\bf v},z)\in X$ is a weak solution to (113)–(115) if

[TABLE]

Next, we replace the integrands (2\mathbb{D}\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}})^{-}_{\tau} and \frac{\partial\rho^{-}}{\partial\mathrm{{\bf n}}_{\mbox{\scriptsize{\Gamma}}}} in $B_{\partial\Omega}$ by means of (114) and employ the identity $\frac{\mathcal{B}}{Re_{\mu}}{=}\mathcal{D}(Ca)^{2}$. This will generate four different forms:

[TABLE]

Let us note that in (118b) and (118d), we applied the integration by parts in the second terms on the right-hand sides. The system is supplemented with initial conditions for the density and the velocity:

[TABLE]

The weak formulation (117) and (118) is discretized in time by a simple $\Theta$-scheme and in space by the Galerkin method. Denoting (finite-element) discrete subspaces of $X$ and $Y$ by $X_{h}$ and $Y_{h}$, respectively, we define the discrete solution at the $n$-th time level as ${\bf u}_{h}^{n}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(\rho_{h}^{n},{\bf v}_{h}^{n},z_{h}^{n})\in X_{h}$ satisfying

[TABLE]

6.3 Numerical solution

The discrete system (119) is implemented by a finite-element method in the software package FEniCS (Alnaes et al., 2015) and for the discrete spaces, we choose continuous piece-wise polynomial approximations $X_{h}=P_{3}{\times}(P_{2})^{3}{\times}P_{3}$, where $P_{N}$ denotes polynomials of order $N$. We apply a structured mesh to a two-dimensional domain $\Omega$ with the aspect height-to-length ratio $1{:}3$. The mesh consists of $N_{1}{\times}N_{2}$ squares, each divided into 4 regular triangles. We apply a scaling of the capillary number $Ca$ based on the refinement methodology proposed by Gomez et al. (2010). Gomez et al. (2010) argue that since the realistic resolution of the diffuse interface zone in the Korteweg models is out of the scope of macroscopic models, it is reasonable to treat the capillary number in such cases as a regularizing parameter; it’s adjustment is based on the given spatial resolution of the model in such a way that the diffuse interface remains reasonably resolved. Based on this idea Gomez et al. (2010) introduce the parameterization $Ca{=}\frac{h}{L_{0}}$, where $h$ is the characteristic length scale of the spatial mesh, here defined as $h{=}\frac{L_{0}}{2\sqrt{N_{1}N_{2}}}$, with the length scale $L_{0}{=}1$ and $N_{1}{=}90$, $N_{2}{=}30$. Moreover, we also adopt the scaling of the Reynolds numbers from Gomez et al. (2010), setting $Re_{\mu}{=}Re_{\lambda}{=}2Ca^{-1}$. Being interested only in qualitative properties of the model, we set all but one of the remaining dimensionless numbers equal to one, i.e., we assign $\mathcal{A}{=}\mathcal{B}{=}\mathcal{D}{=}1$. The exception is the parameter $\mathcal{C}$ (depending on temperature), which governs the equilibrium contact angle, which we want to control quantitatively. For temperature $\vartheta{=}0.85\vartheta_{c}$, we get $\mathcal{C}(0.85){\doteq}\frac{25.5}{9\sqrt{3}Ca}$, see (128) in the Appendix.

Since the boundary conditions corresponding to Models B1, B2, and C represent a subclass of the boundary conditions for Models A1 and A2, we only consider the latter two models. The dimensionless boundary conditions in the considered setting simplify to (we omit tilde symbols for brevity):

[TABLE]

Model A1:

[TABLE]

• 2.

Model A2:

[TABLE]

**Experiment 1

**In this numerical experiment, we study the evolution of the Korteweg - van der Waals fluid in a two-dimensional container $\Omega$ in the absence of body forces, meaning that ${\bf b}{=}{\bf 0}$ in (117). The system is initially at rest (${\bf v}_{0}{=}{\bf 0}$). Consequently, the only driving mechanisms for its evolution are the boundary conditions (120) provided that $\varphi{\neq}\frac{\pi}{2}$ on the upper or lower part of the boundary. In order to isolate the effect of the novel contact angle condition and the generalized Navier slip condition from the traditional Navier slip boundary condition, we set $\alpha{=}0$.

In Fig. 4, we depict the evolution of the density distribution from the initial condition (top left). The system consists of a vapour in the Maxwell state $\rho_{v}^{M}$ in the right part of the domain (white) and liquid in the Maxwell state with density $\rho_{l}^{M}$ in the left part of the panel (grey) separated by a flat interface perpendicular to the boundary. In order to demonstrate that with the value of $\mathcal{C}$ as in (128), the parameter $\varphi$ controls (and equals) the equilibrium value of the contact angle of the fluid-vapor interface, we prescribe in all simulations in Experiment 1 the value of the static contact angle $\varphi$ on the top boundary and $\pi{-}\varphi$ on the bottom boundary. The reason is that in this case the equilibrium interfaces are particularly simple, being linear. In the snaphots of the simulation shown in Fig. 4, we employ only the equilibrium (static) part of the contact angle-condition, which means we consider $\beta{=}0$. The red contour denotes the interface between the liquid and water phases defined here by the density value $\frac{\rho_{v}^{M}+\rho_{l}^{M}}{2}$. The arrows depict the velocity field. All quantities are dimensionless and since (with the exception of the equilibrium contact angle) we are interested only in qualitative behavior of the model, we do not show any scales. In Fig. 5 we show the final states of three simulations, which differ only in the value of $\varphi$, considering $\varphi{=}\frac{\pi}{3}$, $\frac{\pi}{4}$, and $\frac{\pi}{6}$. For comparison, we plot also black dashed lines with the slope corresponding to the prescribed $\varphi$ and we observe very good agreement.

In order to study the effect of the dynamic part of the contact angle condition, in Fig. 6, we depict the time evolution of the interface based on the value of the parameter $\beta$ for Models A1 (top row) and A2 (bottom row) for a given static contact angle $\varphi{=}\frac{\pi}{3}$. Note that the parameter $\beta$ appears both in the contact angle condition and in the generalized Navier slip condition; see eq. (120). The case $\beta{=}0$ corresponds to the solely static contact angle condition while for $\beta{>}0$ additional dissipative surface mechanism is present. In the second case, the evolution of the interface and motion of the contact points lags behind the case with the static contact angle and this dynamic effect is stronger for Model A2 than for Model A1 and depends in both cases on the values of $\beta$. While for Model A1 there appears to be a saturation of the dynamic effect with respect to increasing values of $\beta$, for Model A2 the bigger the value of $\beta$, the stronger the dynamic effect. It is important to note that for all non-zero values of $\beta$, the final equilibrium configuration matches the case with $\beta{=}0$ as expected since the additional terms are of non-equilibrium nature and must vanish in the final equilibrium state.

**Experiment 2

**In the second numerical experiment, we study the spreading of a droplet of a Korteweg - van der Waals liquid (83) in contact with a wall under the action of gravity (i.e. for ${\bf b}{=}-{g}{\bf e}_{z}$; $g$ being the gravity acceleration and ${\bf e}_{z}$ the unit vector in the vertical direction). We again set $\alpha{=}0$. We consider the same set of dynamic and static boundary conditions as is the first experiment and we plot the same quantities. In particular, in Fig. 7, we show the time evolution in the case of a solely static contact angle and in Fig. 8, we depict the zoomed-in evolution of the interface between the liquid and vapor phases for different values of $\beta$ for Models A1 (top row) and A2 (bottom row), respectively. As in Experiment 1, we see that the introduction of the dynamic contact angle condition and associated surface traction term in the momentum balance leads to a delay in the evolution of the interface and the contact point when compared with solely static contact angle conditions. The equilibria are again the same for all models. Again, the dynamic effect is weaker for Model A1 compared to Model A2, and the dependence on the value of $\beta$ is the same as in Experiment 1 - we observe a saturation of the dynamic effect for higher values of $\beta$ for Model A1, while the effect appears to monotonously increase with the increasing value of $\beta$ for Model A2.

**Experiment 3

**In the last numerical experiment, we study the combined effect of the dynamic contact angle condition and generalized Navier-slip at the boundary. We consider the same geometry as in Experiment 2, only the initial condition is such that the droplet is positioned more to the left. The body force is prescribed as

[TABLE]

i.e., we consider a droplet sliding down an inclined slope (with an inclination $\Phi{=}\frac{\pi}{6}$), viewed from a coordinate system rotated such that its horizontal axis is aligned with the slope. The (dimensionless) friction parameter is set to $\alpha{=}50$ in all experiments. In Fig. 9, we plot several time snapshots of the evolution (for Model A2 and $\beta{=}100$). Interestingly, we observe a difference between the values of the contact angles between the advancing side (right) and the receding side (left) of the droplet.

In Fig. 10 we show how this effect depends on the values of the dynamic coefficient $\beta$ for the two Models A1 and A2. We can see that the observed phenomenon is clearly governed by the $\beta$ parameter and is rather insensitive to the type of the Model (A1 vs. A2). The bigger the value of $\beta$, the more pronounced the effect. These results are satisfactory in the sense that they provide a possible explanation of the dynamic contact angle hysteresis observed in nature (see e.g. Bormashenko, 2013), often attributed to pinning of the contact line. Here it results from dissipative processes within the interfacial zone between the phases; such an explanation corresponds to the ideas suggested recently in Makkonen (2017).

7 Conclusions

In this paper we have developed a thermodynamical framework to identify the boundary conditions for a class of Korteweg-type fluids. We exploited the tools of continuum thermodynamics stemming from the balance equations both in the bulk and at the boundary of the domain, which was treated as an interface between the body and its surroundings. Assuming the constitutive equations for the Helmholtz free energy in the bulk and at the boundary, we identified the surface and bulk entropy production mechanisms giving us a starting point for the formulation of the constitutive equations in the bulk and at the boundary.

For three types of surface Helmholtz free energy of various complexity, we derived a hierarchy of corresponding constitutive equations at the boundary. While some of the constitutive relations on the boundary took standard forms, in particular the in-surface Fourier heat flux, and the heat transmission conditions across the surface (Kapitza conditions), we obtained also two novel boundary conditions mutually coupled by a common parameter. The first one represents a nontrivial generalization of the Navier slip condition, relating the traction force at the boundary with the slip velocity and a novel dynamic term - either the trace of the divergence of the bulk velocity field, or the normal derivative of the normal velocity component. The second novel boundary condition was interpreted as a contact angle boundary condition for the Korteweg fluid model and it relates the normal derivative of density with two types of terms. The static terms arise from the surface Helmhotz free energy and characterize the value of the equilibrium contact angle attained after cessation of all motion in the fluid. The other terms are dynamic and dissipative in nature and involve either the trace of the divergence of the bulk velocity or the normal derivative of its normal component. These terms do not affect the equilibrium value of the contact angle and are only active in the dynamic situation when the fluid is flowing. It should be noted that in the literature it is possible to find alternative dynamic contact angle conditions involving the term $\frac{\partial\rho^{-}}{\partial t}$; see Jacqmin (2000). From the perpective of the derivation presented here, such boundary conditions would be recovered provided that we do not replace the term \frac{D_{\mbox{\scriptsize{\Gamma}}}\rho^{-}}{Dt} in (32) using the mass balance in the bulk. The approach presented here thus represents a possible generalization of such models.

Considering isothermal processes at a subcritical temperature admitting coexistence of liquid and gaseous phases, we then made the model explicit. We assumed that the Helmholtz free energy in the bulk corresponds to the Korteweg - van der Waals fluid, and the surface Helmholtz free energy reflects a simple characterization of the static contact angle. For this model, we derived explicit formulae for the contact angle condition and for the generalized Navier slip. The resulting model was implemented in the finite-element software package FEniCS. We studied the qualitative behavior of the Korteweg fluid model with the derived boundary conditions in three numerical experiments. The first two experiments confirmed the interpretation of the novel boundary conditions, namely we observed a time lag in the attainment of the static (equilibrium) value of the contact angle and also a lag in the motion of the contact line with an increasing amplitude of the novel dynamic terms. In the third experiment, we studied the sliding of a liquid droplet over an inclined plane, and we observed the so-called contact angle hysteresis, that is, a difference of the contact angle between the advancing and the receding side of the droplet. This phenomenon is often attributed to pinning of the contact line to irregularities on the surface; in our model it results from a dissipative process localized in the contact zone.

It should be noted that when constructing the constitutive relations, we constrained ourselves to linear relations for simplicity. A nonlinear generalization of our approach is possible. Here one could follow various thermodynamic approaches, such as the construction based on the maximization of the rate of entropy production (Rajagopal and Srinivasa, 2004) or by defining a suitable convex dissipation potential and deriving the constitutive response accordingly in the context of the so-called generalized standard materials (e.g. Halphen and Son Nguyen, 1975). It is also worth noting that while we considered just one particular member of the rich family of the so-called diffuse interface models - a Korteweg fluid - we are positive that the developed methodology could also be applied to other members of this class. The concept of a diffuse interface between two distinct subregions has been used since its origin at the turn of the 19th century for instance in the context of multicomponent materials (Cahn and Hilliard, 1958, 1959), in the classical theory of superconductivity (Landau and Ginzburg, 1965), and, more recently in the modeling of various natural phenomena such as foams (Fonseca et al., 2007), solidification (Kobayashi, 1994), phase transitions in solids (Fried and Gurtin, 1994), and glass formation (Řehoř et al., 2017), to name just a few of the plethora of applications. In all of these applications, generalizations of the boundary conditions and in particular, the dynamic contact angle conditions expressed as conditions for the normal derivative of the particular order parameter, should be possible following the methodology developed in this manuscript. In particular, this approach might play a key role in identification of suitable boundary conditions for viscoelastic rate type models with stress diffusion (Málek et al., 2018). Yet another generalization of our models could be obtained by relaxing the assumptions made on the structure of the surface Cauchy stress tensor; here we assumed that it is spherical (membrane model) and that the surface tension is constant. We conjecture that relaxing these assumptions would lead to the appearance of additional dynamic terms and admitting surface tension gradient would allow one to capture phenomena such as the Marangoni effect (Marangoni, 1871).

Appendix A Evaluation of the static angle function $\gamma_{0}(\vartheta)$

We provide an explicit evaluation of the parameter $\gamma_{0}$, which governs the static (equilibrium) part of the contact angle condition for Models A1 and A2, see (108). Let us recall the definition of $\gamma_{0}$:

[TABLE]

The crucial step is to evaluate the fraction $\frac{\sigma_{{lv}}}{\sigma}$, i.e. to find the relation between the surface tension of the liquid-vapor interface and the parameter $\sigma$ appearing at the gradient term in the bulk Helmholtz free energy of the Korteweg fluid (see (83)). Based on Diehl (2007) and Dreyer and Kraus (2010), for the model described by the Helmholtz free energy (83) and (84), it holds

[TABLE]

with the thermodynamic pressure $p_{vdW}$ and the chemical potential $\mu_{vdW}$ given by (86) and (87), respectively. We introduce dimensionless function $\tilde{c_{0}}(\tilde{\vartheta})$ as in Diehl (2007) through

[TABLE]

where $\rho_{c}$, $p_{c}$ and $\vartheta_{c}$ are the critical density, pressure and temperature, respectively, given by (89), and $\tilde{\vartheta}{=}\frac{\vartheta}{\vartheta_{c}}$ is the dimensionless temperature. The function $\tilde{c_{0}}(\tilde{\vartheta})$ can be approximated by the following expression (Diehl, 2007, p.37):

[TABLE]

which provides a good fit for $\tilde{\vartheta}{\in}\langle 0.6,1\rangle$. Applying the scaling and the definition of the capilary number from Section 6.1, we rewrite finally (108)1 as follows

[TABLE]

where we introduced the dimensionless function $C(\tilde{\vartheta})$ as follows

[TABLE]

Finally, for the value $\tilde{\vartheta}{=}0.85$, considered in our numerical simulations, we evaluate the Maxwell states numerically by solving (90): $\rho_{l}^{M}{=}0.6024$, $\rho_{v}^{M}{\doteq}0.1066$, and, consequently, from (125) and (126), we obtain $\tilde{r}(0.85){\doteq}0.0203$, $\tilde{c}_{0}(0.85){\doteq}0.5172$. This yields the value of $\mathcal{C}$ used in the numerical simulations in Section 6:

[TABLE]

Acknowledgements

J. Málek and O. Souček acknowledge support of the project 18-12719S financed by the Czech Science Foundation. M. Heida is financed by Deutsche Forschungsgemeinschaft (DFG) through Grant CRC 1114 “Scaling Cascades in Complex Systems”, Project C05 Effective models for materials and interfaces with multiple scales.