# Viscous stresses modelling for one-fluid approach in Newtonian two-phase flows

Adrien Toutant

► To cite this version:

Adrien Toutant. Viscous stresses modelling for one-fluid approach in Newtonian two-phase flows. Physics of Fluids, 2025, 37 (2), 10.1063/5.0257217. hal-04959597

HAL Id: hal-04959597

https://hal.science/hal-04959597v1

Submitted on 20 Feb 2025

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
Viscous stresses modelling for one-fluid approach in Newtonian two-phase flows

Adrien Toutant*

University of Perpignan

PROMES-CNRS (UPR 8521)

Tecnosud, Rambla de la Thermodynamique

66100 Perpignan, France

(Dated: January 28, 2025)

A crucial issue in two-phase flows is related to the modelling of viscous stresses. A new general viscous stress model for Newtonian two-phase flows is derived theoretically from the Navier-Stokes equations and jump conditions at the interface. The model is developed in the framework of the one-fluid approach. The only assumptions are incompressibility, no phase change and constant surface tension. It is argued that this new modeling allows a theoretical definition of the equivalent viscous stress in the vicinity of interfaces. Moreover, for computational fluid dynamics, it can be seen as the physical justification for using both harmonic and arithmetic averages to define the one-fluid viscosity.

The theoretical problem consists in modelling the equivalent viscous stress tensor, $\bar{S}$, in a volume of control, $\Omega$, containing the interface. Noting $\bar{D}$ the average strain tensor, we are looking for an operator $\mu$:

$$ \bar{S} = \mu(\bar{D}) \quad (1) $$

The operator $\mu$ can be seen as an equivalent viscosity. Assuming that the operator is linear, it involves a fourth-order tensor with 81 coefficients that are expected to be functions of the volume fraction, $\alpha_k$, and the dynamic viscosity $\mu_k$ of each phase $k$.

We consider incompressible Newtonian flows with no phase change and constant surface tension. $D_k$ and $S_k$ are respectively the strain tensor and the viscous stress tensor of phase $k$ [1]

$$ D_k = \frac{1}{2} (\nabla u_k + \nabla^T u_k) \quad (2) $$

$$ S_k = \mu_k (\nabla u_k + \nabla^T u_k) \quad (3) $$

where $u_k$ is the velocity of phase $k$.

To show that the modelling of viscous stresses corresponds to a subgrid closure problem, let's introduce the volume-average operator:

$$ \bar{\phi} = \frac{1}{\int_{\Omega} d\Omega} \int_{\Omega} \phi d\Omega \quad (4) $$

where $\phi$ is a one-fluid quantity. Because $\Omega$ contains the interface, $\phi$ can be discontinuous. It is a distribution that can be a scalar, a vector or a tensor. For example, the stress tensor is discontinuous at the interface. Due to the regularization of the volume-average operator, the equivalent viscous stress tensor is continuous ($C^0$). We also define the k-phase volume-average operator:

$$ \tilde{\phi}_k = \frac{1}{\int_{\Omega} \chi_k d\Omega} \int_{\Omega} \chi_k \phi_k d\Omega = \frac{\chi_k \tilde{\phi}_k}{\bar{\chi}_k} = \frac{\bar{\chi}_k \tilde{\phi}_k}{\alpha_k} \quad (5) $$

where $\chi_k$ is the phase indicator function (equal to 1 in phase $k$ and to 0 otherwise) and $\alpha_k$ is the volume average of the phase indicator function $\alpha_k = \bar{\chi}_k$.

The average strain tensor and the equivalent viscous stress tensor are respectively defined by:

$$ \begin{aligned} \bar{D} &\triangleq \sum_k \chi_k D_k \\ &= \frac{1}{2} \sum_k \chi_k (\nabla u_k + \nabla^T u_k) \end{aligned} \quad (6) $$

$$ \begin{aligned} \bar{S} &\triangleq \sum_k \chi_k S_k \\ &= \sum_k \chi_k \mu_k (\nabla u_k + \nabla^T u_k) \\ &= \sum_k \chi_k \mu_k \sum_k \chi_k (\nabla u_k + \nabla^T u_k) \\ &= 2\mu \bar{D} \end{aligned} \quad (7) $$

