# Modelling of miscible liquids with the Korteweg stress

Ilya Kostin, Martine Marion, Rozenn Texier-Picard, Vitaly A Volpert

► To cite this version:

Ilya Kostin, Martine Marion, Rozenn Texier-Picard, Vitaly A Volpert. Modelling of miscible liquids with the Korteweg stress. ESAIM: Mathematical Modelling and Numerical Analysis, 2003, 37 (5), pp.741-753. 10.1051/m2an:2003042. hal-04028587

HAL Id: hal-04028587

https://hal.science/hal-04028587

Submitted on 24 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MODELLING OF MISCIBLE LIQUIDS WITH THE KORTEWEG STRESS

ILYA KOSTIN<fnref n="1" />, MARTINE MARION<fnref n="2" />, ROZENN TEXIER-PICARD<fnref n="3" /> AND
VITALY A. VOLPERT<fnref n="3" />

**Abstract.** When two miscible fluids, such as glycerol (glycerin) and water, are brought in contact, they immediately diffuse in each other. However if the diffusion is sufficiently slow, large concentration gradients exist during some time. They can lead to the appearance of an “effective interfacial tension”. To study these phenomena we use the mathematical model consisting of the diffusion equation with convective terms and of the Navier-Stokes equations with the Korteweg stress. We prove the global existence and uniqueness of the solution for the associated initial-boundary value problem in a two-dimensional bounded domain. We study the longtime behavior of the solution and show that it converges to the uniform composition distribution with zero velocity field. We also present numerical simulations of miscible drops and show how transient interfacial phenomena can change their shape.

**Mathematics Subject Classification.** 35K50, 76D05.

Received: November 4, 2002.

## 1. INTRODUCTION

In the case of immiscible liquids nonlocal molecular interaction can lead to appearance of capillary phenomena. One of the most well-known examples is related to the behavior of drops: a drop of oil in water becomes spherical because of the interfacial tension which acts to minimize its surface.

In the case of miscible liquids the nonlocal interaction can also exist and lead to some transient phenomena that, by analogy with the immiscible case, we call transient capillary phenomena. The difficulty to study them is related not only to their transient character but to the fact that they are very weak and can be hidden by other phenomena. For example, contrary to a drop of oil on the surface of water, a drop of water on the surface of glycerin, initially round, changes its shape in time (Fig. 1). Another example of “strange” behavior of miscible liquids is described in [4]. A drop of water with sugar in glycerin with the adjusted density splits in two drops, one of them going up and another one going down. In both cases the gravity can play a role but the explanation of these effects is not known.

To model miscible liquids, as in the examples above, we cannot use the approach often used in the immiscible case where it is supposed that the liquids are separated by an infinitely narrow interface. We have to consider not an interface but a transition zone with a finite width. Already in 1893 Van der Waals wrote a paper about

<footnote>*Keywords and phrases.* Miscible liquids, Korteweg stress, drops.</footnote>

<footnote n="1">¹ Université de Franche-Comté, UMR 6623 CNRS, 25030 Besançon, France. e-mail: kostin@pdmi.ras.ru</footnote>

<footnote n="2">² École Centrale de Lyon, UMR 5585 CNRS, 69134 Ecully Cedex, France. e-mail: Martine.Marion@ec-lyon.fr</footnote>

<footnote n="3">³ Université Lyon 1, UMR 5585 CNRS, 69622 Villeurbanne Cedex, France. e-mail: texier@maply.univ-lyon1.fr, volpert@maply.univ-lyon1.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Evolution of the drop of water on the surface of glycerin.</figcaption></figure>

binary systems with continuous variation of density [13] where he had introduced the free energy functional with the square of the density gradient. This term describes the nonlocal interaction and plays an important role in the theory of phase separation [3].

Korteweg was the first one to propose in 1901 that a nonuniform density (concentration or temperature) distribution can induce stresses and convection in a fluid [5]. To model this behavior he suggested to add a specific stress tensor in the Navier-Stokes equations. In the case of two miscible liquids brought in contact, this tensor acts like an effective interfacial tension that relaxes with time due to the mass diffusion. There is a number of more recent works where various approaches are used to derive or to justify such kind of models [1,2].

Under the assumptions that diffusive fluxes are sufficiently small, that the fluids are incompressible and have the same density and viscosity, following [14], we write the mass and momentum conservation equations in the form:

