On lower-dimensional models of thin film flow,Part C: Derivation of a Reynolds type of equation for fluids with temperature and pressure dependent viscosity

Abstract

This paper constitutes the third part of a series of works on lower-dimensional models in lubrication. In Part A, it was shown that implicit constitutive theory must be used in the modelling of incompressible fluids with pressure-dependent viscosity and that it is not possible to obtain a lower-dimensional model for the pressure just by letting the film thickness go to zero, as in the proof of the classical Reynolds equation. In Part B, a new method for deriving lower-dimensional models of thin-film flow of fluids with pressure-dependent viscosity was presented. Here, in Part C, we also incorporate the energy equation so as to include fluids with both temperature and pressure dependent viscosity. By asymptotic analysis of this system, as the film thickness goes to zero, we derive a simplified model of the flow. We also carry out an asymptotic analysis of the boundary condition, in the case where the normal stress is specified on one part of the boundary and the velocity on the remaining part.

Nomenclature

$u,v,w$	Velocity components	m s ${}^{-1}$
$\mathbf{v}$	Velocity field $\mathbf{v}=u\mathbf{i}+v\mathbf{j}+w\mathbf{k}$	m s ${}^{-1}$
$p$	Mechanical pressure	Pa
$h$	Distance between the surfaces	m
${\mathbf{V}}_l$	Velocity of the lower surface, ${\mathbf{V}}_l=(u_l,v_l,0)$	m s ${}^{-1}$
${\mathbf{V}}_u$	Velocity of the upper surface, ${\mathbf{V}}_u=(0,0,w_u)$	m s ${}^{-1}$
${\mathbf{U}}_l$	Velocity of the lower surface in the $xy$ -plane, ${\mathbf{U}}_l=u_l\mathbf{i}+v_l\mathbf{j}$	m s ${}^{-1}$
$\rho$	Density	kg m ${}^{-3}$
$\mu$	Viscosity, $\mu =\mu (p,\theta )=f(\alpha p,\beta \theta )$	Pa s
$\theta$	Temperature	K
$e$	Internal energy	J
$\mathbf{q}$	Heat flux	W m ${}^{-2}$
$k$	Thermal conductivity	${\mathrm{W}\mathrm{m}}^{-2}{\mathrm{K}}^{-1}$
$c_v$	Specific heat	J K ${}^{-1}$
$\mathbf{T}$	Cauchy stress tensor	Pa
$\mathbf{g}$	Body force	m s ${}^{-2}$

Introduction

The celebrated lower-dimensional model, referred to as the Reynolds equation, that is appropriate to study of problems in lubrication was derived (in 1886) by Osborne Reynolds, ¹ under the assumption that the fluid film is thin, in the case of the incompressible Navier–Stokes fluid whose viscosity is constant. However, when the viscosity depends on the pressure, just the assumption that the film thickness is small is insufficient to obtain a lower-dimensional model similar to that of the Reynolds equation. In Part A, ² of the two-part paper,^2,3 Almqvist et al. have recently shown that by following the standard procedure used to reduce the dimension of the Navier–Stokes equation, one finds that there are additional terms that are related to the body forces and inertia that do not reduce in dimension appropriately. In Part B, ³ Almqvist et al. discussed an alternative procedure to obtain a lower-dimensional equation that is appropriate for lubrication problems by appealing to asymptotic analysis.

As an increase in pressure increases the viscosity and an increase in temperature decreases the viscosity, the competing effects between pressure and temperature help moderate changes in the viscosity. In this study, which should be recognised as Part C to the aforementioned two-part paper,^2,3 we address the competing effect between pressure and temperature as we develop a lower-dimensional model by applying asymptotic analysis when the viscosity depends on both the pressure and temperature. The flow can of course be modelled by the general system of equations which is obtained via balance of mass, balance of linear and angular momentum, and balance of energy, together with a set of constitutive relations. However, in realistic applications, it is difficult to solve this system of equations even with present day computers and software. Therefore, it is important to find simplified models which are accurate enough and can be efficiently handled from a numerical point of view. When the viscosity only depends on the pressure this has been done with success by using the fact that the fluid domain is very thin, i.e. that the gap between the surfaces is very small compared to the size of the surfaces. More precisely, one has by asymptotic analysis, as the film thickness goes to zero, derived lower-dimensional models (2-D models) for the pressure which are significantly easier to solve than the original system of equations, see Part A and B.^2,3

Seminal work in deriving simplified models of thin film flow by asymptotic analysis, as the film thickness goes to zero, was carried out by Osborne Reynolds. ¹ By rewriting the incompressible Navier–Stokes equation in a special dimensionless form and neglecting terms, which are small under the assumption that the lubricant film is very thin, he derived the following simplified model of lubricant flow:

$$\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}=\frac{1}{\mu }\frac{\mathrm{\partial }p}{\mathrm{\partial }x},\frac{{\mathrm{\partial }}^{2}v}{\mathrm{\partial }{y}^2}=\frac{1}{\mu }\frac{\mathrm{\partial }p}{\mathrm{\partial }y},\frac{\mathrm{\partial }p}{\mathrm{\partial }z}=0,$$

(1)

where $u$ and $v$ represent the $x$ and $y$ components of the fluid velocity field $\mathbf{v}=u\mathbf{i}+v\mathbf{j}+w\mathbf{k}$ . Moreover, by integrating the first two equations in (1) with respect to $z$ , inserting the result in the equation for conservation of mass, integrating across the fluid film and assuming no-slip boundary conditions at the boundaries Reynolds obtained a single equation for the mechanical pressure $p$ , namely, the so-called Reynolds equation

$$\frac{\mathrm{\partial }h}{\mathrm{\partial }t}=\mathrm{div}(\frac{h^3}{12\mu }\mathrm{\nabla }p-\frac{h}{2}\mathbf{U}),$$

(2)

where $h$ is the distance between the surfaces, $\mu$ is the viscosity and $\mathbf{U}$ is the sum of the velocities of the lower and upper surfaces in the $xy$ -plane. It is important to note that Reynolds derived (1) and (2) under quite restrictive assumptions. In particular, he assumed that the density $\rho$ is constant, the viscosity $\mu$ is constant and that the surfaces are rigid.

Mathematical models of lubrication are often very complex, but they are almost always built on (1) and (2). Unfortunately, the arguments which lead to these complex models are often conceptually wrong and there are numerous scientific works where the arguments leading to the model should be improved, corrected, or better justified. Let us give an example: In elastohydrodynamic lubrication (EHL) the pressure is extremely high which implies that the viscosity changes by several orders of magnitude in the fluid domain. In addition, the viscosity also depends on the temperature and the surfaces are elastically deformed due to the high pressure. All these effects are commonly taken into account by just replacing the constant viscosity in (2) by a viscosity-pressure-temperature relation $\mu (p,\theta )$ and say that the distance $h$ between the surfaces depends on $p$ , i.e. that $h=h(x,y,p)$ ; see, for example, the studies by Lugt and Morales-Espejel, ⁴ Spikes, ⁵ and Venner and Lubrecht. ⁶ It is obvious that such an approach must be questioned since (2) was derived under the assumption that the viscosity is constant and that the surfaces are rigid. It is not clear at all that one would obtain the same result if one in the full 3-D system assumed that the viscosity depends on the pressure.

In the two previous parts of this paper,^2,3 we pointed out common misinterpretations and incorrect usage of the Reynolds equation when the viscosity depends on the pressure. In particular, it was highlighted that to obtain a system of equations, which is suitable for asymptotic analysis one must start by applying the theory for implicit constitutive relations developed by Rajagopal,^7,8 which leads to a generalized Navier–Stokes equation. The reason is that it is conceptually wrong to say that the viscosity in the Navier–Stokes (Newtonian) constitutive model for an incompressible fluid depends on the pressure (see the studies by Almqvist et al.^2,3 and Rajagopal^7,8), as such a dependence makes the constitutive relation for the stress and the symmetric part of the velocity gradient to be implicit. That is, one no more has the stress defined explicitly in terms of the symmetric part of the velocity gradient.

It is well-known that in many realistic applications it is important to consider the fact that the viscosity depends on the temperature. The purpose of this paper is to generalize the results in Part A and B, by deriving a simplified system of equations governing incompressible thin film flow where the viscosity not only depends on the pressure, but also the temperature. In addition, we show that it is not trivial to determine the appropriate boundary conditions that should be used in the simplified model.

There are numerous works considering the derivation of simplified models of thin film flow. We do not give a complete review of all these results, but we only mention those which we judge to be most relevant for the present work. In the mathematically oriented paper ⁹ a generalized Reynolds equation and an associated energy equation is derived. However, the asymptotic analysis starts from the Stokes equation (i.e. inertial terms are neglected a priori), the viscosity only depends on the temperature and the gap between the surfaces does not change with time. The results in Ciuperca et al. ⁹ are consistent with the more general results presented in this paper. In the well-cited paper by Dowson ¹⁰ a generalized Reynolds equation was also derived. It permits the variation of the viscosity across, as well as along, the lubricant film. The generalized form of the Reynolds equation which was derived by Dowson ¹⁰ is the same as (54), except for that the derivation presented by Dowson does not admit the viscosity to depend on the pressure. The reasons for this are clearly explained in the previous two parts of this paper.^2,3 In Dowson ¹⁰ an approximate energy equation is also discussed. In particular, it is pointed out that the term in the energy equation which describes the viscous dissipation can be significantly reduced with little error. This is consistent with the energy equation (40) which is derived in this paper. The analysis in Dowson ¹⁰ is not restricted to incompressible fluids. However, one should be aware of that the analysis starts by inserting the Stokes assumption into the Navier–Stokes equation and nowadays it is known that Stokes assumption is inapt, see Rajagopal. ¹¹ The derivation of a generalized Reynolds equation for fluids with pressure dependent viscosity has also been considered in previous works,^12–16 with other techniques and assumptions.

The paper is organized in the following way: In Section “The 3-D mathematical model”, we present a mathematical model of isochoric flow for fluids with temperature and pressure dependent viscosity. In particular we show that, implicit constitutive theory must be used since the viscosity in an incompressible Newtonian fluid cannot depend on the pressure. In Section “Derivation of a simplified model of thin film flow”, we derive a simplified model of the flow by asymptotic analysis of the full model as the fluid film thickness goes to zero. From this simplified model we derive, in Section “A system of Reynolds type for fluids with pressure and temperature dependent viscosity”, a system of equations for the pressure and temperature. The system consists of a generalized form of the Reynolds equation and a simplified version of the energy equation. The boundary condition to be used for the pressure in the generalized form of the Reynolds equation depends on the type of boundary condition one has on the lateral surface in the full model. The relation between the boundary condition in the full model and the boundary condition that should be used in the generalized form of the Reynolds equation is presented, for some typical examples, in Section “Asymptotic analysis of the boundary conditions”. Finally, we present the main conclusions and some concluding remarks in Section “Conclusions and concluding remarks”.

