## 4 Numerical Results
### 2D Droplet Impinging on a Thin Liquid Film  
The DG-VDI solver is employed for simulating the impact of a two-dimensional droplet on a thin liquid layer. Since we neglect the effect of the surface tension, only the early stages of the interaction between the droplet and the film are simulated. Water (subscript \(w\)) and air (subscript \(a\)) are subject to the gravitational field. Reynolds and Froude numbers based on water properties, droplet diameter, and droplet initial vertical velocity are set to \(\mathrm{Re}=1000\) and \(\mathrm{Fr}=20\), respectively. As viscosity is different within the two fluids, a suitable function \(\mu(\rho)\) has to be introduced. To this end, we rely on the following definition of the mixture of kinematic viscosity  
\[\
u=\
u_{w}\varphi_{w}+\
u_{a}\varphi_{a},\]  
where \(\varphi\) is the mass fraction. Taking advantage of the definition \(\mu=\rho\
u\) and the property \(\varphi_{w}+\varphi_{a}=1\), we are able to obtain  
\[\mu(\rho,\varphi_{a})=\rho\
u_{w}+\rho\varphi_{a}(\
u_{a}-\
u_{w})=\frac{\rho }{\rho_{w}}\mu_{w}+\frac{\rho}{\rho_{w}}\varphi_{a}\big{(}r_{\rho}\mu_{a}-\mu _{w}\big{)}. \tag{29}\]  
If we now adopt the model of a simple mixture, we are able to write  
\[\frac{1}{\rho}=\frac{1}{\rho_{w}}\varphi_{w}+\frac{1}{\rho_{a}}\varphi_{a},\]  
and, accordingly, we have  
\[\varphi_{a}(\rho)=\frac{\frac{\rho_{w}}{\rho}-1}{r_{\rho}-1}. \tag{30}\]  
Substituting (30) in (29), after some trivial algebraic manipulations, we finally obtain  
\[\mu(\rho)=\mu_{w}\bigg{(}\frac{\rho}{\rho_{w}}\frac{r_{\rho}}{r_{\mu}}\frac{r_ {\mu}-1}{r_{\rho}-1}+\frac{1}{r_{\mu}}\frac{r_{\rho}-r_{\mu}}{r_{\rho}-1}\bigg{)}, \tag{31}\]  
where \(r_{\rho}=\rho_{w}/\rho_{a}\) and \(r_{\mu}=\mu_{w}/\mu_{a}\) are the density and dynamic viscosity ratios, respectively. The expression (31) defines a linear dependence between the dynamic viscosity and the density and ensures that \(\mu>0\) provided that \(\rho\in[\rho_{a},\rho_{w}]\). However, at the discrete level, such density limits are not verified with the current approach, see Section 3.3. Nevertheless, assuming \(\rho>0\) (which is ensured by both exponential (12) and hyperbolic tangent (13) transformations), it is easy to conclude that \(\mu>0\) holds if  
\[1\leq r_{\mu}\leq r_{\rho}.\]  
Since \(r_{\rho}=800\) and \(r_{\mu}=50\), the dynamic viscosity computed via (31) is always positive. Given the square domain \(\Omega=(0,4)\times(0,4)\), the non-dimensional initial field is defined as  
\[p =1,\] \[\mathbf{u} =\left\{ \begin{array}{ll}-j&|\mathbf{x}-\mathbf{x}_{c}|\leq 0.5\\ \mathbf{0}&\text{otherwise}\end{array} \right.,\] \[\rho =\left\{ \begin{array}{ll}1&|\mathbf{x}-\mathbf{x}_{c}|\leq 0.5\text{ or }y \leq\delta\\ r_{\rho}^{-1}&\text{otherwise}\end{array} \right.,\]  
being \(\mathbf{x}_{c}=1.15j\) the droplet center and \(\delta=0.2\) the film thickness. In practice, the smoother function  
\[\text{erf}(\texttt{ad}),\]  
where d is the distance of any point \(\mathbf{x}\in\Omega\) from the (velocity or density) discontinuity, while \(\alpha\in\mathbb{R}^{+}\) is a user-defined scaling parameter and is applied to spread jumps intosmooth variations of the chosen variable and avoid discontinuities within density and velocity initial fields. Notice that the smoother function is applied in the radial direction for the droplet and along the vertical direction for the layer. Top and lateral domain boundaries are treated with symmetric conditions and the bottom boundary with a no-slip wall condition. The simulation is performed for three different computational setups of increasingly high accuracy, denoted as A, B, and C. For each setup, a different irregular Cartesian grid is considered. Such grids are characterized by a uniformly refined region at the bottom-left corner, i.e., for \(x\in(0,2)\times(0,1.8)\), with element dimension \(h=1/15\) (grid A), \(1/30\) (grid B) and \(1/60\) (grid C), and by a coarsening region in the remaining part of the domain \(\Omega\).