$$ \frac{\partial C}{\partial t} + u \cdot \nabla C = d\Delta C, \qquad (1.1) $$

$$ \frac{\partial u}{\partial t} + (u \cdot \nabla)u = -\nabla p + \nu \Delta u + \nabla \cdot T(C), \qquad (1.2) $$

$$ \nabla \cdot u = 0. \qquad (1.3) $$
Here $C$ is the concentration, $u$ is the velocity, $p$ is the pressure, $\nu$ the viscosity, $d$ is the coefficient of mass diffusion, $\Delta$ denotes the Laplace operator, $T(C) = (T_{ij})_{1\le i,j\le 2}$ is an additional stress tensor defined by

$$T_{11} = k \left( \frac{\partial C}{\partial x_2} \right)^2, \quad T_{12} = T_{21} = -k \frac{\partial C}{\partial x_1} \frac{\partial C}{\partial x_2}, \quad T_{22} = k \left( \frac{\partial C}{\partial x_1} \right)^2,$$

where $k$ is a nonnegative constant.

In this work we study some mathematical properties of this model, namely, the existence and uniqueness of solutions, and their longtime behavior. We show that the solution converges to the uniform composition distribution with zero velocity. We carry out numerical simulations to study the behavior of miscible drops. In particular we show that a drop, initially elliptic, becomes round because of the convection arising due to the transient interfacial tension. The contents of the paper are as follows. In Section 2 we prove the global existence and uniqueness of the solution, and study its asymptotic for large time. In Section 3 we present numerical simulations of miscible drops. Finally, conclusions are given in Section 4.

## 2. EXISTENCE AND UNIQUENESS

Let $\Omega$ be a bounded regular open set in $R^2$ with boundary $\Gamma$. We consider the problem: Find $C: \Omega \times R^+ \to R$, $u: \Omega \times R^+ \to R^2$, and $p: \Omega \times R^+ \to R$, such that

$$\frac{\partial C}{\partial t} + u \cdot \nabla C = d\Delta C, \qquad (2.1)$$

$$\frac{\partial u}{\partial t} + (u \cdot \nabla)u = -\nabla p + \nu \Delta u + \nabla \cdot T(C), \qquad (2.2)$$

$$\nabla \cdot u = 0. \qquad (2.3)$$

This system of equations is completed by the boundary conditions

$$\frac{\partial C}{\partial n} = 0, \quad u = 0 \quad \text{on } \Gamma \qquad (2.4)$$

and the initial conditions

$$C(x, 0) = C_0(x), \quad u(x, 0) = u_0(x), \quad x \in \Omega. \qquad (2.5)$$

Here $\nu > 0$ and $d > 0$ are given. The stress tensor $T(C) = (T_{ij})_{1\le i,j\le 2}$ is defined by

$$T_{11} = k \left( \frac{\partial C}{\partial x_2} \right)^2, \quad T_{12} = T_{21} = -k \frac{\partial C}{\partial x_1} \frac{\partial C}{\partial x_2}, \quad T_{22} = k \left( \frac{\partial C}{\partial x_1} \right)^2, \qquad (2.6)$$

and its divergence by

$$\nabla \cdot T(C) = \begin{pmatrix} \frac{\partial T_{11}}{\partial x_1} + \frac{\partial T_{12}}{\partial x_2} \\ \frac{\partial T_{21}}{\partial x_1} + \frac{\partial T_{22}}{\partial x_2} \end{pmatrix}.$$

In this section we aim to show that the above problem possesses a unique solution in some appropriate functional space. For that purpose, we consider

$$H_1 = H^1(\Omega), \quad V_1 = \left\{ B \in H^2(\Omega), \frac{\partial B}{\partial n} |_{\Gamma} = 0 \right\}$$
and also

$$ 
\begin{aligned}
H_2 &= \{v \in L^2(\Omega)^2, \nabla \cdot v = 0, v \cdot n = 0 \text{ on } \Gamma\}, \\
V_2 &= \{v \in H_0^1(\Omega)^2, \nabla \cdot v = 0\}.
\end{aligned}
 $$

We will denote by ( , ) and | · | the usual inner product and norm both in L²(Ω) and in L²(Ω)².

The following existence and uniqueness result holds.

**Theorem 2.1.** For any initial condition $(C_0, u_0)$ from $H_1 \times H_2$ there exists a unique solution of problem (2.1)-(2.5) which satisfies