The 3-D mathematical model

In this section we describe the full $3$ -D model which will be used as a starting point in the asymptotic analysis. In particular, it is highlighted that implicit constitutive theory must be used since it is conceptually wrong to just assume that the viscosity in the incompressible Navier–Stokes constitutive relation depends on the pressure and temperature.

Description of the “thin” flow domain

We consider flows where the fluid is confined between two surfaces. More precisely, the lower surface is a flat surface in the $xy$ -plane and is denoted by $\omega$ . Moreover, the lower surface is moving in the $xy$ -plane with the velocity ${\mathbf{V}}_l=(u_l,v_l,0)$ . The upper surface is $z=\epsilon h(x,y,t)$ , where $\epsilon >0$ is a small parameter and $h(x,y,t)$ is a smooth strictly positive function which is defined for $(x,y)\in \omega$ and for each time $t>0$ . This means that $\epsilon$ controls the thickness of the fluid domain and $h$ describes the shape of the upper surface, which we denote by ${\mathrm{\Gamma }}_u$ , and that the distance between the surfaces is $h_{\epsilon }=\epsilon h$ . The upper surface is only moving in the $z$ -direction with the velocity ${\mathbf{V}}_u=(0,0,w_u)$ . Note that this assumption implies that $w_u=\mathrm{\partial }{h}_{\epsilon }/\mathrm{\partial }t$ . Hence, the fluid domain for a fixed time $t$ is

$$\begin{array}{rl}& {\mathrm{\Omega }}_{\epsilon }(t)=\{(x,y,z)\in {\mathbb{R}}^3:(x,y)\in \omega \subset {\mathbb{R}}^2,\\ & 0<z<\epsilon h(x,y,t),t>0\},\end{array}$$

(3)

see Figure 1.

Figure 1.

An illustration of a typical fluid domain ${\mathrm{\Omega }}_{\epsilon }(t)$ (transparent blue), where the upper surface of the lower body is denoted by $\omega$ , the distance is $h_{\epsilon }=\epsilon h$ between the surfaces, and the velocity of the upper- and the lower bodies are, ${\mathbf{V}}_u$ and ${\mathbf{V}}_l$ , respectively.

Moreover, we define

$$Q_{\epsilon }=\{(x,y,z,t):(x,y,z)\in {\mathrm{\Omega }}_{\epsilon }(t),0<t<T\}.$$

The fluid model

The Navier–Stokes constitutive relation for an incompressible homogeneous fluid is based on the assumption that the Cauchy stress tensor $\mathbf{T}$ has the form

$$\mathbf{T}=-p\mathbf{I}+\mathit{\tau }(\theta ,\mathrm{\nabla }\mathbf{v}),$$

(4)

where $\theta$ is the temperature,

$$\mathrm{\nabla }\mathbf{v}=\left(\begin{array}{ccc}{\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }x}}& {\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }y}}& {\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }z}}\\ {\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }x}}& {\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }y}}& {\displaystyle \frac{\mathrm{\partial }v}{\mathrm{\partial }z}}\\ {\displaystyle \frac{\mathrm{\partial }w}{\mathrm{\partial }x}}& {\displaystyle \frac{\mathrm{\partial }w}{\mathrm{\partial }y}}& {\displaystyle \frac{\mathrm{\partial }w}{\mathrm{\partial }z}}\end{array}\right)$$

is the velocity gradient, and $p$ is the indeterminate part of the stress due to the constraint of incompressibility (cannot be constitutively specified) and $p$ is referred to as the mechanical pressure (or simply the pressure). On requiring $\tau$ to be linear in the components of $\mathrm{\nabla }\mathbf{v}$ and that $\mathit{\tau }$ is isotropic one obtains the classical incompressible Navier–Stokes fluid model, that allows the material moduli to depend arbitrarily on the temperature $\theta$ , namely

$$\mathbf{T}=-p\mathbf{I}+2\mu (\theta )\mathbf{D}(\mathbf{v}),$$

(5)

where $\mathbf{D}(\mathbf{v})$ is the symmetric part of the velocity gradient $\mathrm{\nabla }\mathbf{v}$ , and $\mu$ is the viscosity of the fluid. In an incompressible Navier–Stokes fluid $\mathrm{tr}(\mathbf{T})=-3p$ , since the constraint of incompressibility requires that $\mathrm{tr}(\mathbf{D}(\mathbf{v}))=\mathrm{div}\mathbf{v}=0.$ Fluids which may be modeled by (5) are referred to as incompressible Navier–Stokes fluids (or incompressible Newtonian fluids). The viscosity of a Navier–Stokes fluid cannot depend on the pressure. Because, then we would have that

$$\mathbf{T}=\frac{1}{3}\mathrm{tr}(\mathbf{T}\mathbf{)}\mathbf{I}+2\mu (-\frac{1}{3}\mathrm{tr}\left(\mathbf{T}\right),\theta )\mathbf{D}(\mathbf{v}),$$

(6)

which clearly contradicts that $\mathit{\tau }$ in (4) only depends on $\theta$ and $\mathrm{\nabla }\mathbf{v}$ . Unfortunately, there are many works, where the authors, without any justification, let the viscosity in the incompressible Navier–Stokes constitutive relation (5) depend on the pressure and merely substitute this modified form in the final lubrication approximation, which is obviously conceptually incorrect.

The way to justify constitutive relations for incompressible fluids where the material moduli may depend on the pressure, is to go beyond the class of the Navier–Stokes fluid model or even the more general Stokesian fluid model and consider the fluid as a sub-class of fluids described the implicit constitutive relations. The idea of using implicit constitutive relations was recently introduced by Rajagopal.^7,8 The main idea is to, instead of starting from (4), start with the following implicit constitutive assumption: $\mathbf{F}(\rho ,\theta ,\mathbf{T},\mathbf{D})=\mathbf{0}$ , where $\mathbf{F}$ is a tensor-valued function. Since $\mathbf{F}$ is an isotropic function in the case of an isotropic fluid it follows by standard representation results that:

$$\begin{array}{rl}\mathbf{F}(\mathbf{T},\mathbf{D})=& {\alpha }_0\mathbf{I}+{\alpha }_1\mathbf{T}+{\alpha }_2\mathbf{D}+{\alpha }_3{\mathbf{T}}^2+{\alpha }_4{\mathbf{D}}^2\\ & +{\alpha }_5(\mathbf{T}\mathbf{D}+\mathbf{D}\mathbf{T})+{\alpha }_6({\mathbf{T}}^2\mathbf{D}+{\mathbf{D}}^2\mathbf{T})\\ & +{\alpha }_7(\mathbf{T}{\mathbf{D}}^2+{\mathbf{D}}^2\mathbf{T})+{\alpha }_8({\mathbf{T}}^2{\mathbf{D}}^2+{\mathbf{D}}^2{\mathbf{T}}^2)=\mathbf{0},\end{array}$$

(7)

where the material moduli ${\alpha }_i$ , $i=0,\dots ,8$ , depend on $\rho ,$ $\theta$ and the invariants

$$\begin{array}{rl}& \mathrm{tr}\mathbf{T},\mathrm{tr}\mathbf{D},\mathrm{tr}{\mathbf{T}}^2,\mathrm{tr}{\mathbf{D}}^2,\mathrm{tr}{\mathbf{T}}^3,\mathrm{tr}{\mathbf{D}}^3,\mathrm{tr}\left(\mathbf{T}\mathbf{D}\right),\mathrm{tr}\left({\mathbf{T}}^2\mathbf{D}\right),\\ & \mathrm{tr}\left({\mathbf{D}}^2\mathbf{T}\right),\mathrm{tr}\left({\mathbf{D}}^2{\mathbf{T}}^2\right)\end{array}.$$

This class of implicit constitutive relations is much richer than the class of Navier–Stokes fluids (5). In particular, it captures the constitutive relation

$$\mathbf{T}=-p\mathbf{I}+2\mu (p,\theta )\mathbf{D}(\mathbf{v})$$

(8)

for incompressible fluids.

The viscosity-pressure-temperature relationship in (8) is determined empirically or by molecular modeling. In general, the viscosity increases with the pressure and decreases with the temperature. In our analysis we consider viscosity-pressure-temperature relations of the form

$$\mu =\mu (p,\theta )=f(\alpha p,\beta \theta ),$$

(9)

where $\alpha$ and $\beta$ are nonnegative constants. The function $f$ is assumed to be smooth, strictly increasing in the first argument, strictly decreasing in the second argument and $0<{\mu }_{min}\le f(\alpha p,\beta \theta )\le {\mu }_{max}$ for all $p$ and $\theta$ , where ${\mu }_{min}$ and ${\mu }_{max}$ are constants. Hence, the viscosity-pressure-temperature relationship is very general.

In this paper we consider incompressible homogeneous fluids described by implicit constitutive relations which are modeled by the constitutive relation (8). A generalization would be to consider fluids described by implicit constitutive relations which undergo isochoric motion in isothermal flows, but which could change its volume due to changes in temperature. The pioneering works considering this possibility, for classical linear viscous fluids, were done independently by Oberbeck ¹⁷ and Boussinesq. ¹⁸ Many attempts have been made to rigorously justify the Oberbeck–Boussinesq approximate equations. An overview can be found in Rajagopal, ¹⁹ and the case with pressure dependent viscosity is considered in Rajagopal. ²⁰

Balance equations

In the present work we consider homogeneous fluids which may be considered as incompressible, i.e. the density $\rho$ is constant. In the following, we use the subscript $\epsilon$ to indicate that the thickness of the fluid domain $Q_{\epsilon }$ is $h_{\epsilon }=\epsilon h$ , when we list the balance equations that govern the physics of the fluid flow.

Balance of mass for an incompressible fluid, leads to the velocity field being divergence free, i.e.

$$\mathrm{div}{\mathbf{v}}_{\epsilon }=\frac{\mathrm{\partial }{u}_{\epsilon }}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{v}_{\epsilon }}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{w}_{\epsilon }}{\mathrm{\partial }z}=0,\mathrm{in}{Q}_{\epsilon }.$$