where $\mu$ the one-fluid viscosity and $D$ the one-fluid strain tensor (see the work of Ishii [2] and Kataoka [3] for an introduction to the one-fluid formulation). Because, the volume $\Omega$ is arbitrarily small, a naive approximation of the equivalent viscous stress tensor is $\bar{S} = 2\mu\bar{D} \approx 2\bar{\mu}\bar{D}$ with $\bar{\mu} = \alpha_1\tilde{\mu}_1 + \alpha_2\tilde{\mu}_2$. Unfortunately, because the quantities are not continuous, this is not true. However, the property is true for the phase volume-average operator (because single phase quantities are assumed to be continuous and $\Omega$ is arbitrarily small): $\tilde{\phi}_k\tilde{\phi}_k = \tilde{\phi}_k\tilde{\phi}_k$. Consequently, one gets $\bar{\phi}_k\bar{\phi}_k = \alpha_1^2\tilde{\phi}_1\tilde{\phi}_1 + \alpha_2^2\tilde{\phi}_2\tilde{\phi}_2 + \alpha_1\alpha_2(\tilde{\phi}_1\tilde{\phi}_2 + \tilde{\phi}_2\tilde{\phi}_1)$. Using this equality and the fact that $\alpha_1 + \alpha_2 = 1$, basic algebra gives the following general equation for interfacial subgrid terms:

$$ \bar{\phi}\bar{\phi} - \bar{\phi}\bar{\phi} = \alpha_1\alpha_2(\tilde{\phi}_1 - \tilde{\phi}_2)(\tilde{\phi}_1 - \tilde{\phi}_2) \quad (8) $$

Applying this equation to the subgrid term related to the
equivalent viscous stress tensor, one obtains:

$$ \tau_{SG} \hat{=} \mu \overline{\mathbf{D}} - \overline{\mu} \overline{\mathbf{D}} = \alpha_1 \alpha_2 (\tilde{\mu}_1 - \tilde{\mu}_2) (\tilde{\mathbf{D}}_1 - \tilde{\mathbf{D}}_2) \quad (9) $$

The same equation (excepted that the dynamic viscosities are assumed to be constant) is given in the PhD thesis of Ivan Girault [4]. The jump relation at the interface is written:

$$ [\phi] = \phi_1 - \phi_2 \quad (10) $$

At this step, it is required to determine $\tilde{\mathbf{D}}_1 - \tilde{\mathbf{D}}_2$. Because the volume $\Omega$ is arbitrarily small and crossed by the interface, this quantity corresponds in fact to the jump at the interface of the strain tensor:

$$ \lim_{|\Omega| \to 0} \tilde{\mathbf{D}}_1 - \tilde{\mathbf{D}}_2 = \mathbf{D}_1 - \mathbf{D}_2 = [\mathbf{D}] \quad (11) $$

Obviously, the components of the strain tensor that are continuous (respectively discontinuous) at the interface give a jump condition equal to zero (respectively not equal to zero). Without phase-change, the velocity is continuous at the interface and one has:

$$ [\mathbf{u}] = \mathbf{u}_1 - \mathbf{u}_2 = 0 \quad (12) $$

Consequently, the tangential derivatives of the velocity are continuous at the interface and one gets:

$$ [\mathbf{D} - ((\mathbf{n} \otimes \mathbf{n} \cdot \mathbf{D} + \mathbf{D} \cdot \mathbf{n} \otimes \mathbf{n}) - \mathbf{n} \otimes \mathbf{n} \cdot \mathbf{D} \cdot \mathbf{n} \otimes \mathbf{n})] = 0 \quad (13) $$

where $\mathbf{n}$ is the unit normal to the interface (the demonstration is written in appendix). The incompressibility constraint implies that the trace of the strain tensor is equal to zero in each phase $k$:

$$ \mathrm{tr}(\mathbf{D}_k) = 0 \quad (14) $$

Consequently, the jump of this quantity is also equal to zero:

$$ [\mathrm{tr}(\mathbf{D})] = 0 \quad (15) $$

This relation combined with the continuity of the tangential derivatives of the velocity implies that (the demonstration is written in appendix):

$$ [\mathbf{n} \otimes \mathbf{n} \cdot \mathbf{D} \cdot \mathbf{n} \otimes \mathbf{n}] = 0 \quad (16) $$

Let's introduce the following projection:

$$ P : \mathcal{M}_3(\mathbb{R}) \to \mathcal{M}_3(\mathbb{R}) \\ \mathbf{A} \mapsto (\mathbf{n} \otimes \mathbf{n} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{n} \otimes \mathbf{n}) \\ -2\mathbf{n} \otimes \mathbf{n} \cdot \mathbf{A} \cdot \mathbf{n} \otimes \mathbf{n} \quad (17) $$

The complementary projection is $P^{\perp} = \mathbf{I} - P$ where $\mathbf{I}$ is the identity function. Note that $P$ and $P^{\perp}$ are idempotent

(P ◦ P = P) and that P ◦ P$^{\perp}$ = 0. In particular, one gets the following decomposition for the equivalent stress tensor:

$$ \bar{\mathbf{S}} = P^{\perp}(\bar{\mathbf{S}}) + P(\bar{\mathbf{S}}) \quad (18) $$

Combining eq. 13 and eq. 16, one finds:

$$ [P^{\perp}(\mathbf{D})] = 0 \quad (19) $$

Coming back to the subgrid term related to the equivalent viscous stress tensor, it implies that $P^{\perp}(\tau_{SG}) = 0$ or

$$ P^{\perp}(\bar{\mathbf{S}}) = 2P^{\perp}(\overline{\mu \mathbf{D}}) = 2\overline{\mu}P^{\perp}(\bar{\mathbf{D}}) \quad (20) $$

Concerning the normal derivative of the tangential velocity components, one has to use the momentum jump condition [5]

$$ \sum_k (\mathbf{u}_k \dot{\mathbf{m}}_k - \mathbf{n}_k \cdot (-p_k \mathbf{I} + \mathbf{S}_k)) + (\nabla_s \sigma - (\nabla_s \cdot \mathbf{n})\sigma \mathbf{n})\delta_{\sigma} = 0 \quad (21) $$

where $\mathbf{I}$ is the identity function and $\mathbf{n}$ the unit normal at the interface. The $s$ means that differenciation operators are surface operators. $-\nabla_s \cdot \mathbf{n}$ is the mean curvature of the interface $\kappa$. $\sigma$ is the surface tension. $\delta_{\sigma}$ is the Dirac function associated to the interface. Without phase change $\dot{\mathbf{m}}_k = 0$ and with a constant surface tension $\nabla_s \sigma = 0$, the momentum jump condition becomes:

$$ \sum_k (\mathbf{n}_k \cdot (-p_k \mathbf{I} + \mathbf{S}_k)) + (\nabla_s \cdot \mathbf{n}) \sigma \mathbf{n} \delta_{\sigma} = 0 \quad (22) $$

The projection of this previous equation along each direction tangent to the interface implies:

$$ [P(\mathbf{S})] = 0 \quad (23) $$

That is equivalent to

$$ \tilde{\mu}_1 P(\tilde{\mathbf{D}}_1) = \tilde{\mu}_2 P(\tilde{\mathbf{D}}_2) \quad (24) $$

To come back to the subgrid term related to the equivalent viscous stress tensor, one combines eq. 9 and eq. 24. One finds:

$$ P(\bar{\mathbf{S}}) = 2\mu_h P(\bar{\mathbf{D}}) \quad (25) $$

where

$$ \mu_h = \frac{\tilde{\mu}_1 \tilde{\mu}_2}{\alpha_2 \tilde{\mu}_1 + \alpha_1 \tilde{\mu}_2} \quad (26) $$

Finally, the theoretical equivalent viscous stress tensor is given by eq. 18 combined with eq. 20 and eq. 25:

$$ \bar{\mathbf{S}} = 2 (\mu_a P^{\perp}(\bar{\mathbf{D}}) + \mu_h P(\bar{\mathbf{D}})) \quad (27a) $$
where

$$ \mu_a = \bar{\mu} = \alpha_1 \tilde{\mu}_1 + \alpha_2 \tilde{\mu}_2 \quad (27b) $$

$$ \mu_h = \frac{\tilde{\mu}_1 \tilde{\mu}_2}{\alpha_2 \tilde{\mu}_1 + \alpha_1 \tilde{\mu}_2} \quad (27c) $$

Assuming that the viscous stress tensor is linearly related to the deformation tensor, viscosity is generally a fourth-order tensor. Very general invariance properties, or others linked to the microstructure of the medium, drastically reduce the number of non-zero components of this tensor [6]. Ericksen [7], based on [8], determines the general form of the viscous stress tensor. For this, he assumes, among other things, that the mean fluid is anisotropic with a single preferred direction. In the case of a liquid-vapor interface and using the second gradient theory, Seppecher [9, 10] arrives at the same result

$$ \mathbf{S} = a (\mathbf{D} \cdot \mathbf{n} \otimes \mathbf{n} + \mathbf{n} \otimes \mathbf{n} \cdot \mathbf{D}) + b \mathbf{n} \otimes \mathbf{n} \cdot \mathbf{D} \cdot \mathbf{n} \otimes \mathbf{n} + c \mathbf{D} + (d \operatorname{tr}(\mathbf{D}) + e \mathbf{D} : \mathbf{n} \otimes \mathbf{n}) \mathbf{I} + e \operatorname{tr}(\mathbf{D}) \mathbf{n} \otimes \mathbf{n} \quad (28) $$

The coefficients a, b, c, d and e must satisfy the following conditions:

$$ c \geq 0 \quad (29a) $$

$$ d \geq 0 \quad (29b) $$

$$ a + c \geq 0 \quad (29c) $$

$$ b + c + d + 2e \geq 0 \quad (29d) $$

$$ d(b + c - 2a) - (a + e)^2 \geq 0 \quad (29e) $$

To compare the proposed closure (27) with the general form of the viscous stress tensor, we can rewrite it:

$$ \mathbf{S} = 2(\mu_h - \mu_a)(\mathbf{D} \cdot \mathbf{n} \otimes \mathbf{n} + \mathbf{n} \otimes \mathbf{n} \cdot \mathbf{D}) + 4(\mu_a - \mu_h)\mathbf{n} \otimes \mathbf{n} \cdot \mathbf{D} \cdot \mathbf{n} + 2\mu_a \mathbf{D} \quad (30) $$

Since we are working in incompressible flows, we have

$$ \operatorname{tr}(\mathbf{D}) = 0 \quad (31) $$

and the d coefficient can be any value. Consequently, only conditions (29a) and (29c) need to be considered. By identification, we have:

$$ a = 2(\mu_h - \mu_a) \quad (32a) $$

$$ b = 4(\mu_a - \mu_h) \quad (32b) $$

$$ c = 2\mu_a \quad (32c) $$

$$ e = 0 \quad (32d) $$

The conditions (29a) et (29c) correspond to $\mu_a \geq 0$ and $\mu_h \geq 0$, which is of course verified. Thus, the proposed closure is consistent with the general form of the viscous stress tensor.

The model given by eq. 27 has already been proposed and used by several authors in the case of constant dynamic viscosity (see for example [11–14]). However, for

the best of our knowledge, it was deduced from very simplified 2D cases. The present work provides the theoretical justification of the model in the general case.

In order to test the model, one uses the analytic case proposed by [15]. For a parallelepiped vertical channel in which the interface is parallel to one pair of lateral walls and assuming an infinite channel width in the y direction, the steady-state analytic solution is obtained by integrating the second order differential equation of the vertical velocity component in each phase k:

$$ \mu_k \frac{\partial^2 w_k}{\partial x^2} = \rho_k g \quad (33) $$

with no-slip boundary condition at the vertical channel walls ($w_1(0) = w_2(L) = 0$). Considering the interface is located at $x_I$, the analytical solution is:

$$ w_1 = \frac{\rho_1}{\mu_1} g \frac{x^2}{2} + c_1 x \quad \text{for } x \in [0, x_I] \quad (34a) $$

$$ w_2 = \frac{\rho_2}{\mu_2} g \frac{x^2 - L^2}{2} + c_2(x - L) \quad \text{for } x \in [x_I, L] \quad (34b) $$

with

$$ c_1 = g \frac{2\mu_1 x_I (L - x_I) (\rho_1 - \rho_2) + \rho_1 \mu_2 x_I^2 + \rho_2 \mu_1 (L^2 - x_I^2)}{2\mu_1 (\mu_1 (x_I - L) - \mu_2 x_I)} $$

$$ c_2 = g \frac{\mu_2 x_I^2 (2\rho_2 - \rho_1) + \rho_2 \mu_1 (L^2 - x_I^2)}{2\mu_2 (\mu_1 (x_I - L) - \mu_2 x_I)} $$

In this very simple case, $P^{\perp}(\bar{S}) = 0$ so $\bar{S} = 2\mu_h \bar{D}$. Figure 1 compares the results obtained with six velocity discretization points and a second-order centered scheme. The values of the parameters are $L = 1m$, $x_I = L/2$, $\mu_1 = 1 Pa.s$, $\mu_2 = 500 \mu_1$, $\rho_1 g = 1 N.m^{-3}$ and $\rho_2 = 100 \rho_1$. The model theoretically derived in this paper agrees very well with the analytical solution even if the number of grid points is small and the property jumps are large. The naive approximation ($\bar{S} = 2\mu_a \bar{D}$) fails to reproduce the analytical solutions.

It is worth noting that in the general case, arithmetic and harmonic viscosities must be used depending on the tensor component (see eq. 27).
<figure><img src="image_2.png" /><figcaption>FIG. 1. Velocity profiles obtained with the naive approximation of the equivalent viscous stress tensor (arithmetic viscosity) and the model theoretically derived in this paper (harmonic viscosity). The reference corresponds to the analytical solution.</figcaption></figure>

## APPENDIX

This section is dedicated to the demonstration of eq. 13, 16 and 23. Let's introduce the orthogonal natural basis ($e_1, e_2, n$) where $e_1$ and $e_2$ are the tangential directions associated to the main curvatures of the interface. The coordinates of this basis are noted ($\zeta^1, \zeta^2, \zeta^3$). Using these notations, eq. 12 implies:

$$ \left[ \frac{\partial \mathbf{u}}{\partial \zeta^1} \right] = \left[ \frac{\partial \mathbf{u}}{\partial \zeta^2} \right] = 0 \qquad (36) $$

In the orthogonal natural basis, the deformation tensor is given by:

$$ \mathbf{D} = \frac{1}{2} \begin{pmatrix} \frac{2}{h_1^2} \frac{\partial u_1}{\partial \zeta^1} & \frac{1}{h_1^2} \frac{\partial u_2}{\partial \zeta^1} + \frac{1}{h_2^2} \frac{\partial u_1}{\partial \zeta^2} & \frac{1}{h_1^2} \frac{\partial u_3}{\partial \zeta^1} + \frac{\partial u_1}{\partial \zeta^3} \\ \frac{1}{h_1^2} \frac{\partial u_2}{\partial \zeta^1} + \frac{1}{h_2^2} \frac{\partial u_1}{\partial \zeta^2} & \frac{2}{h_1^2} \frac{\partial u_2}{\partial \zeta^2} & \frac{1}{h_2^2} \frac{\partial u_3}{\partial \zeta^2} + \frac{\partial u_2}{\partial \zeta^3} \\ \frac{1}{h_1^2} \frac{\partial u_3}{\partial \zeta^1} + \frac{\partial u_1}{\partial \zeta^3} & \frac{1}{h_2^2} \frac{\partial u_3}{\partial \zeta^2} + \frac{\partial u_2}{\partial \zeta^3} & 2 \frac{\partial u_1}{\partial \zeta^3} \end{pmatrix} \qquad (37) $$