$$ C \in L^2(0, \tau; V_1) \cap C(R^+, H_1), \quad u \in L^2(0, \tau; V_2) \cap C(R^+, H_2), \quad \forall \tau > 0. $$

*Proof.* The proof relies on the classical Faedo-Galerkin method (see for instance Lions [7], Temam [11, 12], Ladyzhenskaya [6]). We start with estimates for an exact smooth solution of (2.1)-(2.5). As usual, they will also be valid for the Galerkin approximation at least if we use a convenient basis (see below).

(i) We first multiply equation (2.1) by C and integrate over Ω. Using (2.4) and the Green formula, we find

$$ \frac{1}{2} \frac{d}{dt} |C|^2 + d|\nabla C|^2 + \int_{\Omega} (u \cdot \nabla C) C dx = 0. \quad (2.7) $$

Now since $\nabla \cdot u = 0$ and $u = 0$ on $\partial\Omega$, integration by parts gives that the last term in (2.7) vanishes. Therefore (2.7) enables us to say that

$$ C \text{ is bounded in } L^{\infty}(0, \tau; L^2(\Omega)), \quad \forall \tau > 0. \quad (2.8) $$

(ii) Next, let us multiply (2.1) by $-k\Delta C$ and integrate over $\Omega$. We find

$$ \frac{k}{2} \frac{d}{dt} |\nabla C|^2 + dk|\Delta C|^2 - k \int_{\Omega} (u \cdot \nabla C) \Delta C dx = 0. \quad (2.9) $$

Also, we multiply (2.2) by u and integrate over $\Omega$. This gives

$$ \frac{1}{2} \frac{d}{dt} |u|^2 + \nu |\nabla u|^2 + \int_{\Omega} (u \cdot \nabla) u \cdot u \, dx - \int_{\Omega} (\nabla \cdot T(C)) \cdot u \, dx = 0. \quad (2.10) $$

Here, as above,

$$ \int_{\Omega} (u \cdot \nabla) u \cdot u \, dx = 0. $$

Besides, recalling (2.6), we have

$$ 
\begin{aligned}
\frac{\partial T_{11}}{\partial x_1} + \frac{\partial T_{12}}{\partial x_2} &= k \frac{\partial C}{\partial x_2} \frac{\partial^2 C}{\partial x_1 \partial x_2} - k \frac{\partial C}{\partial x_1} \frac{\partial^2 C}{\partial x_2^2} \\
&= k \left( \frac{\partial C}{\partial x_2} \frac{\partial^2 C}{\partial x_1 \partial x_2} + \frac{\partial C}{\partial x_1} \frac{\partial^2 C}{\partial x_1^2} \right) - k \frac{\partial C}{\partial x_1} \Delta C \\
&= \frac{k}{2} \frac{\partial}{\partial x_1} |\nabla C|^2 - k \frac{\partial C}{\partial x_1} \Delta C.
\end{aligned}
 $$

A similar argument for the second component of the divergence term finally gives the equality

$$ \nabla \cdot T(C) = \frac{k}{2} \nabla |\nabla C|^2 - k \Delta C \nabla C. \quad (2.11) $$
Therefore, since $\nabla \cdot u = 0$, we have

$$ \int_{\Omega} (\nabla \cdot T(C)) \cdot u \, dx = -k \int_{\Omega} \Delta C \nabla C \cdot u \, dx. \qquad (2.12) $$

Coming back to (2.10), we conclude that

$$ \frac{1}{2} \frac{d}{dt} |u|^2 + \nu |\nabla u|^2 - k \int_{\Omega} \Delta C \nabla C \cdot u \, dx = 0. \qquad (2.13) $$

Adding (2.9) and (2.13) enables us to write

$$ \frac{1}{2} \frac{d}{dt} (k |\nabla C|^2 + |u|^2) + dk |\Delta C|^2 + \nu |\nabla u|^2 = 0. \qquad (2.14) $$

Combined with (2.8), this equality implies that

$$ C \text{ is bounded in } L^2(0, \tau; V_1) \cap L^\infty(0, \tau; H_1), \forall \tau > 0, \qquad (2.15) $$

$$ u \text{ is bounded in } L^2(0, \tau; V_2) \cap L^\infty(0, \tau; H_2), \forall \tau > 0. \qquad (2.16) $$