(10)

Balance of linear momentum implies that

$$\rho \frac{D{\mathbf{v}}_{\epsilon }}{Dt}=\rho \mathbf{g}+\mathrm{div}{\mathit{T}}_{\epsilon },\mathrm{in}{Q}_{\epsilon },$$

(11)

where $\mathbf{g}=(g_x,g_y,g_z)$ is the body force acting per unit mass, $D/Dt$ denotes the material derivative, i.e.

$$\frac{D}{Dt}=\frac{\mathrm{\partial }}{\mathrm{\partial }t}+u_{\epsilon }\frac{\mathrm{\partial }}{\mathrm{\partial }x}+v_{\epsilon }\frac{\mathrm{\partial }}{\mathrm{\partial }y}+w_{\epsilon }\frac{\mathrm{\partial }}{\mathrm{\partial }z}.$$

Inserting the constitutive relation (8) into the equation (11) representing balance of linear momentum leads to the following generalized incompressible Navier–Stokes equation

$$\begin{array}{r}\rho \frac{D{\mathbf{v}}_{\epsilon }}{Dt}=\rho \mathbf{g}-\mathrm{\nabla }{p}_{\epsilon }+2\mathrm{div}(\mu (p_{\epsilon },{\theta }_{\epsilon })\mathbf{D}({\mathbf{v}}_{\epsilon })),\mathrm{in}{Q}_{\epsilon }\end{array}.$$

(12)

Balance of energy yields the equation

$$\rho \frac{D{e}_{\epsilon }}{Dt}=-\mathrm{div}{\mathbf{q}}_{\epsilon }+\mathit{T}:\mathrm{\nabla }{\mathbf{v}}_{\epsilon }=-\mathrm{div}{\mathbf{q}}_{\epsilon }+2\mu (p_{\epsilon },{\theta }_{\epsilon })\mathbf{D}({\mathbf{v}}_{\epsilon }):\mathbf{D}({\mathbf{v}}_{\epsilon }),$$

(13)

in $Q_{\epsilon }$ , where. $e_{\epsilon }$ is the internal energy per unit mass and ${\mathbf{q}}_{\epsilon }$ is the heat flux. By inserting the constitutive relations

$$e_{\epsilon }=c_v{\theta }_{\epsilon }\mathrm{and}{\mathbf{q}}_{\epsilon }=-k\mathrm{\nabla }{\theta }_{\epsilon },$$

(14)

where $c_v$ is the specific heat and $k$ is the thermal conductivity, into (13) we obtain the energy equation

$$\rho c_v\frac{D{\theta }_{\epsilon }}{Dt}=k\mathrm{\Delta }{\theta }_{\epsilon }+{\mathrm{\Psi }}_{\epsilon },\mathrm{in}{Q}_{\epsilon },$$

(15)

where

$$\begin{array}{rl}{\mathrm{\Psi }}_{\epsilon }& =2\mu (p_{\epsilon },{\theta }_{\epsilon })\mathbf{D}({\mathbf{v}}_{\epsilon }):\mathbf{D}({\mathbf{v}}_{\epsilon })=\\ & 2\mu (p_{\epsilon },{\theta }_{\epsilon })(D_{11}^2+D_{22}^2+D_{33}^2+2D_{12}^2+2D_{13}^2+2D_{23}^2)\\ & =2\mu (p_{\epsilon },{\theta }_{\epsilon })({\left(\frac{\mathrm{\partial }{u}_{\epsilon }}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }{v}_{\epsilon }}{\mathrm{\partial }y}\right)}^2+{\left(\frac{\mathrm{\partial }{w}_{\epsilon }}{\mathrm{\partial }z}\right)}^2\\ & +\frac{\mathrm{\partial }{u}_{\epsilon }}{\mathrm{\partial }y}\frac{\mathrm{\partial }{v}_{\epsilon }}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{u}_{\epsilon }}{\mathrm{\partial }z}\frac{\mathrm{\partial }{w}_{\epsilon }}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{v}_{\epsilon }}{\mathrm{\partial }z}\frac{\mathrm{\partial }{w}_{\epsilon }}{\mathrm{\partial }y})\\ & +\mu (p_{\epsilon },{\theta }_{\epsilon })({\left(\frac{\mathrm{\partial }{u}_{\epsilon }}{\mathrm{\partial }y}\right)}^2+{\left(\frac{\mathrm{\partial }{u}_{\epsilon }}{\mathrm{\partial }z}\right)}^2+{\left(\frac{\mathrm{\partial }{v}_{\epsilon }}{\mathrm{\partial }x}\right)}^2\\ & +{\left(\frac{\mathrm{\partial }{v}_{\epsilon }}{\mathrm{\partial }z}\right)}^2+{\left(\frac{\mathrm{\partial }{w}_{\epsilon }}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }{w}_{\epsilon }}{\mathrm{\partial }y}\right)}^2).\end{array}$$

The system of equations (10), (12) and (15) is closed (it has five equations and five unknowns) and governs the flow.

Boundary and initial conditions

The fluid domain ${\mathrm{\Omega }}_{\epsilon }(t)$ is defined in (3). The boundary, $\mathrm{\partial }{\mathrm{\Omega }}_{\epsilon }(t)$ (illustrated in Figure 2) of the fluid domain consists of three different parts, namely, the lower surface $\omega$ in the $xy$ -plane, the upper surface ${\mathrm{\Gamma }}_u(t)$ where $z=h_{\epsilon }(x,y,t)$ and the wall of the domain

$$\begin{array}{rl}& {\mathrm{\Gamma }}_w(t)=\{(x,y,z)\in {\mathbb{R}}^3:(x,y)\in \mathrm{\partial }\omega ,\\ & 0<\mathit{z}<\epsilon \mathit{h}(\mathit{x}\mathit{,}\mathit{y}\mathit{,}\mathit{t}\mathit{)},\mathit{t}>0\}.\end{array}$$

As described in Subsection “Description of the “thin” flow domain” the lower surface $\omega$ is moving in the $xy$ -plane with the velocity ${\mathbf{V}}_l=(u_l,v_l,0)$ and the upper surface ${\mathrm{\Gamma }}_u$ ( $z=\epsilon h(x,y,t)$ ) is only moving in the $z$ -direction with the velocity ${\mathbf{V}}_u=(0,0,w_u)$ . We consider the flow where no-slip boundary conditions can be assumed at the lower and upper surfaces. Thus ${\mathbf{v}}_{\epsilon }={\mathbf{V}}_l$ on $\omega$ and ${\mathbf{v}}_{\epsilon }={\mathbf{V}}_u$ on ${\mathrm{\Gamma }}_u$ , i.e.

$${\mathbf{v}}_{\epsilon }(x,y,0,t)={\mathbf{V}}_l=(u_l,v_l,0)$$

and

$${\mathbf{v}}_{\epsilon }(x,y,h(x,y,t),t)={\mathbf{V}}_u=(0,0,w_u).$$

We assume that the fluid sticks to the surfaces. This together with the assumption that the upper surface only moves in the $z$ -direction imply that $w_{\epsilon }=\mathrm{\partial }{h}_{\epsilon }/\mathrm{\partial }t$ on ${\mathrm{\Gamma }}_u$ .

Figure 2.

An illustration of a fluid domain (transparent blue), with fluid boundary ${\mathrm{\Gamma }}_w$ exemplified with ${\mathrm{\Gamma }}_{\mathrm{v}}$ (blue arrow) at the inlet and ${\mathrm{\Gamma }}_{\mathit{T}}$ (red arrows) at the other vertical walls, and its solid boundaries ${\mathrm{\Gamma }}_u$ and $\omega$ .

To have a mathematical well-posed problem, we must also impose some type of boundary condition on the wall of the bearing domain. There are several choices for this. However, from a physical point of view it is often not obvious which boundary condition one should choose to best describe the situation under study. In this work we consider the situation where the velocity ${\mathbf{v}}_{\epsilon }$ is given on one part of ${\mathrm{\Gamma }}_w$ and the normal stress ${\mathit{T}}_{\epsilon }\mathbf{n}$ is given on the rest of ${\mathrm{\Gamma }}_w$ . Indeed, let us divide ${\mathrm{\Gamma }}_w$ into two parts ${\mathrm{\Gamma }}_{\mathbf{v}}$ and ${\mathrm{\Gamma }}_{\mathit{T}}$ and assume that

$${\mathbf{v}}_{\epsilon }={\mathbf{v}}_b\mathrm{on}{\mathrm{\Gamma }}_{\mathbf{v}}\mathrm{and}{\mathit{T}}_{\epsilon }\mathbf{n}=-p_b\mathbf{n}\mathrm{on}{\mathrm{\Gamma }}_{\mathit{T}},$$

where ${\mathbf{v}}_b=(u_b,v_b,w_b)$ and $p_b$ are given and do not depend on $z$ . We consider the case when both the measure of ${\mathrm{\Gamma }}_{\mathbf{v}}$ and ${\mathrm{\Gamma }}_{\mathit{T}}$ are nonzero. In this case, both the pressure and velocity are unique. Moreover, we assume that the following initial condition

$$p_{\epsilon }=P_0,{\theta }_{\epsilon }={\mathrm{\Theta }}_0\mathrm{and}{\mathbf{v}}_{\epsilon }={\mathbf{V}}_0,\mathrm{at}t=0,$$

where $P_0,{\mathrm{\Theta }}_0$ and ${\mathbf{V}}_0$ are known smooth functions. We will in our analysis assume that the temperature is given at the surfaces, i.e.

$${\theta }_{\epsilon }(x,y,0,t)={\theta }_l\mathrm{and}{\theta }_{\epsilon }(x,y,h(x,y,t),t)={\theta }_u.$$

While this is just one example of boundary conditions, we want to point out that this is not by any means a universal condition and that the thermal boundary conditions (especially) need to be carefully chosen and matched to the application at hand.

Remark.

If the measure of ${\mathrm{\Gamma }}_{\mathit{T}}$ is zero, then ${\mathbf{v}}_b$ must satisfy the compatibility condition ${\int }_{{\mathrm{\Gamma }}_w}{\mathbf{v}}_b\cdot \mathbf{n}ds=0$ .

Derivation of a simplified model of thin film flow

