1. Governing equations

For a detailed derivation and description of the following overview, see Holey, H. et al., Tribology Letters 70 (2022).

Splitting Approach

The local mass and momentum balances are given by the partial differential equation:

\[ \frac{\partial \mathbf{q}}{\partial t} = - \mathbf{\nabla} \cdot \mathbb{F}(\mathbf{q}) \]

where

\[\begin{split} \mathbf{q} = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \end{pmatrix}, \quad \mathbb{F}(\mathbf{q}) = \begin{pmatrix} \rho u & \rho v & \rho w \\ \rho u^2 + p - \tau_{\mathrm{11}}(\mathbf{q}) & \rho u v - \tau_{\mathrm{12}}(\mathbf{q}) & \rho u w - \tau_{\mathrm{13}}(\mathbf{q}) \\ \rho v u - \tau_{\mathrm{21}}(\mathbf{q}) & \rho v^2 + p - \tau_{\mathrm{22}}(\mathbf{q}) & \rho v w - \tau_{\mathrm{23}}(\mathbf{q}) \\ \rho w u - \tau_{\mathrm{31}}(\mathbf{q}) & \rho w v - \tau_{\mathrm{32}}(\mathbf{q}) & \rho w^2 + p - \tau_{\mathrm{33}}(\mathbf{q}) \end{pmatrix}, \end{split}\]

where \(\rho\) denotes mass density, \(u\), \(v\), \(w\) are the Cartesian components of the velocity field \(\vec{u}\), \(\tau_{ij}\) are components of the viscous stress tensor \(\underline{\tau}\), and \(p\) denotes the fluid pressure.

In the following we treat the mass fluxes \(j_x\equiv\rho u\), \(j_y\equiv\rho v\), and \(j_z\equiv\rho w\) as independent variables, i.e. we work with conservative instead of primitive variables.

The divergence operates on the rows such that the continuum equation states (similarly for the momentum equations):

\[ \frac{\partial \rho}{\partial t} = - \frac{\partial j_x}{\partial x} - \frac{\partial j_y}{\partial y} - \frac{\partial j_z}{\partial z} \]

We now assume that the problem can be divided into a time-dependent macro problem that uses the height-average of the conserved quantities \(\overline{\mathbf{q}}(x,y,t)\) and a steady-state micro problem that resolves the values of the conserved quantities across the film thickness \(\delta \mathbf{q}(z)\).

\[\mathbf{q}(x,y,z,t) = \overline{\mathbf{q}}(x,y,t) + \delta \mathbf{q}(z)\]

Macro Problem

Derivation

To establish the macro problem, we eliminate the \(z\)-direction dependency by integrating the balances across the film thickness

\[ \int_{h_{\mathrm{1}}}^{h_{\mathrm{2}}} \left( \frac{\partial \mathbf{q}}{\partial t} - \mathbf{\nabla} \cdot \mathbb{F}(\mathbf{q}) \right) \, dz = 0. \]

Here, \(h_1\equiv h_1(x,y,t)\) and \(h_2\equiv h_2(x,y,t)\) describe the topographies of the lower and upper wall, respectively. Their difference is denoted with \(h \equiv h(x,y,t) = h_2(x,y,t) - h_1(x,y,t)\). We need to account for the spatial and temporal dependencies of the topography when evaluating the above integral. Let’s look at what we obtain by integrating the individual terms:

---

Integral of the time derivative term on the left side:

\[ \int_{h_{\mathrm{1}}}^{h_{\mathrm{2}}} \frac{\partial \mathbf{q}}{\partial t} \, dz = h \frac{\partial \overline{\mathbf{q}}}{\partial t} + \left( \overline{\mathbf{q}} - \mathbf{q} \big|_{h_{\mathrm{2}}} \right) \frac{d h_{\mathrm{2}}}{d t} - \left( \overline{\mathbf{q}} - \mathbf{q} \big|_{h_{\mathrm{1}}} \right) \frac{d h_{\mathrm{1}}}{d t} \]

Here, we obtain the time derivative of the height averaged conserved quantities \(\frac{\partial \overline{\mathbf{q}}}{\partial t}\) multiplied by \(h\) as well as some source terms.

---

Integral of the flux divergence on the right side (with \(\mathbf{f}_x\), \(\mathbf{f}_y\), \(\mathbf{f}_z\) denoting the columns of \(\mathbb{F}\) from left to right).

\[ \int_{h_{\mathrm{1}}}^{h_{\mathrm{2}}} \left( \mathbf{\nabla} \cdot \mathbb{F}(\mathbf{q}) \right) \, dz = \int_{h_{\mathrm{1}}}^{h_{\mathrm{2}}} \frac{\partial \mathbf{f}_x}{\partial x} \, dz + \int_{h_{\mathrm{1}}}^{h_{\mathrm{2}}} \frac{\partial \mathbf{f}_y}{\partial y} \, dz + \int_{h_{\mathrm{1}}}^{h_{\mathrm{2}}} \frac{\partial \mathbf{f}_z}{\partial z} \, dz \]

For \(x\) (and similarly for \(y\)) we obtain the following terms. Note that we obtain a height-averaged \(\overline{\mathbf{f}}_x\) (and similarly \(\overline{\mathbf{f}}_y\)).

\[ \int_{h_{\mathrm{1}}}^{h_{\mathrm{2}}} \frac{\partial \mathbf{f}_x}{\partial x} \, dz = h \frac{\partial \overline{\mathbf{f}}_x}{\partial x} + \left( \overline{\mathbf{f}}_x - \mathbf{f}_x \big|_{h_{\mathrm{2}}} \right) \frac{d h_{\mathrm{2}}}{d t} - \left( \overline{\mathbf{f}}_x - \mathbf{f}_x \big|_{h_{\mathrm{1}}} \right) \frac{d h_{\mathrm{1}}}{d t} \]