(iii) The last estimate concerns the time derivatives. In view of equation (2.1) and the Sobolev-Gagliardo-Nirenberg inequalities, we have

$$ \begin{aligned} \left| \frac{\partial C}{\partial t} \right| & \le d |\Delta C| + |u \cdot \nabla C| \\ & \le d |\Delta C| + \|u\|_{L^4(\Omega)^2} \|\nabla C\|_{L^4(\Omega)^2} \\ & \le d |\Delta C| + M |u|^{1/2} \|u\|_{H^1(\Omega)^2}^{1/2} \|C\|_{H^1(\Omega)}^{1/2} \|C\|_{H^2(\Omega)}^{1/2} \end{aligned} $$

with some $M > 0$. Next, (2.15) and (2.16) provide that

$$ C' = \frac{\partial C}{\partial t} \text{ is bounded in } L^2(0, \tau; L^2(\Omega)), \forall \tau > 0. \qquad (2.17) $$

For the $u$ equation, the additional term with respect to the Navier-Stokes equations is the last one. Here, $\nabla \cdot T(C)$ is a sum of expressions of the form $\lambda D_i(D_j C D_k C)$. Denoting by $V_2'$ the dual space of $V_2$, we have

$$ \begin{aligned} \|D_i(D_j C D_k C)\|_{V_2'} & \le \|D_j C D_k C\|_{L^2(\Omega)} \\ & \le \|D_j C\|_{L^4(\Omega)} \|D_k C\|_{L^4(\Omega)} \\ & \le M |D_j C|^{\frac{1}{2}} |D_k C|^{\frac{1}{2}} \|D_j C\|_{H^1(\Omega)}^{\frac{1}{2}} \|D_k C\|_{H^1(\Omega)}^{\frac{1}{2}}. \end{aligned} \qquad (2.18) $$

Using also some classical estimates for the Navier-Stokes equations (see for example Temam [12]), we conclude that

$$ u' = \frac{\partial u}{\partial t} \text{ is bounded in } L^2(0, \tau; V_2'), \forall \tau > 0. \qquad (2.19) $$

(iv) We now implement the Galerkin method using a “special” basis of $H_1$ (and of $V_1$) composed with the eigenvectors of the Laplace operator associated to the Neumann boundary condition, together with an arbitrary basis of $V_2$. The Galerkin procedure yields the finite-dimensional solutions $C_m$ and $u_m$ defined on some maximal interval of time $[0, T_m]$. It is easy to see that the estimates (2.15, 2.16, 2.17) and (2.19) apply to $C_m$ and $u_m$. Consequently the solutions are global in time ($T_m = +\infty$).
Also, the estimates provide the existence of subsequences, still denoted by $C_m$ and $u_m$, such that

$C_m \to C$ in $L^2(0, T; V_1)$ weakly,

$C_m \to C$ in $L^\infty(0, T; H_1)$ weak-star,