To compute the solution of the system of balance equations (10), (12) and (15) with associated boundary- and initial conditions is very difficult. It becomes extra difficult when the domain is very thin, i.e. $\epsilon \ll 1$ . Therefore, it is desirable to find a simplified mathematical model, which describes the flow with reasonable accuracy, but is much simpler to solve numerically. We will by asymptotic analysis of (10), (12) and (15) as $\epsilon \to 0$ derive a simplified model of the flow, which may be used when $\epsilon \ll 1$ .

The main idea

The main idea for deriving is similar to the one that was used in Almqvist et al. ³ Indeed, assume that we want to compute the velocity, pressure and temperature for a given $\epsilon \ll 1$ , denoted by $\overline{\epsilon }.$ The strategy to find a simplified model of the flow takes into consideration the following:

We want to compute the solution $({\mathbf{v}}_{\overline{\epsilon }},p_{\overline{\epsilon }},{\theta }_{\overline{\epsilon }})$ to the system (10), (12) and (15) for a given $\overline{\epsilon }\ll 1$ .

Construct an auxiliary sequence $\{({\overline{\mathbf{v}}}_{\epsilon },{\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\}$ so that $({\overline{\mathbf{v}}}_{\overline{\epsilon }},{\overline{p}}_{\overline{\epsilon }},{\overline{\theta }}_{\overline{\epsilon }})$ coincides with the sought solution $({\mathbf{v}}_{\overline{\epsilon }},p_{\overline{\epsilon }},{\theta }_{\overline{\epsilon }})$ for $\epsilon =\overline{\epsilon }$ , i.e. $({\overline{\mathbf{v}}}_{\overline{\epsilon }},{\overline{p}}_{\overline{\epsilon }},{\overline{\theta }}_{\overline{\epsilon }})=({\mathbf{v}}_{\overline{\epsilon }},p_{\overline{\epsilon }},{\theta }_{\overline{\epsilon }})$ .

Show that ${\overline{\mathbf{v}}}_{\epsilon }\to {\overline{\mathbf{v}}}_0$ , ${\epsilon }^2{\overline{p}}_{\epsilon }\to {\overline{p}}_0$ and ${\overline{\theta }}_{\epsilon }\to {\overline{\theta }}_0$ , as $\epsilon \to 0,$ where $({\overline{\mathbf{v}}}_0,{\overline{p}}_0,{\overline{\theta }}_0)$ is the solution of some limit problem. The factor ${\epsilon }^2$ in front of ${\overline{p}}_{\epsilon }$ will be explained in Subsection “Asymptotic expansion of the velocity field, pressure and temperature”.

Use the fact that for small values of $\overline{\epsilon }$ it is reasonable to assume that the following approximations hold: ${\mathbf{v}}_{\overline{\epsilon }}={\overline{\mathbf{v}}}_{\overline{\epsilon }}\approx {\overline{\mathbf{v}}}_0$ , ${\overline{\epsilon }}^2{p}_{\overline{\epsilon }}={\overline{\epsilon }}^2{\overline{p}}_{\overline{\epsilon }}\approx {\overline{p}}_0$ and ${\theta }_{\overline{\epsilon }}={\overline{\theta }}_{\overline{\epsilon }}\approx {\overline{\theta }}_0$ .

Since $({\overline{\mathbf{v}}}_{\mathbf{0}},{\overline{p}}_0,{\overline{\theta }}_0)$ is the solution of the limit problem and $({\mathbf{v}}_{\overline{\epsilon }},{\overline{\epsilon }}^2{p}_{\overline{\epsilon }},{\theta }_{\overline{\epsilon }})\approx ({\overline{\mathbf{v}}}_0,{\overline{p}}_0,{\overline{\theta }}_0)$ the limit problem can be regarded as a simplified model of the flow.

Construction of the auxiliary sequence (scaling of $\mathit{\alpha }$ )

A crucial question is how to construct the auxiliary sequence $\{({\overline{\mathbf{v}}}_{\epsilon },{\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\}$ , described in Subsection “The main idea”. It is tempting to choose $\{({\overline{\mathbf{v}}}_{\epsilon },{\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\}$ as $\{({\mathbf{v}}_{\epsilon },p_{\epsilon },{\theta }_{\epsilon })\}$ . However, this is not possible. Roughly, speaking the problem is that: if the velocity of the lower surface is kept constant, then the pressure $p_{\epsilon }={\overline{p}}_{\epsilon }$ will go to infinity as $\epsilon \to 0$ . Hence, the viscosity will approach ${\mu }_{max}$ and therefore the limit equation will not reflect the fact that the viscosity depends on the pressure. This problem would not have appeared if the pressure-viscosity coefficient $\alpha$ was decreasing in an appropriate way at the same time as $\epsilon \to 0$ . Motivated by this we introduce an auxiliary viscosity-pressure-temperature relationship of the form

$$\overline{\mu }={\overline{\mu }}_{\epsilon }(p,\theta )=f(\overline{\alpha }{\epsilon }^{k}p,\overline{\beta }{\epsilon }^l\theta ),$$

(16)

where $k$ and $l$ are constants to be determined. Moreover, the constants $\overline{\alpha }$ and $\overline{\beta }$ are chosen such that $\overline{\alpha }{\overline{\epsilon }}^k=\alpha$ and $\overline{\beta }{\overline{\epsilon }}^l=\beta$ . In particular, this means that

$${\overline{\mu }}_{\overline{\epsilon }}(p,\theta )=f(\overline{\alpha }{\overline{\epsilon }}^{k}p,\overline{\beta }{\overline{\epsilon }}^l\theta )=f(\alpha p,\beta \theta )=\mu (p,\theta ).$$

Let us now define an auxiliary system of equations by replacing $\mu$ with ${\overline{\mu }}_{\epsilon }$ in the system (10), (12) and (15). We now chose the auxiliary sequence $\{({\overline{\mathbf{v}}}_{\epsilon },{\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\}$ , which we are looking for, see Subsection “The main idea”, as the solution of the auxiliary system. Indeed, $({\overline{\mathbf{v}}}_{\epsilon },{\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })$ is the solution of the auxiliary system:

$$\mathrm{div}{\overline{\mathbf{v}}}_{\epsilon }=0,\mathrm{in}{Q}_{\epsilon },$$

(17)

$$\rho \frac{D{\overline{\mathbf{v}}}_{\epsilon }}{Dt}=\rho \mathbf{g}-\mathrm{\nabla }{\overline{p}}_{\epsilon }+2\mathrm{div}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\mathbf{D}({\overline{\mathbf{v}}}_{\epsilon })),\mathrm{in}{Q}_{\epsilon },$$

(18)

$$\rho c_v\frac{D{\overline{\theta }}_{\epsilon }}{Dt}=k\mathrm{\Delta }{\overline{\theta }}_{\epsilon }+{\overline{\mathrm{\Psi }}}_{\epsilon },\mathrm{in}{Q}_{\epsilon },$$

(19)

where

$$\begin{array}{rl}\overline{\mathrm{\Psi }}& =2{\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\mathbf{D}({\overline{\mathbf{v}}}_{\epsilon }):\mathbf{D}({\overline{\mathbf{v}}}_{\epsilon })=2{\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\\ & (D_{11}^2+D_{22}^2+D_{33}^2+2D_{12}^2+2D_{13}^2+2D_{23}^2)\\ & =2{\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })({\left(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }y}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }z}\right)}^2\\ & +\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y}\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }z}\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }z}\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y})\\ & +{\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })({\left(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }z}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x}\right)}^2\\ & +{\left(\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }z}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y}\right)}^2).\end{array}$$

We note that for $\epsilon =\overline{\epsilon }$ the solution $({\mathbf{v}}_{\overline{\epsilon }},p_{\overline{\epsilon }},{\theta }_{\overline{\epsilon }})$ to the physical system and the solution $({\overline{\mathbf{v}}}_{\overline{\epsilon }},{\overline{p}}_{\overline{\epsilon }},{\overline{\theta }}_{\overline{\epsilon }})$ to the auxiliary system are the same. This means that if we can find a good approximation of $({\overline{\mathbf{v}}}_{\overline{\epsilon }},{\overline{p}}_{\overline{\epsilon }},{\overline{\theta }}_{\overline{\epsilon }})$ , then we have also found a good approximation of the solution $({\mathbf{v}}_{\overline{\epsilon }},p_{\overline{\epsilon }},{\theta }_{\overline{\epsilon }})$ to the physical problem.

Remark.

We deliberately do not discuss the appropriate boundary conditions for the auxiliary system (17)–(19) here. This discussion is postponed to Section “Asymptotic analysis of the boundary conditions”.

Transformation of the $\epsilon$ -dependent domain $Q_{\epsilon }$ into a fixed domain

The fluid domain ${\mathrm{\Omega }}_{\epsilon }(t)$ depends on the parameter $\epsilon$ . This means that ${\mathbf{v}}_{\epsilon }$ and $p_{\epsilon }$ as functions of $(x,y,z,t)$ are defined on different domains $Q_{\epsilon }$ for different values of $\epsilon$ . However, by a simple change of variables it is possible to convert the $\epsilon$ -dependent fluid domain ${\mathrm{\Omega }}_{\epsilon }(t)$ into a fixed domain $\mathrm{\Omega }(t)$ (Figure 3), and equivalently $Q_{\epsilon }$ into a fixed domain $Q$ . Indeed, if we define $Z=z/\epsilon$ , then the new variables are $(x,y,Z,t)$ and for each $\epsilon$ the domain $Q_{\epsilon }$ becomes

$$\begin{array}{rl}& Q=\{(x,y,Z,t):(x,y)\in \omega \subset {\mathbb{R}}^2,\\ & 0<Z<h(x,y,t),t>0\}.\end{array}$$

Figure 3.

An illustration of the $\epsilon$ -dependent fluid domain ${\mathrm{\Omega }}_{\epsilon }(t)$ (left) and corresponding the fixed domain $\mathrm{\Omega }(t)$ (right).

In the new variables, i) the equation (10) describing conservation of mass becomes

$$\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }y}+\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }Z}=0,$$

(20)

ii) the equation for balance of linear momentum (12) in component form becomes

$$\begin{array}{rl}& {\epsilon }^2\rho (\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }t}+{\overline{u}}_{\epsilon }\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }x}+{\overline{v}}_{\epsilon }\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y}+\frac{{\overline{w}}_{\epsilon }}{\epsilon }\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }Z})\\ & ={\epsilon }^2\rho g_x-{\epsilon }^2\frac{\mathrm{\partial }{\overline{p}}_{\epsilon }}{\mathrm{\partial }x}+{\epsilon }^{2}2\frac{\mathrm{\partial }}{\mathrm{\partial }x}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }x})\\ & +{\epsilon }^2\frac{\mathrm{\partial }}{\mathrm{\partial }y}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x}))\\ & +\frac{\mathrm{\partial }}{\mathrm{\partial }Z}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }Z}+\epsilon \frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x})),\end{array}$$