where

$$ h_1 = |1 - \kappa_1 \zeta^3| \qquad (38a) $$

$$ h_2 = |1 - \kappa_2 \zeta^3| \qquad (38b) $$

with $\kappa_1$ and $\kappa_2$ the main curvatures of the interface. Using eq. 36, one gets for the jump relation of the deformation tensor:

$$ [\mathbf{D}] = \frac{1}{2} \begin{pmatrix} 0 & 0 & \left[ \frac{\partial u_1}{\partial \zeta^3} \right] \\ 0 & 0 & \left[ \frac{\partial u_2}{\partial \zeta^3} \right] \\ \left[ \frac{\partial u_1}{\partial \zeta^3} \right] & \left[ \frac{\partial u_2}{\partial \zeta^3} \right] & 2 \left[ \frac{\partial u_3}{\partial \zeta^3} \right] \end{pmatrix} \quad (e_1, e_2, n) $$

Eq. 13 expresses the four components of the deformation tensor -where only tangential derivatives are involved- are continuous at the interface. Taking the trace of eq. 39, one obtains:

$$ [\mathrm{tr}(\mathbf{D})] = \left[ \frac{\partial u_3}{\partial \zeta^3} \right] \qquad (40) $$

Combined with the incompressible constraint, this last equation is equivalent to eq. 16. Using the natural orthogonal basis, the scalar product of the momentum jump relation (eq. 22) with the tangential directions gives:

$$ [\mathbf{n} \cdot \mathbf{S} \cdot \mathbf{e}_1] = 0 \qquad (41a) $$

$$ [\mathbf{n} \cdot \mathbf{S} \cdot \mathbf{e}_2] = 0 \qquad (41b) $$

The previous system is equivalent to eq. 23.

<footnote>* adrien.toutant@univ-perp.fr</footnote>

<footnote>[1] L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Pergamon, Oxford, 1987).</footnote>

<footnote>[2] M. Ishii, Thermo-fluid Dynamics Theory of Two-phase Flow (Eyrolles, 1975).</footnote>

<footnote>[3] I. Kataoka, Int. J. Multiphase Flow 12, 745 (1986).</footnote>

<footnote>[4] I. Girault, Formal and numerical developments towards Particle-Resolved Direct Numerical Simulation of hydrogen combustion in fluidized bed, Ph.D. thesis, INP Toulouse (2024).</footnote>

<footnote>[5] J. Delhaye, Int. J. Multiphase Flow 1, 395 (1974).</footnote>

<footnote>[6] A. C. Eringen, Continuum Physics, Vol. 2 (Academic Press, New York, 1975) pp. 131–172.</footnote>

<footnote>[7] J. L. Ericksen, Arch. Rational Mech. Anal. 4, 231 (1960).</footnote>

<footnote>[8] G. F. Smith and R. S. Rivlin, Quart. Appl. Math. 15, 308 (1957).</footnote>

<footnote>[9] P. Seppecher, Ph.D. thesis (1987).</footnote>

<footnote>[10] P. Seppecher, C. R. Acad. Sci. serie II, 497 (1989).</footnote>

<footnote>[11] D. B. Kothe, in Free Surface Flows, edited by H. C. Kuhlmann and H.-J. Rath (Springer Vienna, Vienna, 1998) pp. 267–331.</footnote>

<footnote>[12] A. Benkenida, Développement et validation d'une méthode de simulation d'écoulements diphasiques sans (e1, e2, n)reconstruction d'inter interfaces : application à la dynamique des bulles de Taylor, Ph.D. thesis, INP Toulouse (1999).</footnote>

<footnote>[13] A. Benkenida and J. Magnaudet, Comptes Rendus de l'Académie des Sciences - Series IIB - Mechanics-Physics-Astronomy 328, 25 (2000).</footnote>

<footnote>[14] A. Toutant, Modélisation physique des interactions entre interfaces et turbulence, Ph.D. thesis, INP Toulouse (2006).</footnote>

<footnote>[15] D. Garajeu and M. Medale, Commun. Comput. Phys. 24, 593 (2018).</footnote>