$C'_m \to C'$ in $L^2(0, T; V'_1)$ weakly,

and

$u_m \to u$ dans $L^2(0, T; V_2)$ weakly,

$u_m \to u$ dans $L^\infty(0, T; H_2)$ weak-star,

$u'_m \to u'$ dans $L^2(0, T; V'_2)$ weakly.

Applying some classical compactness theorems (see for example [7, 12]) we also obtain the strong convergence of the Galerkin approximations $(C_m, u_m)$ in the space $L^2(0, \tau; H^1(\Omega)) \times L^2(0, \tau; L^2(\Omega)^2)$. We can now pass to the limit in the finite-dimensional problems and obtain the existence of solutions for the system (2.1)-(2.5). The arguments are mostly classical, and we only consider the term related to $\nabla \cdot T(C)$. For the approximate solutions, similarly to (2.18), we see that $D_j C_m D_k C_m$ is bounded in $L^2(0, \tau \times \Omega)$, so that some subsequence converges weakly in $L^2(0, \tau \times \Omega)$ to some function $\chi_{jk} \in L^2(0, \tau \times \Omega)$. Using the strong convergence of $C_m$ in the space $L^2(0, \tau; H^1(\Omega))$, we infer immediately that $\chi_{jk} = D_j C D_k C$. From this we conclude easily that we can pass to the limit in the approximate problem, and we obtain the existence of a solution of (2.1)-(2.5).

(v) We now aim to prove the uniqueness of such a solution. Let $(C^1, u^1)$ and $(C^2, u^2)$ be two solutions of (2.1)-(2.5). We set $C = C^1 - C^2$ and $u = u^1 - u^2$. Using the identity (2.11), we see that $C$ and $u$ are solutions of the problem

$$ \frac{\partial C}{\partial t} - d\Delta C + u^1 \cdot \nabla C^1 - u^2 \cdot \nabla C^2 = 0 \qquad (2.20) $$

$$ \begin{aligned} \frac{\partial u}{\partial t} - \nu \Delta u + (u^1 \cdot \nabla) u^1 - (u^2 \cdot \nabla) u^2 + \nabla p^1 - \nabla p^2 = & k \nabla (|\nabla C^1|^2 - |\nabla C^2|^2) \\ & - k (\Delta C^1 \nabla C^1 - \Delta C^2 \nabla C^2). \end{aligned} \qquad (2.21) $$

Next, we multiply (2.20) by $-k\Delta C$, (2.21) by $u$, integrate over $\Omega$ and add the corresponding equalities. We obtain

$$ \frac{1}{2} \frac{d}{dt} (|u|^2 + k |\nabla C|^2) + kd |\Delta C|^2 + \nu |\nabla u|^2 = - \int_{\Omega} ((u \cdot \nabla) u^2) \cdot u \, dx + k \int_{\Omega} (-\Delta C^1 \nabla C \cdot u + \Delta C \nabla C \cdot u^1) \, dx. \quad (2.22) $$

The integrals in the right-hand side can be estimated thanks to interpolation inequalities and the Young inequality. We find:

$$ \left| \int_{\Omega} ((u \cdot \nabla) u^2) \cdot u \, dx \right| \le \sqrt{2} |u| |\nabla u| |\nabla u^2| \le \epsilon |\nabla u|^2 + C_{\epsilon} |u|^2 |\nabla u^2|^2. $$

$$ \begin{aligned} \left| \int_{\Omega} \Delta C^1 \nabla C \cdot u \, dx \right| & \le M |\Delta C^1| |\Delta C|^{\frac{1}{2}} |\nabla C|^{\frac{1}{2}} |\nabla u|^{\frac{1}{2}} |u|^{\frac{1}{2}} \\ & \le \epsilon (|\Delta C|^2 + |\nabla u|^2) + C_{\epsilon} |\Delta C^1|^2 |\nabla C| |u|, \end{aligned} $$

where $M$ and $\epsilon$ are positive numbers, and a constant $C_\epsilon > 0$ depends on $\epsilon$. Similarly,

$$ \begin{aligned} \left| \int_{\Omega} \Delta C \nabla C \cdot u^1 \, dx \right| & \le M |\Delta C|^{\frac{3}{2}} |\nabla C|^{\frac{1}{2}} |\nabla u^1|^{\frac{1}{2}} |u^1|^{\frac{1}{2}} \\ & \le \epsilon |\Delta C|^2 + C_{\epsilon} |\nabla C|^2 |\nabla u^1|^2 |u^1|^2. \end{aligned} $$
Now take $\epsilon$ such that $2\epsilon < \min(\nu, kd)$. We obtain:

$$ \begin{aligned} \frac{d}{dt} [|u|^2 + k|\nabla C|^2] &\le C_\epsilon (|u|^2|\nabla u|^2 + |\Delta C^1|^2 |\nabla C| |u| + |\nabla C|^2 |u^1|^2 |\nabla u^1|^2) \\ &\le C_\epsilon (|\nabla u^2|^2 + |\Delta C^1|^2 + |u^1|^2 |\nabla u^1|^2) [|u|^2 + k|\nabla C|^2]. \end{aligned} \quad (2.23) $$

Let us denote

$$ \phi(t) = |\nabla u^2(t)|^2 + |\Delta C^1(t)|^2 + |u^1(t)|^2 |\nabla u^1(t)|^2. $$

Since $u^1, u^2 \in L^2(0, \tau; H_0^1(\Omega)^2) \cap L^\infty(0, \tau; L^2(\Omega)^2)$ and $C^1, C^2 \in L^2(0, \tau; H^2(\Omega)) \cap L^\infty(0, \tau; H^1(\Omega))$ we have $\phi \in L^1([0, \tau])$. Hence the inequality (2.23) can be rewritten in the form

$$ \forall t \ge 0, \quad \frac{d}{dt} \left[ \exp \left( -C_{\epsilon} \int_{0}^{t} \phi(s) ds \right) (|u|^2 + k|\nabla C|^2) \right] \le 0. $$

It follows that

$$ \forall t \ge 0, \quad \exp \left( -C_{\epsilon} \int_{0}^{t} \phi(s) ds \right) (|u|^2 + k|\nabla C|^2) \le |u(0)|^2 + k|\nabla C(0)|^2. $$

Since the right-hand side vanishes, we have $C = u = 0$ which provides the uniqueness result. The proof of Theorem 2.1 is complete. $\square$

**Remark 2.2.** 1. From the proof of uniqueness follows the continuous dependence of the solution on the initial condition.

2. The equation for the concentration satisfies the maximum principle. Therefore, from the estimate $a \le C_0(x) \le b$ for the initial condition follows the same estimate for the solution.

At the end of this section we study the longtime behavior of solutions of problem (2.1)–(2.3). First of all we note that the mean value of C is conserved:

$$ \int_{\Omega} C(x, t) dx = \int_{\Omega} C_0(x) dx, \quad \forall t \ge 0. $$

This equality can be obtained by integration of (2.1) over $\Omega$.

**Theorem 2.3.** The solution $(C, u)$ of (2.1)–(2.3) converges to a constant solution as $t \to \infty$:

$$ \begin{aligned} C(x,t) &\to \int_{\Omega} C_0(x) dx \quad \text{in } H_1, \\ u(x,t) &\to 0 \quad \text{in } H_2. \end{aligned} $$

The convergence is exponential in time.

*Proof.* Let

$$ \bar{C}(x,t) = C(x,t) - \int_{\Omega} C(x,t) dx. $$

Then $\bar{C}$ satisfies the equation

$$ \frac{\partial \bar{C}}{\partial t} + u \cdot \nabla \bar{C} = d \Delta \bar{C}. $$

Multiplying it by $\bar{C}$ and integrating over $\Omega$, we obtain

$$ \frac{1}{2} \frac{d}{dt} |\bar{C}|^2 + d |\nabla \bar{C}|^2 = 0. $$
From this and (2.14), we conclude that

$$ \frac{1}{2} \frac{d}{dt} (|\bar{C}|^2 + k|\nabla \bar{C}|^2 + |u|^2) + d|\nabla \bar{C}|^2 + dk|\Delta \bar{C}|^2 + \nu|\nabla u|^2 = 0. \quad (2.24) $$

Now $u = 0$ on $\Gamma$, so that the Poincaré inequality yields that

$$ |\nabla u| \geq \alpha_1 |u| $$

with some positive $\alpha_1$. Also since $\int_{\Omega} \bar{C}(x, t)dx = 0$, it is known that there exists some $\alpha_2 > 0$ such that

$$ |\nabla \bar{C}| \geq \alpha_2 |\bar{C}|. $$

Finally we also have

$$ |\Delta \bar{C}| \geq \alpha_3 |\nabla \bar{C}|, $$

where $\alpha_3$ denotes the smallest eigenvalue of $-\Delta$ associated with the Neumann boundary condition for zero mean value functions. Let $\alpha = \min(\alpha_1, \alpha_2, \alpha_3)$. We infer from (2.24) that

$$ \frac{1}{2} \frac{d}{dt} (|\bar{C}|^2 + k|\nabla \bar{C}|^2 + |u|^2) + \alpha \min(d, \nu) (|\nabla \bar{C}|^2 + k|\Delta \bar{C}|^2 + |\nabla u|^2) \leq 0. $$

Thus, by integration in time,

$$ |\bar{C}(t)|^2 + k|\nabla \bar{C}(t)|^2 + |u(t)|^2 \leq e^{-2\alpha \min(d, \nu)t} (|\bar{C}(0)|^2 + k|\nabla \bar{C}(0)|^2 + |u(0)|^2). $$

Theorem 2.3 is proved. □

<figure><img src="image_3.png" /></figure>

## 3. NUMERICAL RESULTS

In this section we present numerical simulations of problem (2.1)–(2.5). We have proved in Theorem 2.2 that the composition $C$ converges to a uniform distribution in space, and the velocity $u$ to 0 as time goes to infinity. However, during some transient time, while the concentration gradients remain large they can induce convection and change the shape of the drop.

### 3.1. Numerical method

For the numerical study, we introduce a stream function/vorticity formulation of the problem. We define the stream function $\psi$ and the vorticity $\omega$ by:

$$ u_1 = \frac{\partial \psi}{\partial x_2}, \quad u_2 = -\frac{\partial \psi}{\partial x_1} $$

and

$$ \omega = -\Delta \psi = \frac{\partial u_2}{\partial x_1} - \frac{\partial u_1}{\partial x_2}. $$

The equations of motion take the form

$$ \frac{\partial \omega}{\partial t} = \nu \Delta \omega - \frac{\partial \psi}{\partial x_2} \frac{\partial \omega}{\partial x_1} + \frac{\partial \psi}{\partial x_1} \frac{\partial \omega}{\partial x_2} + \frac{1}{\rho} \left[ \frac{\partial}{\partial x_1} \left( \frac{\partial T_{21}}{\partial x_1} + \frac{\partial T_{22}}{\partial x_2} \right) - \frac{\partial}{\partial x_2} \left( \frac{\partial T_{11}}{\partial x_1} + \frac{\partial T_{12}}{\partial x_2} \right) \right], \quad (3.1) $$

$$ 0 = \Delta \psi + \omega. \quad (3.2) $$

The boundary condition $u = 0$ can be written as

$$ \nabla \psi|_{\partial \Omega} = 0. $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Evolution of the spherical drop: a) initial composition and b) difference between the composition after 2 s and the initial composition.</figcaption></figure>