(21)

$$\begin{array}{rl}& {\epsilon }^2\rho (\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }t}+{\overline{u}}_{\epsilon }\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x}+{\overline{v}}_{\epsilon }\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }y}+\frac{{\overline{w}}_{\epsilon }}{\epsilon }\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }Z})\\ & ={\epsilon }^2\rho g_y-{\epsilon }^2\frac{\mathrm{\partial }{\overline{p}}_{\epsilon }}{\mathrm{\partial }y}+{\epsilon }^2\frac{\mathrm{\partial }}{\mathrm{\partial }x}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x}))\\ & +{\epsilon }^{2}2\frac{\mathrm{\partial }}{\mathrm{\partial }y}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }y})\\ & +\frac{\mathrm{\partial }}{\mathrm{\partial }Z}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })(\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }Z}+\epsilon \frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y})),\end{array}$$

(22)

$$\begin{array}{rl}& {\epsilon }^2\rho (\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }t}+{\overline{u}}_{\epsilon }\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x}+{\overline{v}}_{\epsilon }\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y}+\frac{{\overline{w}}_{\epsilon }}{\epsilon }\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }Z})\\ & ={\epsilon }^2\rho g_z-\epsilon \frac{\mathrm{\partial }{\overline{p}}_{\epsilon }}{\mathrm{\partial }Z}\\ & +\frac{\mathrm{\partial }}{\mathrm{\partial }x}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })(\epsilon \frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }Z}+{\epsilon }^2\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x}))\\ & +\frac{\mathrm{\partial }}{\mathrm{\partial }y}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })(\epsilon \frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }Z}+{\epsilon }^2\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y}))\\ & +2\frac{\mathrm{\partial }}{\mathrm{\partial }Z}({\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }Z}),\end{array}$$

(23)

and, iii) the energy equation (15) becomes

$$\rho c_v(\frac{\mathrm{\partial }{\overline{\theta }}_{\epsilon }}{\mathrm{\partial }t}+{\overline{u}}_{\epsilon }\frac{\mathrm{\partial }{\overline{\theta }}_{\epsilon }}{\mathrm{\partial }x}+{\overline{v}}_{\epsilon }\frac{\mathrm{\partial }{\overline{\theta }}_{\epsilon }}{\mathrm{\partial }y}+\frac{{\overline{w}}_{\epsilon }}{\epsilon }\frac{\mathrm{\partial }{\overline{\theta }}_{\epsilon }}{\mathrm{\partial }Z})=k(\frac{{\mathrm{\partial }}^2{\overline{\theta }}_{\epsilon }}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{\overline{\theta }}_{\epsilon }}{\mathrm{\partial }{y}^2}+\frac{1}{{\epsilon }^2}\frac{{\mathrm{\partial }}^2{\overline{\theta }}_{\epsilon }}{\mathrm{\partial }{Z}^2})+{\overline{\mathrm{\Psi }}}_{\epsilon },$$

(24)

where

$$\begin{array}{rl}& {\overline{\mathrm{\Psi }}}_{\epsilon }=2{\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })({\left(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }y}\right)}^2+\frac{1}{{\epsilon }^2}{\left(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }Z}\right)}^2\\ & +\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y}\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x}+\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }Z}\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x}+\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }Z}\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y})\\ & +{\overline{\mu }}_{\epsilon }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })({\left(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y}\right)}^2+\frac{1}{{\epsilon }^2}{\left(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }Z}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x}\right)}^2\\ & +\frac{1}{{\epsilon }^2}{\left(\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }Z}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y}\right)}^2).\end{array}$$

(25)

Asymptotic expansion of the velocity field, pressure and temperature

Let us assume that ${\overline{\mathbf{v}}}_{\epsilon }=({\overline{u}}_{\epsilon }$ , ${\overline{v}}_{\epsilon }$ , ${\overline{w}}_{\epsilon })$ , ${\overline{p}}_{\epsilon }$ and ${\overline{\theta }}_{\epsilon }$ have asymptotic expansions of the form

$$\begin{array}{rl}{\overline{u}}_{\epsilon }(x,y,Z,t)& ={\epsilon }^{m_u}({\overline{u}}_0(x,y,Z,t)+\epsilon {\overline{u}}_1(x,y,Z,t)\\ & +{\epsilon }^2{\overline{u}}_2(x,y,Z,t)+\cdots ),\end{array}$$

(26)

$$\begin{array}{rl}{\overline{v}}_{\epsilon }(x,y,Z,t)& ={\epsilon }^{m_v}({\overline{v}}_0(x,y,Z,t)+\epsilon {\overline{v}}_1(x,y,Z,t)\\ & +{\epsilon }^2{\overline{v}}_2(x,y,Z,t)+\cdots ),\end{array}$$

(27)

$$\begin{array}{rl}{\overline{w}}_{\epsilon }(x,y,Z,t)& ={\epsilon }^{m_w}({\overline{w}}_0(x,y,Z,t)+\epsilon {\overline{w}}_1(x,y,Z,t)\\ & +{\epsilon }^2{\overline{w}}_2(x,y,Z,t)+\cdots ),\end{array}$$

(28)

$$\begin{array}{rl}{\overline{p}}_{\epsilon }(x,y,Z,t)& ={\epsilon }^{m_p}({\overline{p}}_0(x,y,Z,t)+\epsilon {\overline{p}}_1(x,y,Z,t)\\ & +{\epsilon }^2{\overline{p}}_2(x,y,Z,t)+\cdots ),\end{array}$$

(29)

$$\begin{array}{rl}{\overline{\theta }}_{\epsilon }(x,y,Z,t)& ={\epsilon }^{m_{\theta }}({\overline{\theta }}_0(x,y,Z,t)+\epsilon {\overline{\theta }}_1(x,y,Z,t)\\ & +{\epsilon }^2{\overline{\theta }}_2(x,y,Z,t)+\cdots ),\end{array}$$

(30)

where none of the leading terms ${\overline{u}}_0,{\overline{v}}_0,{\overline{w}}_0,{\overline{p}}_0$ , and ${\overline{\theta }}_0$ are identically equal to zero and $m_u,m_v,m_w,m_p$ and $m_{\theta }$ are to be determined.

The no-slip boundary condition at the upper and lower surfaces motivates the assumption that $m_u=m_v=0$ . In the same way, the condition that the temperature is given at the surfaces justifies that $m_{\theta }=0$ . To determine $m_w$ we will use the equation (20) representing balance of mass. Indeed, by inserting the expansions (26)–(28) into (20) gives

$$\begin{array}{rl}& (\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }x}+\epsilon \frac{\mathrm{\partial }{\overline{u}}_1}{\mathrm{\partial }x}+\cdots )+(\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }y}+\epsilon \frac{\mathrm{\partial }{\overline{v}}_1}{\mathrm{\partial }y}+\cdots )\\ & +{\epsilon }^{m_w-1}(\frac{\mathrm{\partial }{\overline{w}}_0}{\mathrm{\partial }Z}+\epsilon \frac{\mathrm{\partial }{\overline{w}}_1}{\mathrm{\partial }Z}+\cdots )=0.\end{array}$$

(31)

If $m_w<1$ , then it must hold that $\mathrm{\partial }{\overline{w}}_0/\mathrm{\partial }Z=0$ , i.e. ${\overline{w}}_0$ does not depend on $Z$ . This together with the boundary condition ${\overline{w}}_0(x,y,0,t)=0$ implies that ${\overline{w}}_0\equiv 0$ , which contradicts the assumption that ${\overline{w}}_0$ is not identically equal to zero. If $m_w>1$ , then we have that

$$\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }y}=0,$$

which is an unrealistic requirement on the leading terms in the velocity field. Hence $m_w=1$ , i.e. ${\overline{w}}_e$ is of order $\mathcal{O}(\epsilon )$ . By letting $\epsilon$ go to zero in (31) we obtain that

$$\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\overline{w}}_0}{\mathrm{\partial }Z}=0.$$

(32)

The no-slip boundary condition on the upper surface together with the fact the the upper surface only moves in the $z$ -direction imply that at the upper surface

$${\overline{w}}_{\epsilon }=\epsilon {\overline{w}}_0+{\epsilon }^2\overline{w}{}_1+\cdots =\frac{\mathrm{\partial }{h}_{\epsilon }}{\mathrm{\partial }t}=\epsilon \frac{\mathrm{\partial }h}{\mathrm{\partial }t},$$

and at the lower surface ${\overline{w}}_0=0,$ i.e.

$$\begin{array}{r}{\overline{w}}_0(x,y,h,t)=\frac{\mathrm{\partial }h}{\mathrm{\partial }t}(x,y,t)\mathrm{and}{\overline{w}}_0(x,y,0,t)=0.\end{array}$$

(33)

It remains to determine $m_p$ and the constants $k$ and $l$ , which are included in the expression (16) for the viscosity. Let us start by inserting the expansions in the expression (16) for the viscosity and then consider the limit as the film thickness goes to zero, i.e. the limit

$$\underset{\epsilon \to 0}{lim}f(\overline{\alpha }{\epsilon }^{k+m_p}({\overline{p}}_0+\epsilon {\overline{p}}_1+\cdots ),\overline{\beta }{\epsilon }^l({\overline{\theta }}_0+\epsilon {\overline{\theta }}_1+\cdots )).$$

(34)

We observe that, in order to obtain a useful limit problem, i.e. a limit problem in which the viscosity depends on the pressure and the temperature, then it must hold that $k+m_p=0$ and $l=0$ . Indeed, if $k+m_p<0$ , then the viscosity in the limit problem will correspond to the maximal viscosity ${\mu }_{max}$ . If $k+m_p>0$ , then the viscosity in the limit problem will be the viscosity at zero pressure. Similar arguments lead to the conclusion that $l=0$ . Let us now insert the expansions for ${\overline{p}}_{\epsilon }$ and ${\overline{\theta }}_{\epsilon }$ into (21) and (22), remember that $k+m_p=0$ and $l=0$ and let $\epsilon$ go to zero we obtain