Here, we obtain the spatial derivative of the height averaged fluxes \(\frac{\partial \overline{\mathbf{f}_x}}{\partial x}\) and \(\frac{\partial \overline{\mathbf{f}_y}}{\partial y}\) as well as - again - some source terms.

In contrast, for the \(\mathbf{f}_z\) term, the integration direction is the same as that of the derivative. The term can be expressed by evaluating the function values of \(\mathbf{f}_z\) at the bottom and top:

\[ \int_{h_{\mathrm{1}}}^{h_{\mathrm{2}}} \frac{\partial \mathbf{f}_z}{\partial z} \, dz = \mathbf{f}_z \big|_{h_{\mathrm{2}}} - \mathbf{f}_z \big|_{h_{\mathrm{1}}} \]

These are only source terms.

---

Overall we obtain the following equation with all the source terms occurring in the previous equations collected in \(\mathbf{s}\). Note that we have reduced the problem from three to two dimensions.

\[ \frac{\partial \overline{\mathbf{q}}}{\partial t} = - \frac{\partial \overline{\mathbf{f}}_x}{\partial x} - \frac{\partial \overline{\mathbf{f}}_y}{\partial y} - \mathbf{s} \]

with four variables \(\overline{\rho}(x,y,t)\), \(\overline{j}_x(x,y,t)\), \(\overline{j}_y(x,y,t)\), and \(\overline{j}_z(x,y,t)\) and four equations, which define the macro problem to be solved.

Further unknowns (within \(\overline{\mathbf{f}}_x\), \(\overline{\mathbf{f}}_y\), \(\mathbf{s}\)) are \(\tau(\mathbf{q})\), \(p(\rho)\), and the \(\mathbf{q}\) and \(\mathbf{f}\) values at the bottom and top boundaries. As you might guess, this is the task for the micro problem.

Numerical Solution of the Macro Problem

For numerical solution, we use the MacCormack scheme. It is a two-step explicit predictor–corrector method for solving hyperbolic conservation laws and provides second-order accuracy in space and time even when including source terms as we need to do.

For illustration, let’s assume this simplified 1D problem:

\[ \frac{\partial \mathbf{q}}{\partial t} + \frac{\partial \mathbf{f}(\mathbf{q})}{\partial x} = 0 \]

where:

• \(\mathbf{q}\) — vector of conserved variables

• \(\mathbf{f}(\mathbf{q})\) — flux vector

In a first (predictor) step, we calculate an intermediate solution by forward difference

\[ \mathbf{q}_i^* = \mathbf{q}_i^n - \frac{\Delta t}{\Delta x} \left[ \mathbf{f}(\mathbf{q}_{i+1}^n) - \mathbf{f}(\mathbf{q}_i^n) \right] \]

followed by the corrector step, where we compute the backward difference based on the intermediate solution.

\[ \mathbf{q}_i^{n+1} = \frac{1}{2} \left[ \mathbf{q}_i^n + \mathbf{q}_i^* - \frac{\Delta t}{\Delta x} \left( \mathbf{f}(\mathbf{q}_i^*) - \mathbf{f}(\mathbf{q}_{i-1}^*) \right) \right] \]

Note that it is possible to switch the order of forward- and backward-difference between each timestep, obtaining a more symmetric solution scheme.

An important concept that comes with explicit schemes is the CFL (Courant–Friedrichs–Lewy) number, which determines the stability of the numerical method. The CFL number is defined as

\[ \text{CFL} = \frac{\max(|u|)\, \Delta t}{\Delta x} \]

and represents the ratio of the physical wave propagation speed to the numerical propagation speed. For stability in explicit schemes, the CFL condition must be satisfied:

\[ \text{CFL} \le 1 \]

Note that depending on the problem you might want to select a lower CFL number for improved convergence.

Micro Problem

Remember the general mass and momentum balances and the splitting approach:

\[ \frac{\partial \mathbf{q}}{\partial t} = - \mathbf{\nabla} \cdot \mathbb{F}(\mathbf{q}) \]

\[\mathbf{q}(x,y,z,t) = \overline{\mathbf{q}}(x,y,t) + \delta \mathbf{q}(z)\]

Now if we put the splitting approach into the general balance we get:

\[ \frac{\partial \overline{\mathbf{q}}}{\partial t} = - \frac{\partial \overline{\mathbf{f}}_x}{\partial x} - \frac{\partial \overline{\mathbf{f}}_y}{\partial y} - \frac{\partial \mathbf{f}_z}{\partial \delta \mathbf{q}} \frac{\partial \delta \mathbf{q}}{\partial z} \]

Comparing this to the resulting equation of the macro problem, we see that

\[ \frac{\partial \mathbf{f}_z}{\partial \delta \mathbf{q}} \frac{\partial \delta \mathbf{q}}{\partial z} = \mathbf{s} \]

meaning that this expression has to be a constant source term at a given \((x,y,t)\). We can evaluate and simplify the above equation to:

\[\begin{split} \frac{\partial}{\partial z} \begin{pmatrix} \tau_{\text{xz}} \\ \tau_{\text{yz}} \\ p + \tau_{\text{zz}} \\ \end{pmatrix} = \text{const} \end{split}\]

This information now allows us to determine the profile of the conserved quantities across the film and is the basis for the solution of the micro problem. In the following tutorial, we will determine the velocity profile across the film using these conditions (remember that \(\tau_{\text{xz}} \propto \frac{\partial u}{\partial z}\) and \(\tau_{\text{yz}} \propto \frac{\partial v}{\partial z}\)).