In particular, $\psi$ is constant on the boundary of $\Omega$, and we suppose that this constant equals zero. For the vorticity, we introduce an artificial second order boundary condition, obtained from a Taylor expansion of the stream function in the neighbourhood of the boundary. For example, at the left boundary

$$ \omega_{0,j}^n = -\frac{1}{2}\omega_{1,j}^n - \frac{3\psi_{1,j}^n}{h^2}, $$

where the superscript $n$ corresponds to the time discretization, the subscripts $i, j$ to the space discretization, $h$ is the space step. To solve the problem numerically we use the finite difference alternating direction method with the implicit-explicit scheme. Numerical accuracy is verified by decreasing the time and space steps.

## 3.2. Evolution of a spherical drop

We first consider a spherical drop in a square domain with the dimensions $L_1 = L_2 = 3$ cm. The parameter $k$ is estimated experimentally for the system water/glycerol, using a spinning drop tensiometer [9, 10]. We take the following values of the parameters:

$$ \nu = 10^{-4} \text{ m}^2/\text{s}; \quad d = 2.5 \times 10^{-8} \text{ m}^2/\text{s}; \quad k = 10^{-9} \text{ m}^4/\text{s}^2. $$

The simulations show that the composition remains radially symmetric (Fig. 2), and there is a weak convection. If we decrease the space step, the convective motion weakens (Fig. 3). Therefore it is a pure numerical effect that becomes negligible if the space discretization is sufficiently good.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Maximum of the stream function in the 2.5 first seconds. From top to bottom the number of discretization points is 100, 200 and 300.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Maximum of the stream function in the first second for an elliptic drop. The plain curve corresponds to 100 × 200 discretization points, and $ht = 5 \cdot 10^{-5}$ s. The crosses correspond to 150 × 300 discretization points, and $ht = 2 \cdot 10^{-5}$ s.</figcaption></figure>