$$0=\underset{\epsilon \to 0}{lim}(-{\epsilon }^{2+m_p}\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}+\frac{\mathrm{\partial }}{\mathrm{\partial }Z}(f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }Z})),$$

(35)

$$0=\underset{\epsilon \to 0}{lim}(-{\epsilon }^{2+m_p}\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }y}+\frac{\mathrm{\partial }}{\mathrm{\partial }Z}(f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }Z})).$$

(36)

Let us now determine $m_p$ . Indeed, if $m_p>-2$ , then the first term in both (35) and (36) will disappear in the limit, i.e., the limit problem will not include the pressure. Such a limit problem is of little use, since it enforces undue restrictions on the dominating terms ${\overline{u}}_0$ and ${\overline{v}}_0$ in the velocity field. For example, if $\mu$ is constant, then the dominating terms can only represent Couette flow. If $m_p<-2$ , then the second term in both (35) and (36) will vanish in the limit, i.e. we need ${\mathrm{\partial }}_x{\overline{p}}_0={\mathrm{\partial }}_y{\overline{p}}_0=0$ , which is not realistic except when the surfaces are parallel. Hence, $m_p=-2$ , which implies that $k=2.$ From (35) and (36) it now follows that

$$0=-\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}+\frac{\mathrm{\partial }}{\mathrm{\partial }Z}(f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }Z}),$$

(37)

$$0=-\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }y}+\frac{\mathrm{\partial }}{\mathrm{\partial }Z}(f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }Z}).$$

(38)

If we multiply both sides of (23) by $\epsilon$ , insert the expansions (26)-(30), with $m_u=0$ , $m_v=0$ , $m_w=1$ and $m_p=-2$ , into (23), have in mind that $k+m_p=0$ and $l=0$ , and let $\epsilon$ go to zero, then we obtain

$$0=\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }Z}.$$

(39)

By inserting the expansions into the scaled auxiliary energy equation (24), multiplying both sides with ${\epsilon }^2$ and passing to the limit we obtain

$$0=k\frac{{\mathrm{\partial }}^2{\overline{\theta }}_0}{\mathrm{\partial }{Z}^2}+f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)({\left(\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }Z}\right)}^2+{\left(\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }Z}\right)}^2).$$

(40)

We have now derived a system of equations (32), (37), (38), (39) and (40) for the lowest order terms $({\overline{\mathbf{v}}}_0,{\overline{p}}_0,{\overline{\theta }}_0)$ in the expansions of $({\overline{\mathbf{v}}}_{\epsilon },{\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })$ , respectively. This system is a first order approximation of the generalized incompressible Navier–Stokes system and can be used as a simplified model governing incompressible thin film flow of fluids with pressure and temperature dependent viscosity.

A system of Reynolds type for fluids with pressure and temperature dependent viscosity

In this section, we will derive a closed system of equations for the pressure and temperature developing in thin-film flows of incompressible fluids. The system consists of a generalized form of the Reynolds equation and a type of energy equation.

A generalized Reynolds equation

From (39) it follows that ${\overline{p}}_0$ does not depend on $Z$ . This implies that if we integrate both sides in (37) and (38) with respect to $Z$ , then we obtain that

$$\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }Z}=\frac{1}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)}\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}Z+\frac{A(x,y,t)}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)},$$

(41)

$$\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }Z}=\frac{1}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)}\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }y}Z+\frac{B(x,y,t)}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)}.$$

(42)

Integrating both sides of (41) across the fluid film, i.e. over $[0,h(x,y,t)]$ yields

$$\begin{array}{rl}& {\overline{u}}_0(x,y,h,t)-{\overline{u}}_0(x,y,0,t)\\ & =\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}{\int }_0^h\frac{ZdZ}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)}+A{\int }_0^h\frac{dZ}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)}.\end{array}$$

(43)

Let us define the functions

$$\mathcal{A}(x,y,Z,t)={\int }_0^Z\frac{\zeta d\zeta }{f(\overline{\alpha }{\overline{p}}_0(x,y,t),\overline{\beta }{\overline{\theta }}_0(x,y,\zeta ,t))},$$

(44)

$$\mathcal{B}(x,y,Z,t)={\int }_0^Z\frac{d\zeta }{f(\overline{\alpha }{\overline{p}}_0(x,y,t),\overline{\beta }{\overline{\theta }}_0(x,y,\zeta ,t))}.$$

(45)

By using the notation (44) and (45) and taking the no-slip boundary conditions at the lower and upper surface (assumed to be moving with velocities ${\mathbf{V}}_l=(u_l,v_l,0)$ and ${\mathbf{V}}_u=(0,0,w_u)$ , respectively) into account, (43) can be written as

$$-u_l(x,y,t)=\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}\mathcal{A}(x,y,h(x,y,t),t)+A(x,y,t)\mathcal{B}(x,y,h(x,y,t),t).$$

Hence

$$A(x,y,t)=-\frac{u_l}{\mathcal{B}(x,y,h,t)}-\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}\frac{\mathcal{A}(x,y,h,t)}{\mathcal{B}(x,y,h,t)}.$$

(46)

If we integrate (41), with respect to $Z$ , over the interval $[0,\zeta ]$ , then we obtain that

$$\begin{array}{rl}& {\overline{u}}_0(x,y,\zeta ,t)-u_l=\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}\mathcal{A}(x,y,\zeta ,t)\\ & +A(x,y,t)\mathcal{B}(x,y,\zeta ,t)\\ & =\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}(\mathcal{A}(x,y,\zeta ,t)-\frac{\mathcal{A}(x,y,h,t)\mathcal{B}(x,y,\zeta ,t)}{\mathcal{B}(x,y,h,t)})\\ & -\frac{u_l\mathcal{B}(x,y,\zeta ,t)}{\mathcal{B}(x,y,h,t)}.\end{array}$$

Thus (with $\zeta$ replaced by $Z$ )

$${\overline{u}}_0(x,y,Z,t)=\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}(\mathcal{A}(x,y,Z,t)-\frac{\mathcal{A}(x,y,h,t)\mathcal{B}(x,y,Z,t)}{\mathcal{B}(x,y,h,t)})+u_l(1-\frac{\mathcal{B}(x,y,Z,t)}{\mathcal{B}(x,y,h,t)}).$$

In the same way we obtain that

$${\overline{v}}_0(x,y,Z,t)=\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }y}(\mathcal{A}(x,y,Z,t)-\frac{\mathcal{A}(x,y,h,t)\mathcal{B}(x,y,Z,t)}{\mathcal{B}(x,y,h,t)})+v_l(1-\frac{\mathcal{B}(x,y,Z,t)}{\mathcal{B}(x,y,h,t)}).$$

The conservation of mass (32) yields

$$\frac{\mathrm{\partial }{\overline{w}}_0}{\mathrm{\partial }Z}=-\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }x}-\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }y}.$$

By integrating with respect to the $Z$ variable over the interval $[0,h]$ we obtain

$${\overline{w}}_0(x,y,h,t)-{\overline{w}}_0(x,y,0,t)=-{\int }_0^h\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }x}dZ-{\int }_0^h\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }y}dZ.$$

(47)

Let us recall the Leibniz integral rule

$$\frac{d}{dx}{\int }_{a(x)}^{b(x)}f(x,t)dt=f(x,b(x))\frac{db}{dx}(x)-f(x,a(x))\frac{da}{dx}(x)+{\int }_{a(x)}^{b(x)}\frac{\mathrm{\partial }f}{\mathrm{\partial }x}(x,t)dt.$$

By applying the Leibniz integral rule and taking the no-slip boundary conditions at the surfaces into account we can write (47) as

$${\overline{w}}_0(x,y,h,t)-{\overline{w}}_0(x,y,0,t)=-\frac{\mathrm{\partial }}{\mathrm{\partial }x}{\int }_0^h{\overline{u}}_{0}dZ-\frac{\mathrm{\partial }}{\mathrm{\partial }y}{\int }_0^h{\overline{v}}_{0}dZ.$$

(48)

Let us now rewrite the first integral on the right hand side of (48), i.e.,

$$\begin{array}{rl}{\int }_0^h{\overline{u}}_{0}dZ=& \frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}({\int }_0^h\mathcal{A}(x,y,Z,t)dZ\\ & -\frac{\mathcal{A}(x,y,h,t)}{\mathcal{B}(x,y,h,t)}{\int }_0^h\mathcal{B}(x,y,Z,t)dZ)\\ & +u_l(h-\frac{1}{\mathcal{B}(x,y,h,t)}{\int }_0^h\mathcal{B}(x,y,Z,t)dZ).\end{array}$$

(49)

Let us define

$$\mathcal{C}(x,y,t)={\int }_0^h\mathcal{A}(x,y,Z,t)dZ-\frac{\mathcal{A}(x,y,h,t)}{\mathcal{B}(x,y,h,t)}{\int }_0^h\mathcal{B}(x,y,Z,t)dZ,$$

(50)

$$\begin{array}{r}\mathcal{D}(x,y,t)=h-\frac{1}{\mathcal{B}(x,y,h,t)}{\int }_0^h\mathcal{B}(x,y,Z,t)dZ.\end{array}$$

(51)

This together with (49) gives

$${\int }_0^h{\overline{u}}_{0}dZ=\mathcal{C}(x,y,t)\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}+u_l\mathcal{D}(x,y,t).$$

(52)

In the same way we obtain

$${\int }_0^h{\overline{v}}_{0}dZ=\mathcal{C}(x,y,t)\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }y}+v_l\mathcal{D}(x,y,t).$$

(53)

The left hand side of (47) can in our case be rewritten with help of the boundary conditions (33) for ${\overline{w}}_0$ . The right hand side of (47) can be rewritten by using (52) and (53). After doing this we obtain

$$\frac{\mathrm{\partial }h}{\mathrm{\partial }t}=-\frac{\mathrm{\partial }}{\mathrm{\partial }x}(\mathcal{C}(x,y,t)\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}+u_l\mathcal{D}(x,y,t))-\frac{\mathrm{\partial }}{\mathrm{\partial }y}(\mathcal{C}(x,y,t)\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }y}+v_l\mathcal{D}(x,y,t))$$

or

$$\frac{\mathrm{\partial }h}{\mathrm{\partial }t}=-\mathrm{div}(\mathcal{C}(x,y,t)\mathrm{\nabla }{\overline{p}}_0+\mathcal{D}(x,y,t){\mathbf{U}}_l),$$

(54)

where ${\mathbf{U}}_l=u_l\mathbf{i}+v_l\mathbf{j}$ .