### 3.3. Evolution of an elliptic drop

Here we consider a rectangular domain with the dimensions $L_1 = L_2/2 = 3$ cm. The drop is initially elliptic. We observe that the convection induced by the concentration gradients is now significant (see Fig. 4). Moreover, it does not depend any longer on the discretization step.

As a consequence of these flows, the shape of the drop evolves to a spherical shape (see Fig. 5). The variation of concentration is much higher than in the previous case. A comparison with a purely diffusive case (with no Korteweg stress) shows that this evolution is mainly a consequence of this stress.
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Evolution of the elliptic drop: a) initial composition and b) composition after 10 s.</figcaption></figure>

## 3.4. Evolution of a long drop

Finally we consider the case of a long linear drop (jet). In the case of immiscible liquids and in three space dimensions, the drop would break into two or more droplets. Here the behaviour is somehow similar: after a few seconds, there is a constriction of the drop in its middle, and it seems to break into two parts. But as time grows, its shape becomes more spherical, as in the previous case (see Fig. 6).

## 4. CONCLUSIONS

To study behavior of miscible liquids we use the model which includes the diffusion equation with convective terms and the Navier-Stokes equations with the Korteweg stress describing the nonlocal interaction and resulting in a volume force similar to the capillary force for immiscible liquids. This model admits an energy functional that decreases in time if the Korteweg constant $k$ is positive. The existence of the Lyapunov functional determines some mathematical properties of the problem. It allows us to prove the global existence and uniqueness of the solution, and to study its longtime asymptotics.