We note that in the stationary case, i.e. when $\mathrm{\partial }h/\mathrm{\partial }t=0$ , the generalized Reynolds equation (54) is the same that was derived from the Stokes equation in Ciuperca et al. ⁹

Remark.

To help the reader to compare the present result with the result in Dowson, ¹⁰ we just point out that by integration by parts we have

$$\begin{array}{rl}\mathcal{C}(x,y,t)& ={[Z\mathcal{A}(x,y,Z,t)]}_{Z=0}^{Z=h}\\ & -{\int }_0^{h}Z\frac{\mathrm{\partial }}{\mathrm{\partial }Z}\mathcal{A}(x,y,Z,t)dZ\\ & -\frac{\mathcal{A}(x,y,h,t)}{\mathcal{B}(x,y,h,t)}({[Z\mathcal{B}(x,y,Z,t)]}_{Z=0}^{Z=h}\\ & -{\int }_0^{h}Z\frac{\mathrm{\partial }}{\mathrm{\partial }Z}\mathcal{B}(x,y,Z,t)dZ)\\ & =-{\int }_0^{h}Z\frac{ZdZ}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)}\\ & +\frac{\mathcal{A}(x,y,h,t)}{\mathcal{B}(x,y,h,t)}{\int }_0^h\frac{ZdZ}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)}\\ & =-{\int }_0^h\frac{Z^{2}dZ}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)}+\frac{{\mathcal{A}}^2(x,y,h,t)}{\mathcal{B}(x,y,h,t)}\end{array}$$

and

$$\begin{array}{rl}& \mathcal{D}(x,y,t)=h(x,y,t)\\ & -\frac{1}{\mathcal{B}(x,y,h,t)}({[Z\mathcal{B}(x,y,Z,t)]}_{Z=0}^{Z=h}-{\int }_0^{h}Z\frac{\mathrm{\partial }}{\mathrm{\partial }z}\mathcal{B}(x,y,Z,t)dZ)\\ & =\frac{1}{\mathcal{B}(x,y,h,t)}{\int }_0^h\frac{ZdZ}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)}\\ & =\frac{\mathcal{A}(x,y,h,t)}{\mathcal{B}(x,y,h,t)}.\end{array}$$

The energy equation associated with the generalized Reynolds equation

Equation (40) may be written as

$$-k\frac{{\mathrm{\partial }}^2{\theta }_0}{\mathrm{\partial }{Z}^2}=f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0){\left|\frac{\mathrm{\partial }}{\mathrm{\partial }Z}{\overline{\mathbf{U}}}_0\right|}^2,$$

(55)

where ${\overline{\mathbf{U}}}_0={\overline{u}}_0\mathbf{i}+{\overline{v}}_0\mathbf{j}$ . By inserting (46) into (41) and (42) we obtain

$$\frac{\mathrm{\partial }{\overline{u}}_0}{\mathrm{\partial }Z}=\frac{(\mathcal{B}(x,y,h,t)Z-\mathcal{A}(x,y,h,t))\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}-u_l}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)\mathcal{B}(x,y,h,t)},$$

(56)

$$\frac{\mathrm{\partial }{\overline{v}}_0}{\mathrm{\partial }Z}=\frac{(\mathcal{B}(x,y,h,t)Z-\mathcal{A}(x,y,h,t))\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }y}-v_l}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0)\mathcal{B}(x,y,h,t)}.$$

(57)

From (56) and (57) it follows that (55) may be written as

$$-k\frac{{\mathrm{\partial }}^2{\theta }_0}{\mathrm{\partial }{Z}^2}=\frac{{|(\mathcal{B}(x,y,h,t)Z-\mathcal{A}(x,y,h,t))\mathrm{\nabla }{\overline{p}}_0-{\mathbf{U}}_l|}^2}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0){\mathcal{B}}^2(x,y,h,t)}.$$

(58)

The main result

Let us just sum up the main results in this section. Indeed, in Subsection “A generalized Reynolds equation” we derived the generalized equation (54) and in Subsection “The energy equation associated with the generalized Reynolds equation” we derived the equation (58) which is a special form of the energy equation. These two equation forms a system of Reynolds type. More precisely, we have derived the following model of incompressible thin film flow for fluids with pressure and temperature dependent viscosity:

$$\begin{array}{r}\frac{\mathrm{\partial }h}{\mathrm{\partial }t}=-\mathrm{div}(\mathcal{C}(x,y,t)\mathrm{\nabla }{\overline{p}}_0+\mathcal{D}(x,y,t){\mathbf{U}}_l)\mathrm{in}\omega ,\end{array}$$

(59)

$$\begin{array}{rl}& -k\frac{{\mathrm{\partial }}^2{\theta }_0}{\mathrm{\partial }{Z}^2}=\frac{{|(\mathcal{B}(x,y,h,t)Z-\mathcal{A}(x,y,h,t))\mathrm{\nabla }{\overline{p}}_0-{\mathbf{U}}_l|}^2}{f(\overline{\alpha }{\overline{p}}_0,\overline{\beta }{\overline{\theta }}_0){\mathcal{B}}^2(x,y,h,t)}\\ & & \mathrm{in}Q,\end{array}$$

(60 )

where $\mathcal{A},\mathcal{B},\mathcal{C}$ and $\mathcal{D}$ are defined as in (44), (45), (50) and (51), respectively.

Particular cases

Let us consider two examples:

Example 1. (The classical Reynolds equation)

If $\mu$ is constant, then

$$\mathcal{A}=\mathcal{A}(Z)={\int }_0^Z\frac{\xi }{\mu }d\xi =\frac{Z^2}{2\mu },\mathcal{B}=\mathcal{B}(Z)={\int }_0^Z\frac{1}{\mu }d\xi =\frac{Z}{\mu }.$$

We obtain

$$\begin{array}{rl}\mathcal{C}(x,y,t)& ={\int }_0^h\frac{Z^2}{2\mu }dZ-\frac{h}{2}{\int }_0^h\frac{Z}{\mu }dZ\\ & =\frac{h^3}{6\mu }-\frac{h^3}{4\mu }=-\frac{h^3}{12\mu },\\ \mathcal{D}(x,y,t)& =h-\frac{\mu }{h}{\int }_0^h\frac{Z}{\mu }dZ=h-\frac{\mu }{h}\frac{h^2}{2\mu }=\frac{h}{2}.\end{array}$$

This leads to the decoupled system

$$\{\begin{array}{cc}{\displaystyle \frac{\mathrm{\partial }h}{\mathrm{\partial }t}=\mathrm{div}(\frac{h^3}{12\mu }\mathrm{\nabla }{\overline{p}}_0-\frac{h}{2}{\mathbf{U}}_l)}& \mathrm{in}\omega ,\\ {\displaystyle -k\frac{{\mathrm{\partial }}^2{\theta }_0}{\mathrm{\partial }{Z}^2}=\frac{{|(\frac{hZ}{\mu }-\frac{h^2}{2\mu })\mathrm{\nabla }{\overline{p}}_0-{\mathbf{U}}_l|}^2}{\frac{h^2}{\mu }}}& \mathrm{in}Q.\end{array}$$

Example 2. (The generalized equation derived in Almqvist et al.^2,3)

If $\mu$ only depends on the pressure, then, by taking into account that ${\overline{p}}_0$ does not depend on $Z,$ it follows that

$$\begin{array}{rl}\mathcal{A}(x,y,Z,t)& ={\int }_0^Z\frac{\xi }{f\left(\overline{\alpha }{\overline{p}}_0\right)}d\xi =\frac{Z^2}{2f\left(\overline{\alpha }{\overline{p}}_0\right)},\\ \mathcal{B}(x,y,Z,t)& ={\int }_0^Z\frac{1}{f\left(\overline{\alpha }{\overline{p}}_0\right)}d\xi =\frac{Z}{f\left(\overline{\alpha }{\overline{p}}_0\right)},\end{array}$$

which implies that

$$\begin{array}{rl}\mathcal{C}(x,y,t)=& {\int }_0^h\frac{Z^2}{2f\left(\overline{\alpha }{\overline{p}}_0\right)}dZ-\frac{h}{2}{\int }_0^h\frac{Z}{f\left(\overline{\alpha }{\overline{p}}_0\right)}dZ\\ & =-\frac{h^3}{12f\left(\overline{\alpha }{\overline{p}}_0\right)},\\ \mathcal{D}(x,y,t)=& h-\frac{f\left(\overline{\alpha }{\overline{p}}_0\right)}{h}{\int }_0^h\frac{Z}{f\left(\overline{\alpha }{\overline{p}}_0\right)}dZ=\frac{h}{2}.\end{array}$$

This leads to the decoupled system

$$\{\begin{array}{cc}{\displaystyle \frac{\mathrm{\partial }h}{\mathrm{\partial }t}=\mathrm{div}(\frac{h^3}{12f\left(\overline{\alpha }{\overline{p}}_0\right)}\mathrm{\nabla }{\overline{p}}_0-\frac{h}{2}{\mathbf{U}}_l)}& \mathrm{in}\omega ,\\ {\displaystyle -k\frac{{\mathrm{\partial }}^2{\theta }_0}{\mathrm{\partial }{Z}^2}=\frac{{|(\frac{hZ}{f\left(\overline{\alpha }{\overline{p}}_0\right)}-\frac{h^2}{2f\left(\overline{\alpha }{\overline{p}}_0\right)})\mathrm{\nabla }{\overline{p}}_0-{\mathbf{U}}_l|}^2}{\frac{h^2}{f\left(\overline{\alpha }{\overline{p}}_0\right)}}}& \mathrm{in}Q.\end{array}$$

Remark: In the above examples the viscosity does not depend on the temperature and it possible to further simplify the system (59)–(60) by using that the pressure does not depend on Z. If the viscosity depends on the temperature, then it is not possible to do further simplifications since the temperature depends on $Z$ .

Asymptotic analysis of the boundary conditions

In this subsection we discuss the boundary conditions that should be imposed in the auxiliary system (17)–(19) and study the asymptotic behavior of these boundary conditions.

Let us start by considering the boundary condition at ${\mathrm{\Gamma }}_{\mathit{T}}$ , i.e. the part of lateral boundary ${\mathrm{\Gamma }}_w$ where the normal stress ${\mathit{T}}_{\epsilon }\mathbf{n}$ is prescribed. Indeed, for the physical problem we have that

$$(-p_{\epsilon }\mathbf{I}+2\mu (p_{\epsilon },{\theta }_{\epsilon })\mathbf{D}({\mathbf{v}}_{\epsilon }))\mathbf{n}=-p_b\mathbf{n}\mathrm{on}{\mathrm{\Gamma }}_{\mathit{T}}.$$

(61)

Let us impose the same boundary conditions in the auxiliary system (17)–(19). After transformation to a fixed domain, see Subsection “Transformation of the $\epsilon$ -dependent domain $Q_{\epsilon }$ into a fixed domain”, and multiplying with ${\epsilon }^2$ we obtain that on ${\mathrm{\Gamma }}_{\mathit{T}}$

$$\begin{array}{rl}& {\epsilon }^2{\overline{p}}_{\epsilon }{n}_1+{\epsilon }^{2}2{\overline{\mu }}_{\epsilon }(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }x}{n}_1+\frac{1}{2}(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x})n_2\\ & +\frac{1}{2}(\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }Z}+\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x})n_3)\\ & =-{\epsilon }^2{p}_b{n}_1,\\ & {\epsilon }^2{\overline{p}}_{\epsilon }{n}_2+{\epsilon }^{2}2{\overline{\mu }}_{\epsilon }(\frac{1}{2}(\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y})n_1+\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }y}{n}_2\\ & +\frac{1}{2}(\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }Z}+\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y})n_3)\\ & =-{\epsilon }^2{p}_b{n}_2,\\ & {\epsilon }^2{\overline{p}}_{\epsilon }{n}_3+{\epsilon }^{2}2{\overline{\mu }}_{\epsilon }(\frac{1}{2}(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x}+\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }Z})n_1\\ & +\frac{1}{2}(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y}+\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }Z})n_2+\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }Z}{n}_3)\\ & =-{\epsilon }^2{p}_b{n}_3.\end{array}$$