We carry out numerical simulations of the problem to investigate behavior of miscible drops. For the physically realistic values of the parameters a two-dimensional drop initially elliptic becomes round because of the convection. This first stage of the process takes about several seconds of the physical time. After that the
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Evolution of the composition for a long drop: a) $t = 0$ s, b) $t = 60$ s and c) $t = 90$ s; left – composition, right – drop profile (level lines of the composition).</figcaption></figure>

convection decays, and the drop diffuses during several hours remaining round. Finally the system converges to the homogeneous in space equilibrium.

It is easy to verify that radially symmetric and plane composition distributions can exist without convection. For example, if the initial composition distribution is radially symmetric and the initial velocity is zero, then the composition will diffuse being radially symmetric, and the velocity will remain zero. Our conjecture is that any other composition distributions except for radially symmetric and plane ones cannot exist without convection.

An interesting behavior is observed for long drops or jets. In the beginning they form two elliptic drops which touch each other. This can recall the capillary instability where a jet splits into drops. However it is just a transient stage before the next one where they have a tendency to form a single drop.
As it is already mentioned in the introduction, transient capillary phenomena in miscible liquids are very weak. To study them, experiments in the International Space Station are planned [8, 9].

REFERENCES

[1] D.M. Anderson, G.B. McFadden and A.A. Wheeler, Diffuse interface methods in fluid mechanics. *Annu. Rev. Fluid Mech.* **30** (1998) 139–165.

[2] L.K. Antanovskii, Microscale theory of surface tension. *Phys. Rev. E* **54** (1996) 6285–6290.

[3] J.W. Cahn and J.E. Hilliard, Free energy of a nonuniform system. I. Interfacial Free Energy. *J. Chem. Phys.* **28** (1958) 258–267.

[4] D. Joseph and M. Renardy, *Fundamentals of two-fluid dynamics, Vol. II*. Springer, New York (1992).

[5] D.J. Korteweg, Sur la forme que prennent les équations du mouvement des fluides si l'on tient compte des forces capillaires causées par des variations de densité considérables mais connues et sur la théorie de la capillarité dans l'hypothèse d'une variation continue de la densité. *Arch. Néerl. Sci. Exactes Nat. Ser. II* **6** (1901) 1–24.

[6] O.A. Ladyzhenskaya, *Mathematical theory of viscous incompressible flow*. Gordon and Breach (1963).

[7] J.L. Lions, *Quelques méthodes de résolution des problèmes aux limites non linéaires*. Gauthier-Villars, Paris (1969).

[8] J. Pojman, N. Bessonov, R. Texier, V. Volpert and H. Wilke, Numerical simulations of transient interfacial phenomena in miscible fluids, in *Proceedings AIAA*, Reno, USA (January 2002).

[9] J. Pojman, Y. Chekanov, J. Masere, V. Volpert, T. Dumont and H. Wilke, Effective interfacial tension induced convection in miscible fluids, in *Proceedings of the 39th AIAA Aerospace Science Meeting*, Reno, USA (January 2001).

[10] P. Petitjeans, Une tension de surface pour les fluides miscibles. *C. R. Acad. Sci. Paris Sér. I Math.* **322** (1996) 673–679.

[11] R. Temam, *Navier-Stokes equations. Theory and numerical analysis*. North-Holland Publishing Co., Amsterdam–New York, Stud. Math. Appl. 2 (1979).

[12] R. Temam, *Navier-Stokes equations and nonlinear functional analysis*. SIAM (1983).

[13] J.S. Rowlinson, Translation of J.D. van der Waals’ “The thermodynamic theory of capillarity under hypothesis of a continuous variation of density”. *J. Statist. Phys.* **20** (1979) 197.

[14] V. Volpert, J. Pojman and R. Texier-Picard, Convection induced by composition gradients in miscible liquids. *C. R. Acad. Sci. Paris Sér. I Math.* **330** (2002) 353–358.