If we let $\epsilon$ go to zero, then we obtain that ${\overline{p}}_0=0$ on ${\mathrm{\Gamma }}_{\mathit{T}}$ whatever the constant $p_b$ is. This is not realistic! In order to get a limit problem which reflects the normal stress boundary condition we must scale the right hand side. Indeed, the appropriate boundary condition for the auxiliary system (17)–(19) is

$$(-{\overline{p}}_{\epsilon }\mathbf{I}+2\overline{\mu }({\overline{p}}_{\epsilon },{\overline{\theta }}_{\epsilon })\mathbf{D}({\overline{\mathbf{v}}}_{\epsilon }))\mathbf{n}=-\frac{{\overline{\epsilon }}^2}{{\epsilon }^2}{p}_b\mathbf{n}\mathrm{on}{\mathrm{\Gamma }}_{\mathit{T}},$$

(62)

which after multiplying by ${\epsilon }^2$ reads in component form

$$\begin{array}{rl}& {\epsilon }^2{\overline{p}}_{\epsilon }{n}_1+{\epsilon }^{2}2{\overline{\mu }}_{\epsilon }(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }x}{n}_1+\frac{1}{2}(\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x})n_2\\ & +\frac{1}{2}(\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }Z}+\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x})n_3)\\ & =-{\overline{\epsilon }}^2{p}_b{n}_1,\\ & {\epsilon }^2{\overline{p}}_{\epsilon }{n}_2+{\epsilon }^{2}2{\overline{\mu }}_{\epsilon }(\frac{1}{2}(\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }y})n_1\\ & +\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }y}{n}_2+\frac{1}{2}(\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }Z}+\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y})n_3)\\ & =-{\overline{\epsilon }}^2{p}_b{n}_2,\\ & {\epsilon }^2{\overline{p}}_{\epsilon }{n}_3+{\epsilon }^{2}2{\overline{\mu }}_{\epsilon }(\frac{1}{2}(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }x}+\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{u}}_{\epsilon }}{\mathrm{\partial }Z})n_1\\ & +\frac{1}{2}(\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }y}+\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{v}}_{\epsilon }}{\mathrm{\partial }Z})n_2+\frac{1}{\epsilon }\frac{\mathrm{\partial }{\overline{w}}_{\epsilon }}{\mathrm{\partial }Z}{n}_3)\\ & =-{\overline{\epsilon }}^2{p}_b{n}_3.\end{array}$$

In the limit, as $\epsilon$ tends to zero, we obtain that

$${\overline{p}}_0=-{\overline{\epsilon }}^2{p}_b\mathrm{on}{\mathrm{\Gamma }}_{\mathit{T}}.$$

This means that ${\mathit{p}}_{\overline{\epsilon }}={\overline{p}}_{\overline{\epsilon }}\approx \frac{1}{{\overline{\epsilon }}^2}{\overline{p}}_0=-p_b$ on ${\mathrm{\Gamma }}_{\mathit{T}}$ . Note that for $\epsilon =\overline{\epsilon }$ the boundary condition (62) for the auxiliary problem and the boundary condition (61) for the physical problem are the same.

Let us consider the boundary condition on ${\mathrm{\Gamma }}_{\mathbf{v}}$ , i.e. the part of lateral boundary ${\mathrm{\Gamma }}_w$ where the the velocity ${\mathbf{v}}_{\epsilon }$ is given. In the auxiliary problem we impose the boundary condition $({\overline{u}}_{\epsilon },{\overline{v}}_{\epsilon },{\overline{w}}_{\epsilon })=(u_b,v_b,\epsilon w_b)$ . The third component is scaled to be in balance with ${\overline{w}}_{\epsilon }$ . Integrating the first two components in both sides with respect to $Z$ over $[0,h]$ and using (52) and (53) yields

$$\begin{array}{rl}\mathcal{C}(x,y,t)\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }x}+u_l\mathcal{D}(x,y,t)& ={\int }_0^h{u}_{b}dZ=u_{b}h,\\ \mathcal{C}(x,y,t)\frac{\mathrm{\partial }{\overline{p}}_0}{\mathrm{\partial }y}+v_l\mathcal{D}(x,y,t)& ={\int }_0^h{v}_{b}dZ=v_{b}h.\end{array}$$

This implies that

$$\mathrm{\nabla }{\overline{p}}_0\cdot \mathbf{n}=(\frac{h(x,y,t)}{\mathcal{C}(x,y,t)}{\mathbf{U}}_b-\frac{\mathcal{D}(x,y,t)}{\mathcal{C}(x,y,t)}{\mathbf{U}}_l)\cdot \mathbf{n}\mathrm{on}{\mathrm{\Gamma }}_{\mathbf{v}},$$

(63)

where ${\mathbf{U}}_l=u_l\mathbf{i}+v_l\mathbf{j}$ and ${\mathbf{U}}_b=u_b\mathbf{i}+v_b\mathbf{j}$ . Thus a Dirichlet boundary condition for the velocity in the generalized Navier–Stokes equation yields a Neumann boundary condition for the pressure in the lower-dimensional equation (59) for the pressure.

Example 3.

If $\mu$ is constant, then $\mathcal{C}=-h^3/(12\mu )$ and $\mathcal{D}=h/2$ according to Example 1. From (63) it follows that

$$\mathrm{\nabla }{\overline{p}}_0\cdot \mathbf{n}=(-\frac{12\mu }{h^2}{\mathbf{U}}_b+\frac{6\mu }{h^2}{\mathbf{U}}_l)\cdot \mathbf{n}\mathrm{on}{\mathrm{\Gamma }}_{\mathbf{v}}.$$

Example 4.

If $\mu$ depends on the pressure, then $\mathcal{C}=-h^3/(12f(\overline{\alpha }{\overline{p}}_0))$ and $\mathcal{D}=h/2$ according to Example 2. From (63) it follows that

$$\mathrm{\nabla }{\overline{p}}_0\cdot \mathbf{n}=(-\frac{12f\left(\overline{\alpha }{\overline{p}}_0\right)}{h^2}{\mathbf{U}}_b+\frac{6f\left(\overline{\alpha }{\overline{p}}_0\right)}{h^2}{\mathbf{U}}_l)\cdot \mathbf{n}\mathrm{on}{\mathrm{\Gamma }}_{\mathbf{v}}.$$

Conclusions and concluding remarks

In this paper we have derived a model of incompressible thin film flow, where the viscosity of the fluid depends on pressure and temperature. It, thereby, constitutes an addition, i.e. Part C, to the previously published two-part paper.^2,3 We have also pointed out shortcomings in the derivations of some of the models that are widely used to model such flows. The main conclusions and results are:

The Navier–Stokes model of an incompressible fluid can not be used since it does not allow the viscosity to depend on the pressure. Hence, to derive equations governing the motion of incompressible fluids with pressure and temperature dependent viscosity one must start by deriving a constitutive relation where the material moduli may depend on the pressure and the temperature. This is done by applying the implicit constitutive theory.

After inserting the derived constitutive relation into the equations representing balance of linear momentum and balance of energy we observe that we must scale the viscosity-pressure coefficient $\alpha$ in a precise way to obtain a simplified system of equations which reflects that the viscosity depends on the pressure. Thus it is not possible to derive the limit problem just by letting the film thickness go to zero as in the classical case with constant viscosity.

From the simplified system, we derived a system of equations for the pressure and the temperature. The system consists of a generalized form of the Reynolds equation (59) and an energy equation (60). The generalized form of the Reynolds equation is a lower-dimensional model, i.e. the domain is two-dimensional. The domain associated with the equation which is related to the energy is three-dimensional. This is because the temperature varies across the film.

A Dirichlet boundary condition for the velocity in the generalized Navier–Stokes equation yields a Neumann boundary condition for the pressure in the lower-dimensional equation (59). Moreover, a normal stress condition in the generalized Navier–Stokes equation yields a Dirichlet condition for the pressure in the lower-dimensional equation (59).

This work presents, by no means, the final solution to thin-film flow modeling. There are various points of improvement that could be addressed in future work, drawing from the results presented herein. The list can be made nearly infinitely long but, we can only mention a few. For instance, it is necessary to rethink and work on adapting suitable boundary conditions (when applying the present modeling approach) to a particular application. Relaxing the assumption of fluid incompressibility (in a physically sound way) is another one. Rigorously including a rheology which is more complex than the constitutive relationship given by (8) is a third. It is also essential that future research is devoted to possibly including the fluid’s phase change from liquid to solid and the associated plastic flow of the solid within the contact region. Another important aspect is to develop numerical methods for solving the lower-dimensional problems, associated with various application. However, it is well-known that even in the simple isothermal situation, it is a very challenging task to get the numerical scheme to converge.